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

1、基于MATLAB的偏微分方程差分解法學院：核工程與地球物理學院專業(yè)：勘查技術與工程班級：學號：姓名:2014/6/11在科學技術各領域中，有很多問題都可以歸結(jié)為偏微分方程問題。在物理專業(yè)的力學、熱學、電學、光學、近代物理課程中都可遇見偏微分方程。偏微分方程，再加上邊界條件、初始條件構(gòu)成的數(shù)學模型，只有在很特殊情況下才可求得解析解。隨著計算機技術的發(fā)展，采用數(shù)值計算方法，可以得到其數(shù)值解。近些年來，求解偏微分方程的數(shù)值方法取得進展，特別是有限差分區(qū)域分解算法，此類算法的特點是在內(nèi)邊界處設計不同于整體的格式， 將全局的隱式計算化為局部的分段隱式計算。使人從感覺上認為這樣得到的解會比全局隱式得到的解

2、的精度差，但大量的數(shù)值實驗表明事實正好相反，用區(qū)域分解算法求得的解的精度更好。差分方法又稱為有限差分方法或網(wǎng)格法，是求偏微分方程定解問題的數(shù)值解中應用最廣泛的方法之一。它的基本思想是：先對求解區(qū)域作網(wǎng)格剖分，將自變量的連續(xù)變化區(qū)域用有限離散點（網(wǎng)格點）集代替；將問題中出現(xiàn)的連續(xù)變量的函數(shù)用定義在網(wǎng)格點上離散變量的函數(shù)代替；通過用網(wǎng)格點上函數(shù)的差商代替導數(shù)，將含連續(xù)變量的偏微分方程定解問題化成只含有限個未知數(shù)的代數(shù)方程組（稱為差分格式）。如果差分格式有解，且當網(wǎng)格無限變小時其解收斂于原微分方程定解問題的解，則差分格式的解就作為原問題的近似解（數(shù)值解）。因此，用差分方法求偏微分方程定解問題一般需要

3、解決以下問題：（i）選取網(wǎng)格；（ii）對微分方程及定解條件選擇差分近似，列出差分格式；（iii）求解差分格式；（iv）討論差分格式解對于微分方程解的收斂性及誤差估計。下面對偏微分方程具體例題的差分解法作一簡要的介紹。1 雙曲型方程中波動方程的有限差分解法。1.1 雙曲型的差分方程通過建立網(wǎng)格并求解中心差分方程結(jié)果為：其中為了保證上式的穩(wěn)定性，必須使1. 2 初始值通過聯(lián)立初始值及邊界條件可以得到代入，可簡化并得到一個改進的對行2的近似值差分方程：1.3 雙曲型方程中波動方程例題的差分解法結(jié)果及程序。題：，其中，邊界條件為：設解：第一步：通過聯(lián)立(1)、(2)編寫MATLAB程序如下：%二維雙曲

4、型偏微分方程，使用DAlembert方法function U=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=r2;r22=r2/2;s1=1-r2;s2=2-2*r2;U=zeros(n,m);for i=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);

5、%P116（13） end%差分方程 for j=3:m; for i=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) end end U=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)；第三步：得出結(jié)果并畫圖入下1.等值線結(jié)果圖2.三位結(jié)果圖1.4 通過MATLAB語言提供了pdepe()函數(shù)，可以直接求解一般偏微分方程(組)，它

6、的調(diào)用格式為：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)，具體程序見附錄得出的結(jié)果為：1.等值線結(jié)果圖2.三維結(jié)果圖1.5 結(jié)果對比通過編寫MATLAB的差分方程程序求取結(jié)果和MATLAB自帶函數(shù)求取結(jié)果進行對比，發(fā)現(xiàn)這兩種方法求得到的結(jié)果是非常理想的。2 拋物線方程中熱傳導方程的有限差分解法。2.1 拋物線方程的差分方程通過建立網(wǎng)格并求解顯示前向差分方程結(jié)果為：其中為了保證上式前向差分方程穩(wěn)定性，當且僅當滿足時。這意味著步長必須滿足。2.2 拋物線方程中熱傳導方程例題的差分解法結(jié)果及程序。題：其中初始條件為其中邊界條件為：解：第一步，分析并帶入（3）并編寫MATL

7、AB求解程序如下：function U=forwdif(c1,c2,a,b,c,n,m)clch=a/(n-1);k=b/(m-1);r=(c2*k)/(h2);s=1-2*r;U=zeros(n,m);U(1,1:m)=c1;U(n,1:m)=c2;for i=2:n-1U(I,1)=1-abs(2*(i-1)*h-1);% U(I,1)=4*(i-1)*h-4*(i-1)*h)2;endfor j=2:mfor i=2:n-1 U(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);

8、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)；第三步：得出結(jié)果并畫圖入下1.等值線結(jié)果圖2.三位結(jié)果圖2.3 通過MATLAB語言提供了pdepe()函數(shù)，可以直接求解一般偏微分方程(組)，它的調(diào)用格式為：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)，具體程序見附錄得出的結(jié)果為：1.等值線結(jié)果圖2.三維結(jié)果圖2.4 結(jié)果對比通過編寫MATLAB的差分方程程序求取結(jié)果和MATLAB自帶函數(shù)求取結(jié)果進行對比，發(fā)現(xiàn)這兩種方法求得到的結(jié)果

9、是非常相似的，差距不大證明程序編寫是成功的。3 橢圓型方程的有限差分解法。3.1 建立線型方程組3. 2 導數(shù)邊界條件Neumann邊界條件確定了邊的法線的方向?qū)?shù)。這里使用零法線導數(shù)條件：對于熱傳導而言，這表示邊是熱絕緣的而且經(jīng)過邊的熱通量為零，從而得到：得出點的Laplace差分方程為：利用差分方程：在（4）上式使用（5）式這個近似值，這結(jié)果為：這個公式講函數(shù)聯(lián)系起來。則Neumann計算空格樣板的4種情況如下所示：3. 3 迭代方法根據(jù)（4）式和（5）式以如下形式進行迭代處理：式中：這里通過迭代公式（6）右邊對所有的內(nèi)部結(jié)點連續(xù)執(zhí)行迭代，直到該式右邊余項 “減少到零”（即，）。可以利用逐

10、次超松弛法（SOR）提高所有余項減少到零的收斂速度。其中逐次超松弛法使用迭代公式：這里參數(shù)位于范圍內(nèi)。對所有網(wǎng)格點應用上式（8）直到。其中：3.4 橢圓型方程中波動方程例題的差分解法結(jié)果及程序。題：用迭代法求解在區(qū)域 Lpalace方程的近似解，這里邊界值為：解：第一步：通過聯(lián)立差分方程及迭代方程編寫MATLAB程序如下：function U=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

11、(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(errtol)&(cnt=maxl) err=0; for j=2:m-1 for i=2:n-1 relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i

12、-1,j)-4*U(i,j)/4; U(i,j)=U(i,j)+relx; if(err=abs(relx) err=abs(relx); end end end cnt=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)；第三步：得出結(jié)果并畫圖入下1.等值線結(jié)果圖2.三位結(jié)果圖3 附錄雙曲型及拋物線方程MATLAB語言提供的pdepe()函數(shù)程序。%寫主函數(shù)clcx=

13、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,PDE Demoby Matlabsky)surf(x,t,u)title(The Solution of u_1)xlabel(X)ylabel(T)zlabel(U)figure(2);contour(u,40);% 目標PDE函數(shù)function c,f,s=pdefun (x,t,u,du)%p141-4c=1;f=du;s=0;% %p139-4% c=4;% f=du;% s=0;