版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜

返回列表

当前主题已经存档。

gengbiaolu

铜虫 (正式写手)

应助: 0 (幼儿园)
金币: 8.7
散金: 90
红花: 2
帖子: 458
在线: 70.2小时
虫号: 762835
注册: 2009-05-04
专业: 原子和分子物理

[交流] 【求助】我现在有一个C++的程序，但是我看不懂，求一高人帮我将它转成Matlab程序。已有4人参与

我现在有一个C++的程序，但是我看不懂，求一高人帮我将它转成Matlab程序。
/***********************************************************
*                薛订鄂方程数值模拟                   *
*                2009.10.25                            *
***********************************************************/
#include //要使用cin或cout进行输入输出;数据流输入／输出
#include    //要使用sin、cos、squrt、fabs等其中的某一个函数
#include    //文件输入／输出
#include    //调用标准库函数rand()返回0-32767间的一个int型随机数;定义杂项函数及内存分配函数
#include
#include //格式控制符
#include    //包括一系列日期和时间处理函数

using namespace std;

#define    pi  3.1415926
/*----------*
|  常量定义  |
*-----------*/

const double A=0.02;             //振幅---自己根据情况定

const int L=60;                //格子长度---自己根据情况定
const double beita=1.0;          //特征非线性因子---自己根据情况定

const double C=1.0;
const double Lamada=1.6*C;       //Lamada与C的关系---自己根据情况定

const int n=L;                //Total Lattice site numbers

const double dt=1.0e-3;             //时间步长选取(不清楚是无量纲还是有量纲)---自己根据情况定

const double eps=1.0e-5;       //判断达到近似稳定条件--上下两层误差范围
/*------------*
|  子函数声明  |
*-------------*/
void initialization1(double *a,int n,double number0);

/*--------*
|  主函数  |
*---------*/

