计算物理课程论文

1、微分方程的数值模拟及应用本文介绍了 matlab、Mathematica等软件在微分方程数值模拟 上的应用。作为基础论文首先介绍了用库塔-龙格方法和有限元差分 方法求解一阶微分方程组及高阶微分方程的方法并给出了实现的 matlab代码，在了解解微分方程的基本原理之后，本文用Mathematica 软件研究了一维深势阱、谐振子的波函数以及有心力场下的量子力学 现像，如原子轨道、分子轨道。接着介绍了一类特殊的微分方程一非 线性薛定谔方程NLSE，这类方程不同与其他微分方程之处在于它存 在孤子解，比较复杂。本文介绍了求这类方程数值解得有限元差分方 法及分步傅里叶方法，并给出了一个后者的matlab实

2、例代码。最后 用mathematica对其进行了数值模拟，研究了其在光波导和光孤子中 的应用。求解一阶微分方程组及高阶微分方程的方法。(1)亚当斯预测核正法求一阶常微分方程。function k,X,Y,wucha,P=dAdamsyx(funfcn,x0,b,y0,h) x=x0;y=y0;p=128; n=fix(b-x0)/h); if n5, return, end; X=zeros(p,1); Y=zeros(p,length(y); f=zeros(p,1);k=1; X(k)=x; Y(k,:)=y; for k=2:4 c1=1/6;c2=2/6;c3=2/6; c4=1/6;

3、a2=1/2; a3=1/2; a4=1;b21=1/2;b31=0;b32=1/2; b41=0;b42=0;b43=1; x1=x+a2*h; x2=x+a3*h; x3=x+a4*h; k1=feval(funfcn,x,y); y1=y+b21*h*k1; x=x+h; k2=feval(funfcn,x1,y1); y2=y+b31*h*k1+b32*h*k2; k3=feval(funfcn,x2,y2); y3=y+b41*h*k1+b42*h*k2+b43*h*k3; k4=feval(funfcn,x3,y3);y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4);

