基于MATLAB的偏微分方程差分解法

基于MATLAB的偏微分方程差分解法學院：核工程與地球物理學院專業(yè)：勘查技術與工程班級：1120203學號：姓名:2014/6/11在科學技術各領域中，有很多問題都可以歸結為偏微分方程問題。在物理專業(yè)的力學、熱學、電學、光學、近代物理課程中都可遇見偏微分方程。偏微分方程，再加上邊界條件、初始條件構成的數(shù)學模型，只有在很特殊情況下才可求得解析解。隨著計算機技術的發(fā)展，采用數(shù)值計算方法，可以得到其數(shù)值解。近些年來，求解偏微分方程的數(shù)值方法取得進展，特別是有限差分區(qū)域分解算法，此類算法的特點是在內邊界處設計不同于整體的格式，將全局的隱式計算化為局部的分段隱式計算。使人從感覺上認為這樣得到的解會比全局隱式得到的解的精度差，但大量的數(shù)值實驗表明事實正好相反，用區(qū)域分解算法求得的解的精度更好。差分方法又稱為有限差分方法或網(wǎng)格法（網(wǎng)格點（稱為差分格式）（數(shù)值解。因此，用差分方法求偏微分方程定解問題一般需要解決以下問題：選取網(wǎng)格；對微分方程及定解條件選擇差分近似，列出差分格式；求解差分格式；下面對偏微分方程具體例題的差分解法作一簡要的介紹。§1雙曲型方程中波動方程的有限差分解法。1.1雙曲型的差分方程通過建立網(wǎng)格并求解中心差分方程結果為：u (22r2)u r2(u u )u ,i,ji,j i1,j i1,j i,j其中為了保證上式的穩(wěn)定性，必須使rck/h1.2初始值通過聯(lián)立初始值及邊界條件可以得到

i2,3,Ln。u(x,k)fi i

kgi

c2k22h2