void main()
{//0
int i,j;
ofstream result("psi-norm1.txt"

;                   //psi-1范数随n变化输出文件

double *E;                                           //格子点处能量赋值---自己根据情况定
E=new double[n];
for(i=0;i {
// E=(4.0/(n-1))*i-2.0;
// E=1.0;
   E=beita+Lamada*cos(2.0*pi*(i-29)*0.5*(sqrt(5.0)+1));

}

double *a,*b;                                        //psi-实部和虚部变量分别声称
a=new double[n];                                     //指针动态申请
b=new double[n];
initialization1(a,n,0.006);                         //变量初始化-也为0时刻初始值---自己根据情况定
initialization1(b,n,0.006);

double *psi2;                                        //psi-2范数变量声称
psi2=new double[n];
initialization1(psi2,n,0.0);

double T=1.0e+10;                                     //时间限制
double time=0.0;                                     //开始计时
int num=0;                                           //推进步数

while(time {//1
num=num+1;
time=time+dt;

double *Da,*Db;                               //改进Euler法中求f(xn,yn)
Da=new double[n];
Db=new double[n];
      initialization1(Da,n,0.0);
      initialization1(Db,n,0.0);

      double *az,*bz;                               //改进Euler法中求y(n+1)-中间过渡量
az=new double[n];
bz=new double[n];
      initialization1(az,n,0.0);
      initialization1(bz,n,0.0);

double *Daz,*Dbz;                               //改进Euler法中求f(x(n+1),y(n+1))
Daz=new double[n];
Dbz=new double[n];
      initialization1(Daz,n,0.0);
      initialization1(Dbz,n,0.0);

for(i=1;i {//2
      Da=Da+E*b+beita*psi2*b+b[i+1]+b[i-1];       //文献1离散方式
   Db=Db-E*a-beita*psi2*a-a[i+1]-a[i-1];

//    Da=Da-E*b-C*(b[i+1]+b[i-1]);       //文献2离散方式
//    Db=Db+E*a+C*(a[i+1]+a[i-1]);

   az=a+dt*Da;
   bz=b+dt*Db;
}//2

// az[0]=2.0*A+bz[1];                               //文献1边界条件处理---自己根据情况定
      az[0]=bz[1];                                     //自己给定边界条件
bz[0]=-az[1];
az[n-1]=bz[n-2];
bz[n-1]=-az[n-2];

for(i=1;i {//2
   psi2=az*az+bz*bz;

         Daz=Daz+E*bz+beita*psi2*bz+bz[i+1]+bz[i-1];
      Dbz=Dbz-E*az-beita*psi2*az-az[i+1]-az[i-1];

//    Daz=Daz-E*bz-C*(bz[i+1]+bz[i-1]);       //文献2离散方式
//    Dbz=Dbz+E*az+C*(az[i+1]+az[i-1]);
}//2

double *aaz,*bbz;                                  //物理量更新时中间传值变量声称
aaz=new double[n];
bbz=new double[n];
      initialization1(aaz,n,0.0);
initialization1(bbz,n,0.0);

for(i=1;i {//2
aaz=a+0.5*dt*(Da+Daz);
         bbz=b+0.5*dt*(Db+Dbz);
}//2

   // aaz[0]=2.0*A+bbz[1];                               //文献1边界条件处理---自己根据情况定
      aaz[0]=bbz[1];                                     //自己给定边界条件
bbz[0]=-aaz[1];
aaz[n-1]=bbz[n-2];
bbz[n-1]=-aaz[n-2];

double error=0.0;                                  //求下一时刻与上一时刻物理量绝对误差
for(i=0;i {
error=error+sqrt((aaz-a)*(aaz-a)+(bbz-b)*(bbz-b));
}

   for(i=0;i {
a=aaz;
b=bbz;
}

      for(i=0;i {psi2=a*a+b*b;}

delete[] Da;delete[] Db;                            //指针删除
delete[] az;delete[] bz;
delete[] Daz;delete[] Dbz;
delete[] aaz;delete[] bbz;

cout<<"推进步数为"<
if(num%10000==0)                                        //每隔10000步输出|PSI|n与粒子位置关系
{//2
         result<<"zone f=point"<    for(i=0;i    {//3
               result<)<    }//3                                           //程序运行结束
}//2

if(error<=eps||num%20000000==0)                                        //达到近似稳定状态或运行到2000万步后结果输出、程序运行结束
{//2
         result<<"zone f=point"<    for(i=0;i    {//3
               result<<    }//3
   cout<<"达到近似稳定后需要推进步数为-------"<    cout<<"达到近似稳定后需要时间为***********"<          break;                                              //程序运行结束
}//2

}//1

}//0

/*------------*
|  子函数定义  |
*-------------*/

//赋初始值
void initialization1(double *a,int n,double number0)
{
int i;
for(i=0;i *(a+i)=number0;
}

回复此楼

» 猜你喜欢

论文终于录用啦！满足毕业条件了已经有27人回复
假如你的研究生提出不合理要求已经有3人回复
所感已经有3人回复
要不要辞职读博？已经有7人回复
不自信的我已经有11人回复
北核录用已经有3人回复
实验室接单子已经有3人回复
磺酰氟产物，毕不了业了！已经有8人回复
求助:我三月中下旬出站，青基依托单位怎么办？已经有10人回复
26申博（荧光探针方向，有机合成）已经有4人回复

1楼 2010-04-05 23:59:31

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

mopsite

木虫 (著名写手)

应助: 2 (幼儿园)
金币: 8677.7
散金: 424
红花: 2
帖子: 2544
在线: 615小时
虫号: 172296
注册: 2006-01-19
专业: 传热传质学

自己看下 c语言书，很好懂的。

赞一下

回复此楼

2楼2010-04-06 09:46:27

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

gengbiaolu

铜虫 (正式写手)

应助: 0 (幼儿园)
金币: 8.7
散金: 90
红花: 2
帖子: 458
在线: 70.2小时
虫号: 762835
注册: 2009-05-04
专业: 原子和分子物理

kuhailangyu:那你连程序的意思都不懂的话，转成M文件后还是出错怎么办？接着再求助。。。 2010-04-06 18:49

我从未学过C语言，完全不懂呀，帮帮我吧！

赞一下(1人)

回复此楼

3楼2010-04-06 17:34:14

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者

仿真EPI: 5
应助: 0 (幼儿园)
贵宾: 0.208
金币: 467.3
红花: 3
帖子: 416
在线: 89.3小时
虫号: 943880
注册: 2010-01-17
性别: GG
专业: 凝聚态物性 II ：电子结构

★ ★
adu886886(金币+2):辛苦了，谢谢交流 2010-04-12 08:11
gengbiaolu(金币+15):辛苦了，谢谢！ 2010-04-14 08:53

代码如下：

引用回帖:

%/*--------*
%|  主函数  |
%*---------*/
function SchrodingerFunction

%%PSL@CSU
%%QQ:547423688
%%Email:anyuezhiji@qq.com
%%Edit @ 2010.4.12

%/*----------*
%|  常量定义  |
%*-----------*/

A=0.02;             %振幅---自己根据情况定
L=60;                %格子长度---自己根据情况定
beita=1.0;          %特征非线性因子---自己根据情况定
C=1.0;
Lamada=1.6*C;       %Lamada与C的关系---自己根据情况定
n=L;                %Total Lattice site numbers
dt=1.0e-3;          %时间步长选取(不清楚是无量纲还是有量纲)---自己根据情况定
MyEps=1.0e-5;       %判断达到近似稳定条件--上下两层误差范围

int i j;
fidout=fopen('psi-norm1.txt','w');       % 创建psi-1范数随n变化输出文件


%E=(4.0/(n-1))*[0:n-1]-2.0;
%E=1.0*ones(1,n);
E=beita+Lamada*cos(2.0*pi*([0:n-1]-29)*0.5*(sqrt(5.0)+1)); %格子点处能量赋值---自己根据情况定

a=0.006*ones(1,n);
b=0.006*ones(1,n);                                        %psi-实部和虚部变量分别声称
                                                         %变量初始化-也为0时刻初始值---自己根据情况定


psi2=0.0*ones(1,n);                                     %psi-2范数变量声称

T=1.0e+10;                                  %时间限制
time=0.0;                                     %开始计时
num=0;                                        %推进步数

while(time
   num=num+1;
   time=time+dt;

   Da=0.0*ones(1,n);
   Db=0.0*ones(1,n);                         %改进Euler法中求f(xn,yn)
   az=0.0*ones(1,n);
   bz=0.0*ones(1,n);                         %改进Euler法中求y(n+1)-中间过渡量

   Daz=0.0*ones(1,n);
   Dbz=0.0*ones(1,n);                      %改进Euler法中求f(x(n+1),y(n+1))



   Da([2:n-1])=Da([2:n-1])+E([2:n-1]).*b([2:n-1])+beita*psi2([2:n-1]).*b([2:n-1])+b([3:n])+b([1:n-2]); %文献1离散方式
   Db([2:n-1])=Db([2:n-1])-E([2:n-1]).*a([2:n-1])-beita*psi2([2:n-1]).*a([2:n-1])-a([3:n])-a([1:n-2]);

% Da([2:n-1])=Da([2:n-1])-E([2:n-1]).*b([2:n-1])-C*(b([3:n])+b([1:n-2]));                         %文献2离散方式
% Db([2:n-1])=Db([2:n-1])+E([2:n-1]).*a([2:n-1])+C*(a([3:n])+a([1:n-2]));

   az([2:n-1])=a([2:n-1])+dt*Da([2:n-1]);
   bz([2:n-1])=b([2:n-1])+dt*Db([2:n-1]);

% az(1)=2.0*A+bz(2);                               %文献1边界条件处理---自己根据情况定
   az(1)=bz(2);                                     %自己给定边界条件
   bz(1)=-az(2);
   az(n)=bz(n-1);
   bz(n)=-az(n-1);


   psi2([2:n-1])=az([2:n-1]).*az([2:n-1])+bz([2:n-1]).*bz([2:n-1]);

   Daz([2:n-1])=Daz([2:n-1])+E([2:n-1]).*bz([2:n-1])+beita*psi2([2:n-1]).*bz([2:n-1])+bz([3:n])+bz([1:n-2]);
   Dbz([2:n-1])=Dbz([2:n-1])-E([2:n-1]).*az([2:n-1])-beita*psi2([2:n-1]).*az([2:n-1])-az([3:n])-az([1:n-2]);

% Daz([2:n-1])=Daz([2:n-1])-E([2:n-1]).*bz([2:n-1])-C*(bz([3:n])+bz([1:n-2]));       %文献2离散方式
% Dbz([2:n-1])=Dbz([2:n-1])+E([2:n-1]).*az([2:n-1])+C*(az([3:n])+az([1:n-2]));

   aaz=0.0*ones(1,n);                                           %物理量更新时中间传值变量声称
   bbz=0.0*ones(1,n);



   aaz([2:n-1])=a([2:n-1])+0.5*dt*(Da([2:n-1])+Daz([2:n-1]));    %改进Euler法一个时间步后物理量更新
   bbz([2:n-1])=b([2:n-1])+0.5*dt*(Db([2:n-1])+Dbz([2:n-1]));

% aaz(1)=2.0*A+bbz(2);                                           %文献1边界条件处理---自己根据情况定
   aaz(1)=bbz(2);                                                 %自己给定边界条件
   bbz(1)=-aaz(2);
   aaz(n)=bbz(n-1);
   bbz(n)=-aaz(n-1);

   error=0.0;                                                    %求下一时刻与上一时刻物理量绝对误差

   error=error+sum(sqrt((aaz(1:n)-a(1:n)).*(aaz(1:n)-a(1:n))+(bbz(1:n)-b(1:n)).*(bbz(1:n)-b(1:n))));



a(1:n)=aaz(1:n);                                           %物理量更新赋值
b(1:n)=bbz(1:n);


psi2(1:n)=a(1:n).*a(1:n)+b(1:n).*b(1:n);                   %物理量更新后Psi-2范数值更新

Da=[]; Db=[];                                              %变量清空
az=[]; bz=[];
Daz=[]; Dbz=[];
aaz=[]; bbz=[];

disp(['推进步数为',num2str(num)]);                         %推进步数输出到屏幕上

if(mod(num,10000)==0)                                     %每隔10000步输出|PSI|n与粒子位置关系
      fprintf(fidout,'zone f=point\n');
      for i=1:n
         fprintf(fidout,'%5d %f\n',i-29,sqrt(psi2(i)));
      end
end

if(error<=MyEps||mod(num,20000000)==0)                      %达到近似稳定状态或运行到2000万步后结果输出、程序运行结束
      fprintf(fidout,'zone f=point\n');
      for i=1:n
         fprintf(fidout,'%5d %f\n',i-29,sqrt(psi2(i)));
      end                                                 %程序运行结束
      fclose(fidout);
      disp(['达到近似稳定后需要推进步数为-------:',num2str(num)]);
      disp(['达到近似稳定后需要时间为***********:',num2str(time)]);
      break; %程序运行结束
end
end

[ Last edited by anyuezhiji on 2010-4-12 at 00:43 ]

回复此楼

暗月下没有留下风的痕迹，但它已经寂然飘逝。。By&amp;lt;暗月之寂&amp;gt;:tiger38: