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

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

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

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

4、LAB程序如下：%二維雙曲型偏微分方程，使用D'Alembert方法function U=hyperbolic(a,b,c,n,m)% a為x的取值范圍% b為t的取值范圍% c為系數% n為x方向上的節(jié)點數% m為t上的節(jié)點數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-

5、2)+sin(pi*h.*(i); %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);第二步：輸入數值并計算a=1;b=0.5;c=4n=11m=11執(zhí)行hyperbolic(1,0.5,4,11,11)；第三步：得出結果并畫圖入下1.等值線結果圖2.三位結果圖§1.4 通過MATLAB語言提供了

6、pdepe()函數，可以直接求解一般偏微分方程(組)，它的調用格式為：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)，具體程序見附錄得出的結果為：1.等值線結果圖2.三維結果圖1.5 結果對比通過編寫MATLAB的差分方程程序求取結果和MATLAB自帶函數求取結果進行對比，發(fā)現這兩種方法求得到的結果是非常理想的。§2 拋物線方程中熱傳導方程的有限差分解法。2.1 拋物線方程的差分方程通過建立網格并求解顯示前向差分方程結果為：其中為了保證上式前向差分方程穩(wěn)定性，當且僅當滿足時。這意味著步長必須滿足。2.2 拋物線方程中熱傳導方程例題的差分解法結果及程序。題：其中

7、初始條件為其中邊界條件為：解：第一步，分析并帶入（3）并編寫MATLAB求解程序如下：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);enden

8、dU=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.3 通過MATLAB語言提供了pdepe()函數，可以直接求解一般偏微分方程(組)，它的調用格式為：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)，具體程序見附錄得出的結果為：1.等值線結果圖2.三維結果圖2.4 結果對比通過編寫MATLAB的差分方

9、程程序求取結果和MATLAB自帶函數求取結果進行對比，發(fā)現這兩種方法求得到的結果是非常相似的，差距不大證明程序編寫是成功的。§3 橢圓型方程的有限差分解法。3.1 建立線型方程組3. 2 導數邊界條件Neumann邊界條件確定了邊的法線的方向導數。這里使用零法線導數條件：對于熱傳導而言，這表示邊是熱絕緣的而且經過邊的熱通量為零，從而得到：得出點的Laplace差分方程為：利用差分方程：在（4）上式使用（5）式這個近似值，這結果為：這個公式講函數聯系起來。則Neumann計算空格樣板的4種情況如下所示：3. 3 迭代方法根據（4）式和（5）式以如下形式進行迭代處理：式中：這里通過迭代公

10、式（6）右邊對所有的內部結點連續(xù)執(zhí)行迭代，直到該式右邊余項 “減少到零”（即，）?？梢岳弥鸫纬沙诜ǎ⊿OR）提高所有余項減少到零的收斂速度。其中逐次超松弛法使用迭代公式：這里參數位于范圍內。對所有網格點應用上式（8）直到。其中：3.4 橢圓型方程中波動方程例題的差分解法結果及程序。題：用迭代法求解在區(qū)域 Lpalace方程的近似解，這里邊界值為：解：第一步：通過聯立差分方程及迭代方程編寫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*

11、(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; for

12、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-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);第二步：輸入數值并計算a=4; b=4; h=0.5; tol=0.1; maxl=20;執(zhí)行dirich(4,4,0.5,0.1,20)；第三步：得出結

13、果并畫圖入下1.等值線結果圖2.三位結果圖§3 附錄雙曲型及拋物線方程MATLAB語言提供的pdepe()函數程序。%寫主函數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','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函數function c,f,s=pdefun (x,t,u,du)%p141-4c=1;f=du;s=0;% %p139-4% c=4;% f=du