(fi1

2fi

fi1

)O(h2)O(k2)O(k3

(1)代入rck/h，可簡化并得到一個改進的對行2的近似值差分方程：ui,2

(1r2)fi

kgi

r2(2

i1

f i1

(2)1.3雙曲型方程中波動方程例題的差分解法結果及程序。題：u(x,t)4u (x,t)，其中0x且0t0.5，邊界條件為：tt xx設h0.2,k0.1,r解：

u(0,t)且u(1,t0,u(x,0)f(x)sin(x),u(x,0)g(x)0,t

0t0.50x10x1第一步：通過聯(lián)立(1)、(2)編寫MATLAB程序如下：%二維雙曲型偏微分方程，使用D'Alembert方法functionU=hyperbolic(a,b,c,n,m)%a為x的取值范圍%b為t的取值范圍%c為系數(shù)%n為x方向上的節(jié)點數(shù)%m為t上的節(jié)點數(shù)h=a/(n-1);%x方向上的步長k=b/(m-1);%t上的步長r=c*k/h;r2=r^2;r22=r^2/2;s1=1-r^2;s2=2-2*r^2;U=zeros(n,m);fori=2:(n-1);U(i,1)=sin(pi*h*(i-1));U(i,2)=s1.*sin(pi*h*(i-1))+k*0...+r22.*(sin(pi*h.*(i-2))+sin(pi*h.*(i))); end%差分方程forj=3:m;fori=2:(n-1);U(i,j)=s2*U(i,j-1)+r2*(U(i+1,j-1)+U(i-1,j-1))-U(i,j-2);%P115(7)endendU=U';figure(1);surf(U);figure(2);contour(U,40);第二步：輸入數(shù)值并計算a=1;b=0.5;c=4n=11m=11執(zhí)行hyperbolic(1,0.5,4,11,11)；第三步：得出結果并畫圖入下等值線結果圖三位結果圖§1.4MATLABpdepe(的調用格式為：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)，具體程序見附錄得出的結果為：等值線結果圖三維結果圖1.5結果對比通過編寫MATLAB的差分方程程序求取結果和MATLAB自帶函數(shù)求取結果進行對比，發(fā)現(xiàn)這兩種方法求得到的結果是非常理想的。§2拋物線方程中熱傳導方程的有限差分解法。拋物線方程的差分方程通過建立網(wǎng)格并求解顯示前向差分方程結果為：u (12r)u r(u u ) (3)i,j1 i,j i-1,j i+1,j2其中為了保證上式前向差分方程穩(wěn)定性，當且僅當r滿足0r1時。這意味著步長k必須2滿足kh2 2c2 。拋物線方程中熱傳導方程例題的差分解法結果及程序。xtt(xt其中0xt初始條件為ux,0f(x)12x1其中t0,0x1,邊界條件為：t1u0,tct1

x0且0t0.1ut

c2

x1且0t0.1解：第一步，分析并帶入（3）并編寫MATLAB求解程序如下：functionU=forwdif(c1,c2,a,b,c,n,m)clch=a/(n-1);k=b/(m-1);r=(c^2*k)/(h^2);s=1-2*r;U=zeros(n,m);U(1,1:m)=c1;U(n,1:m)=c2;fori=2:n-1U(I,1)=1-abs(2*(i-1)*h-1);% U(I,1)=4*(i-1)*h-4*((i-1)*h)^2;endforj=2:mfori=2:n-1U(I,j)=s*U(I,j-1)+r*(U(i-1,j-1)+U(i+1,j-1));endendU=U’;figure(1);surf(U);figure(2);contour(U,30);第二步，代入初始條件以及邊界條件：c1=0；c2=0；a=1；b=0.5；c=1；n=6；m=11；執(zhí)行forwdif(0,0,1,0.1,1,11,11)；第三步：得出結果并畫圖入下1.等值線結果圖2.三位結果圖§2.3MATLABpdepe(的調用格式為：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)，具體程序見附錄得出的結果為：等值線結果圖三維結果圖2.4結果對比通過編寫MATLAB的差分方程程序求取結果和MATLAB自帶函數(shù)求取結果進行對比，發(fā)現(xiàn)這兩種方法求得到的結果是非常相似的，差距不大證明程序編寫是成功的?！?橢圓型方程的有限差分解法。3.1建立線型方程組1u(x,y)u ,1

2jm1

(在左邊)j1,jjux,yu ,

2in

(在底邊)i 1

i,1u x,yn u x,yi

u n,ju i,m

2jm12in1

(在右邊)(在頂邊)2導數(shù)邊界條件Neumann邊界條件確定了u(x,y)邊的法線的方向導數(shù)。這里使用零法線導數(shù)條件：Nu(x,y)0對于熱傳導而言，這表示邊是熱絕緣的而且經過邊的熱通量為零，從而得到：xu(xn,yj)ux(xn,yj)0得出點(x,y的Laplace差分方程為：n j利用差分方程：

u u u u n1,j n1,j n,j1 n,j

(4)un1,j

un1,j2h

u(x,yx n

)0 (5)在（4）上式使用（5）式這個近似值，這結果為：2u u u -4u n1,j n,j+1 n,j-1 n,j

(4)n1,jn,j+1n,j1這個公式講函數(shù)u ，u n1,jn,j+1n,j14種情況如下所示：2ui,2

ui

ui

-4u

i,1

(底部)2ui,m1

ui1,m

ui1,m

-4u

i,m

(頂部)

(5)2u2,j2u

u1,ju

u1,ju

-4u 01,j-4u

(左部)(右部)3.3迭代方法

n1,j

n,j

n,j

n,j根據(jù)（4）式和（5）式以如下形式進行迭代處理：式中：

ui,j

ui,j

ri,

(6)ri,j

ui1,

ui1,

ui,j4

ui,j

i,

(7)這里2xn且2jm1.i,i,i,i,

“減少到零”（|

2in1且2jm1）?？梢岳弥鸫纬沙诜ǎ⊿OR）提高所有余項

ri,j

減少到零的收斂速度。其中逐次超松弛法使用迭代公式：ui,j

ui,ju

r

i1,j

ui1,j

ui,j4

ui,j

i,j

(8)i,j i,j這里參數(shù)位于12范圍內。對所有網(wǎng)格點應用上直到r2 42 4cos cosn12m1

。其中：i,i,j3.4橢圓型方程中波動方程例題的差分解法結果及程序。題：用迭代法求解在區(qū)域rx,y:0x4,0x4Lpalace方程的20近似解，這里邊界值為：ux,020;ux,4解：

0x4第一步：通過聯(lián)立差分方程及迭代方程編寫MATLAB程序如下：functionU=dirich(a,b,h,tol,maxl)%設DX=DY=h,且存在m，是的a=nh和b=mhn=fix(a/h)+1;m=fix(b/h)+1;ave=(a*(20+180)+b*(80+0))/(2*a+2*b);U=ave*ones(n,m);U(1,1:m)=80;U(n,1:m)=0;U(1:n,1)=20;U(1:n,m)=180;U(1,1)=(U(1,2)+U(2,1))/2;U(1,m)=(U(1,m-1)+U(2,m))/2;U(n,1)=(U(n-1,1)+U(n,2))/2;U(n,m)=(U(n-1,m)+U(n,m-1))/2;w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));err=1;cnt=0;while((err>tol)&(cnt<=maxl))err=0;forj=2:m-1fori=2:n-1relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j))/4;U(i,j)=U(i,j)+relx;if(err<=abs(relx))err=abs(relx);end

end

endcnt=cnt+1;relx;endU=flipud(U');figure(1);surf(U);figure(2);contour(U,30);第二步：輸入數(shù)值并計算a=4;b=4;h=0.5;tol=0.1;maxl=20;執(zhí)行dirich(4,4,0.5,0.1,20)；第三步：得出結果并畫圖入下1.等值線結果圖2.三位結果圖§3附錄雙曲型及拋物線方程MATLAB語言提供的pdepe()函數(shù)程序。%寫主函數(shù)clcx=linspace(0,1,50);t=linspace(0,0.1,50);m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);u=sol(:,:,1);figure('numbertitle','off','name','PDEMatlabsky')surf(x,t,u)title('TheSolutionofxlabel('X')ylabel('T')zlabel('U')figure(2);contour(u,40);%目標PDE函數(shù)function[c,f,s]=pdefun(x,t,u,