計算物理課程論文

1、 微分方程的數(shù)值模擬及應用 微分方程的數(shù)值模擬及應用本文介紹了 matlab、Mathematica 等軟件在微分方程數(shù)值模擬上的應用。作為基礎論文首先介紹了用庫塔-龍格方法和有限元差分方法求解一階微分方程組及高階微分方程的方法并給出了實現(xiàn)的matlab 代碼，在了解解微分方程的基本原理之后，本文用Mathematica 軟件研究了一維深勢阱、諧振子的波函數(shù)以及有心力場下的量子力學現(xiàn)像，如原子軌道、分

2、子軌道。接著介紹了一類特殊的微分方程—非線性薛定諤方程 NLSE，這類方程不同與其他微分方程之處在于它存在孤子解，比較復雜。本文介紹了求這類方程數(shù)值解得有限元差分方法及分步傅里葉方法，并給出了一個后者的matlab 實例代碼。最后用 mathematica 對其進行了數(shù)值模擬，研究了其在光波導和光孤子中的應用。1. 求解一階微分方程組及高階微分方程的方法。（1）亞當斯預測 亞當斯預測-校正法求一階常微分方程。 校正法求一階常微分方程。f

3、unction [k,X,Y,wucha,P]=dAdamsyx(funfcn,x0,b,y0,h)x=x0;y=y0;p=128; n=fix((b-x0)/h); if n<5,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:4c1=1/6;c2=2/6;c3=2/6; c4=1

4、/6;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);

5、 y3=y+b41*h*k1+b42*h*k2+b43*h*k3; k4=feval(funfcn,x3,y3);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(k)+h,Y1(k,:)'+k3*h); Y1(k+1,:)=Y1(k,:)

6、+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,

7、 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(:,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;

8、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('是邊值問題的數(shù)值解y(x)的曲線','是初值問題1

9、的數(shù)值解u(x)的曲線', '是初值問題2的數(shù)值解v(x)的曲線')title('用線性打靶法求線性邊值問題的數(shù)值解的圖形')(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=ze

10、ros(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)=-2-(h.^2)*k2(k); A3(k,k+1)= 1-h*k1(k)/2

11、; 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)=(h.^2)*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=A

12、\B1;Y1=Y'; y=[alpha;Y;beta]; 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,1); k=1:n+1; wucha=wucha(1:k,:); plot(X,y(:,1),'mp')xlabel('軸\it x'); ylabel('軸\it y')