4、X(k)=x; Y(k,:)=y; endX;Y;f=feval(funfcn,X(1:4),Y(1:4);f=f, for k=4:nX(k+1)=X(1)+h*k;f(k)=feval(funfcn,X(k),Y(k);P=Y(k)+(h/24)*(f(k-3:k)*-9 37 -59 55);f=f(2) f(3) f(4) feval(funfcn,X(k+1),P), Y(k+1)=Y(k)+(h/24)*(f*1 -5 19 9);f(4)= feval(funfcn,X(k+1),Y(k+1);k=k+1; end for k=1:nwucha(k+1)=norm(Y(k+1)-

5、Y(k);endX=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1,wucha=wucha(1:n,:);P=n,X,Y,wucha;四阶库塔龙格方法解一阶微分方程组function k,X,Y,wucha,P=RK4z(dydx,a,b,CT,h) n=fix(b-a)/h);X=zeros(n+1,1);Y=zeros(n+1,length(CT);X=a:h:b;Y(1,:)= CT;for k=1:nk1=feval(dydx,X(k),Y(k,:) x2=X(k)+h/2;y2=Y(k,:)+k1*h/2; k2=feval(dydx,x2,y2);k3=feval(dy

6、dx,x2,Y(k,:)+k2*h/2);k4=feval(dydx, X(k)+h,Y(k,:)+k3*h);Y(k+1,:)=Y(k,:)+h*(k1+2*k2+2*k3+k4)/6;k=k+1; end for k=2:n+1 wucha(k)=norm(Y(k)-Y(k-1); k=k+1;endX=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:); P=k,X,Y,wucha;求解高阶微分方程线性边值问题的线性打靶法functionk,X,Y,wucha,P=xxdb(dydx1,dydx2,a,b,alpha,beta,h) n=fi

7、x(b-a)/h); X=zeros(n+1,1); CT1=alpha,0; Y=zeros(n+1,length(CT1); Y1=zeros(n+1,length(CT1); Y2=zeros(n+1,length(CT1);X=a:h:b;Y1(1,:)= CT1;CT2=0,1;Y2(1,:)= CT2; for k=1:nk1=feval(dydx1,X(k),Y1(k,:) x2=X(k)+h/2;y2=Y1(k,:)+k1*h/2; k2=feval(dydx1,x2,y2);k3=feval(dydx1,x2,Y1(k,:)+k2*h/2);k4=feval(dydx1, X

8、(k)+h,Y1(k,:)+k3*h);Y1(k+1,:)=Y1(k,:)+h*(k1+2*k2+2*k3+k4)/6,k=k+1; end u=Y1(:,1) for k=1:n k1=feval(dydx2,X(k),Y2(k,:) x2=X(k)+h/2;y2=Y2(k,:)+k1*h/2; k2=feval(dydx2,x2,y2); k3=feval(dydx2,x2,Y2(k,:)+k2*h/2); k4=feval(dydx2, X(k)+h,Y2(k,:)+k3*h);Y2(k+1,:)=Y2(k,:)+h*(k1+2*k2+2*k3+k4)/6,k=k+1; end v=Y2

9、(:,1) Y=u+(beta-u(n+1)*v/v(n+1) for k=2:n+1 wucha(k)=norm(Y(k)-Y(k-1); k=k+1; end X=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:); P=k,X,Y,wucha; plot(X,Y(:,1),ro,X,Y1(:,1),g*,X,Y2(:,1),mp) xlabel(轴it x); ylabel(轴it y)legend(是边值问题的数值解y(x)的曲线，是初值问题1的数值解u(x)的曲 线，是初值问题2的数值解v(x)的曲线)title(用线性打靶法求线性边值问

10、题的数值解的图形)(4)求解高阶微分方程的有限差分方法。function k,A,B1,X,Y,y,wucha,p=yxcf(q1,q2,q3,a,b,alpha,beta,h) n=fix(b-a)/h); X=zeros(n+1,1); Y=zeros(n+1,1); A1=zeros(n,n); A2=zeros(n,n); A3=zeros(n,n); A=zeros(n,n);B= zeros(n,1); for k=1:n X=a:h:b; k1(k)=feval(q1,X(k); A1(k+1,k)=1+h*k1(k)/2; k2(k)=feval(q2,X(k);A2(k,k)

11、=-2-(h2)*k2(k);A3(k,k+1)= 1-h*k1(k)/2; k3(k)=feval(q3,X(k); end for k=2:n B(k,1)=(h.八2)*k3(k); end B(1,1)=(h.八2)*k3(1)-(1+h*k1(1)/2)*alpha; B(n-1,1) = (h2)*k3(n-1)-(1+h*k1(n-1)/2)*beta; A=A1(1:n-1,1:n-1)+A2(1:n-1,1:n-1)+A3(1:n-1,1:n-1); B1=B(1:n-1,1); Y=AB1;Y1=Y; y=alpha;Y;beta; for k=2:n+1 wucha(k)

12、=norm(y(k)-y(k-1); k=k+1; end X=X(1:n+1); y=y(1:n+1,1); k=1:n+1; wucha=wucha(1:k,:); plot(X,y(:,1),mp) xlabel(轴 it x); ylabel(轴 it y),legend(是边值问题的数值 解y(x)的曲线)title （-用有限差分法求线性边值问题的数值解的图形）， p=k,X,y,wucha；(5)椭圆型偏微分方程有限差分法function FD_PDE(fun,gun,a,b,c,d)%用有限差分法求解矩形域上的Poisson方程tol=10八(-6);%误差界N=1000;%最

13、大迭代次数n=20;%轴方向的网格数m=20;%轴方向的网格数h=(b-a)/n; %聋由方向的步长l=(d-c)/m; % y轴方向的步长 for i=1:n-1 x(i)=a+i*h;end %定义网格点坐标for j=1:m-1y(j)=c+j*l;end %定义网格点坐标u=zeros(n-1,m-1); 对 u 赋初值%下面定义几个参数r=h八2/l八2;s=2*(1+r);k=1;%应用Gauss-Seidel法求解差分方程while knorm;norm=abs(u(i,m-1)-z); end u(i,m-1)=z;end对右上角的网格点进行处理z=(-h八2*fun(x(n-

14、1),y(m-1)+gun(b,y(mT)+r*gun(x(nT),d)+r*u(nT, m-2)+u(n-m-1)/s;if abs(u(n-1,m-1)-z)normnorm=abs(u(n-1,m-1)-z); end u(n-1,m-1)=z;%对不靠近上下边界的网格点进行处理 for j=m-2:-1:2对靠近左边界的网格点进行处理z=(-h八2*fun(x(1),y(j)+gun(a,y(j)+r*u(1,j+1)+r*u(1,j-1)+u(2,j)/ s;if abs(u(1,j)-z)normnorm=abs(u(1,j)-z);endu(1,j)=z;对不靠近左右边界的网格点

15、进行处理 for i=2:n-2z=(-h八2*fun(x(i),y(j)+u(i-1,j)+r*u(i,j+1)+r*u(i,j-1)+u(i+1,j)/s ;if abs(u(i,j)-z)normnorm=abs(u(i,j)-z);endu(i,j)=z;end对靠近右边界的网格点进行处理z=(-h八2*fun(x(n-1),y(j)+gun(b,y(j)+r*u(n-1,j+1)+r*u(nT,jT)+u( n-2,j)/s;if abs(u(n-1,j)-z)normnorm=abs(u(n-1,j)-z);endu(n-1,j)=z;end汛寸靠近下边界的网格点进行处理对左下角的

16、网格点进行处理z=(-h八2*fun(x(1),y(1)+gun(a,y(1)+r*gun(x(1),c)+r*u(1,2)+u(2,1) /s;if abs(u(1,1)-z)normnorm=abs(u(1,1)-z);endu(1,1)=z;对靠近下边界的除第一点和最后点外网格点进行处理 for i=2:n-2z=(-h八2*fun(x(i),y(1)+r*gun(x(i),c)+r*u(i,2)+u(i+1,1)+u(i-1,1)/ s;if abs(u(i,1)-z)normnorm=abs(u(i,1)-z);endu(i,1)=z;end对右下角的网格点进行处理z=(-h八2*f

17、un(x(n-1),y(1)+gun(b,y(1)+r*gun(x(n-1),c)+r*u(n-1,2)+u (n-2,1)/s;if abs(u(n-1,1)-z)normnorm=abs(u(n-1,1)-z);endu(n-1,1)=z;结果输出if norm CT=1;1;1;h=0.25;k,X,Y,wucha,P=RK4z(dydx,0.1,60,CT,h),plot(X,Y(:,1),g-,X,Y(:,2),b*,X,Y(:,3),mp)xlabel(轴it x); ylabel(轴it y)legend(是方程解z1的曲线,是方程解z2的曲线，是方程解z3的曲线)是方程解rl是

18、方程解rl的曲线 未是方程解口的由线 女 是方程解W的拊姓由于高阶微分方程可以写成一阶微分方程组的形式，故解一阶微分方 程组的库塔-龙格函数可用于解高阶微分方程。例如上例可写成z =x -1 z 一 3 x -2 z + 2 x( -3) z + 9 x3 sin( x) iiii2.在了解库塔龙格方法和有限差分方法的基础之上，2.在了解库塔龙格方法和有限差分方法的基础之上，用matlab、Mathematica求解薛定谔方程1. 一维深势阱。y (1. 一维深势阱。y (x)J力2 d2、2 日 dx2j0, 0 x aV (r )= R (r)Y (9 ,。), R (r )= lmr3,

19、 otherw iseMathematica 代码:x -s=DSolve x,x-k八2*dx口0,d,x / Flatten结果：dTFunctionx, k x C1+-k x C2一维谐振子。d叩 . 八-+ (入一 E2 )v = 0d匕2Mathematics 代码:s=DSolve x,x X +(X-x2)*Dx0,d,x / Flatten结果： dTFunctionx,C2 ParabolicCylinderD1/2(-1-X )严 x+C1 ParabolicCylinderD1/2 (-1$), x3.氢原子轨道计算及模拟：/ l (l + 1)力 2 ) + V (r

20、)+2 日 r2 Ju (r )= Eu (r ),0 r Complexa,-b;AngularDistributionl ,m ,t :=yl,m ccyl,mMathematica 运行：AngularDistributionl_,m_,t_:=yl,m ccyl,mPolarPlot%,t,0,2*Pi2.混态波函数w n l(r)= mCmWn l m(r)2.混态波函数w n l(r)= mCmWn l m(r)Cmyim|2a e ibm Y_imI2Mathematics 代码:(Apply#1 ExpI,l ,m ,t ,p :=Fa定 义 ： (Apply#1 ExpI,l

21、 ,m ,t ,p :=Fa#2&,c,l(yl,#&/m);AngularDistributionMixc ctorComplexExpandExpandmixedc,l,m ccmixedc,l,m ;运行：Pa_,b_=AngularDistributionMix1,0,a,b,1,1,-1,t,p结果：(3 (1+a2-2 a Cosb-2 p) Sint2)/(8 兀)运行：SphericalPlot3DP0.3,0,t,0,Pi,p,0,2Pi结果：氢原子波函数：定义：AnglePlotl_,m_:=Blocktheta,phi, pl=AbsSphericalHarmonicYl

22、,m,theta,phi人2; ParametricPlot3D-pl Sintheta Cosphi,-pl SinthetaSinphi,pl Costheta,phi,0,2Pi,theta,0,Pi,PlotRange-All,PlotPoints-40,40 运行：AnglePlot1,1,AnglePlot1,0,AnglePlot1,-11 File Edi t Ins er t Format Cell Graphics Evaluat i oii Palettes Window Helpn Quantu issues. nb hydrogen atomic p - orbita

23、l :(AnglePlotl, 1, AnglePlotl, 0, AnglePlotl, -1(AnglePlot2/ 2, AnglePlot2, 1, AnglePlot2, -1 , AnglePlot2z 0, AnglePlot2, -2hydrogen, atomic f - orbital :(AnglePlottS, 2, AnglePlot3, 1, AnglePlot3, -1, AnglePlot3, 0, AnglePlot3, -2, AnglePlot3, 3, AnglePlot3z -3定义：Fal_,alm_,a3_,a3m_,blm_,b3_,b3m_:=

24、AngularDistributionMixal,0 ,alm,blm,a3,b3,a3m,b3m,3,1,-1,3,-3,t,p;运行SphericalPlot3DEvaluateF0.3,1.0,0.8,0.2,Pi,Pi/5,Pi/4,t,0,Pi ,p,0,2Pi,PlotRangeAll,PlotPoints40,DisplayFunction-Identit y结果：0.44.分子轨道(成键和反键)定义：GL_,K_,S_,tmod_:=GL,K=1/(2S+1) Sum(2J1+1) (2J+1)*SixJSymbolL,J1,S,J,L,K人2 *Cos(J1-J) *tmod

25、,J1,-AbsL-S,L+S,J,-AbsL-S,L+S;StL_,K_,q_:=StL,K,q=Sum(-1)八(L-(m+q)Sqrt2K+1*ThreeJSymbol L,m+q,L,-m,K,qcm+qcccm,m,-L,L TotalDistributionn_,l_,m_,c_,r_,t_,p_:=RadialDistributionn,l,r AngularDistributionMixc,l,m,t,pAtomwFn_,l_，m_,x_,y_,z_：=Sqrtx人2+y人2+z人2* Radialn,l,Sqrtx人2+y人2+z人2 *yl,m/PolarToCartesi

26、an;BondingDensitynl_,l1_,ml_,n2_ ,l2_,m2_: = (AtomwFnl,l1,ml,x-2,y,z+AtomwFn2,l2,m2,x+2, y,z)人2；AntiBondingDensitynl_,l1_,ml_,n2_,l2_,m2_: = (AtomwFnl,l1,ml,x-2 ,y,z-AtomwFn2,l2,m2,x+2,y,z)八2;运行：ContourPlot3DEvaluateAbsBondingDensity3,0,0,3,0,0,x,-5,5 ,y,-5,5,z,-5,5,Contours 0.005运行：ContourPlot3DEva

27、luateAbsAntiBondingDensity3,0,0,3,0,0,x,- 5,5,y,-5,5,z,-5,5,Contours 0.005 ContourPlot3DEvaluateAbsAntiBondingDensity3,1,1,3,2,2,x,-5,5,y,-5,5,z,-10,10,Contours 0.005 3.非线性薛定谔方程NLSE这类方程在物理研究领域有广泛的应用。非线性光学：光脉冲在色散与非线性介质中的传输非线性光学的自陷现象.光纤孤子与光纤孤子通讯凝聚态物理：热脉冲的传播激光束中原子的Bose-Einstein凝聚效应电磁学.超导电子在电磁场中的运动方程基本形

28、式如：i u = u + 2 lu ，u + + I q I2 q = 0 X、6Z 2 a T 2解非线性薛定谔方程有有限差分法和分步傅里叶算法。有限差分法基本思想与步骤：1。采用一定网格划分方式离散化场域2。基于差分原理，对场域内的偏微分方程以及定解条件进行差分离散化，对每一个离散的节点列出差分方程。3。利用迭代法求解差分方程组4。利用插值法，从离散解得到定解问题在场域内的近似解。q +1 一 q I 1 q +1 一 2 q，+1 + q i q ，一 2 q，+ q ，1i j + （ 7二+ 7 ） + （I q i I2 q i + I q +1 I2 q +1） = 02 A z

29、22 A122 A122 ., ., .,.,分步傅里叶算法：(D -f N)U - = NU（二 t） = E （二=二丁） exp（沁 ry/T!I f-rOC （二仞）=u（二 T） = I U （z. co） exp （-?7）下面是一个该算法的实例：描述一束光在纤维中运动由于色散作用和 非线性效应的相互作用形成孤子的模型。d Aa b i d2 A . 2d z24 d T 2 京a为光纤损耗系数，b是色散系数，g是非线性系数;代码：程序运行结果：代码：程序运行结果：a=0.0clc; clear all; close all; clf;cputime=0;tic;ln=1;i=sq

30、rt(-1);Po=.00064;alpha=0;alph=alpha/(4.343);gamma=0.003;to=125e-12;C=-2;b2=-20e-27;Ld=(to八2)/(abs(b2);pi=3.1415926535;Ao=sqrt(Po);%tau =- 4096e-12:1e-12: 4095e-12; % dt=t/todt=1e-12;rel_error=1e-5;h=1000;for ii=0.1:0.1:1.5z=ii*Ld;u=Ao*exp(-(1+i*(-C)/2)*(tau/to).八2);figure(1)plot(abs(u),r);title( Inp

31、ut Pulse ) ; xlabel( Time ) ; ylabel( Amplitude); grid on; hold on;l=max(size(u);%fwhm1=find(abs(u)abs(max(u)/2);fwhm1=length(fwhm1);dw=1/l/dt*2火pi;w=(-1*l/2:1:l/2-1)*dw;u=fftshift(u);w=fftshift(w);spectrum=fft(fftshift(u);for jj=h:h:zspectrum=spectrum.*exp(-alph*(h/2)+i*b2/2火w.八2*(h/2); f=ifft(spec

32、trum);f=f.*exp(i*gamma*(abs(f).八2)*(h);spectrum=fft(f);spectrum=spectrum.*exp(-alph*(h/2)+i*b2/2火w.八2*(h/2); endf=ifft(spectrum);op_pulse(ln,:)=abs(f);fwhm=find(abs(f)abs(max(f)/2);fwhm=length(fwhm);ratio=fwhm/fwhm1;pbratio(ln)=ratio;dd=atand(abs(imag(f)/(abs(real(f);phadisp(ln)=dd;ln=ln+1;endtoc;cp

33、utime=toc;figure(2);mesh(op_pulse(1:1:ln-1,:);title(Pulse Evolution);xlabel(Time); ylabel(distance); zlabel(amplitude);grid on;hold on;disp(CPU time:), disp(cputime);运行结果:Input PulseTime分步傅里叶算法对于解这类非线性微分方程是一种通用的方法，但是 对于一些方程可直接用mathematica的命令：1.五u/8t-0.5p 82u/8x2+2|u|2 u=0.1(1-cos(兀 x/L)运行：L=50;Table

34、NDSolveIDut,x,t-0.5栉*Dut,x,x,x+2Absut,x八2 ut,x口0.1 (1-CosPi x/L),u0,xSechx ExpI x,ut,-L口ut,L,u,t,0,10,x,-L,L,p,1,5,12.含耗散项的非线性薛定谔方程：L=50;T=0.05;TableNDSolveI Dut,x,t+0.5栉*Dut,x,x,x+2Absut,x八2ut,x+I*T*ut,x0,u0,xSechxExpIx,ut,-L口ut,L,u,t,0,10,x,-L,L,p,1,5,1TablePlot3DEvaluateAbsut,x/.%,t,0,10,x,-L,L,P

35、lotPoint s60,MaxRecursionT 3,MeshNone,p,1,5,13 .耦合的非线性薛定谬方程：代码：8=0.5 ,K12=0.5, coll=l, a 1=1, col2=l, a2=2,K22=0.5, cd22=1 , co21=1) L=15;s=NDSolve 於DQ1 &,!； ,g +於g,T ,c -K12*DQ1,t,t - (coll*Exp -al*g * (Abs QI)八2+col2*Expa2*g * (Abs Q2 g,c)八2)g,tD0, I*DQ2 g,c ,g -於b*DQ2 g,c ,c -K22*DQ2,t,t - (co22*Exp-a2柘* (Abs Q2 g,们)八2+co21*Exp-al*g * (Abs QI g,c ) A2)