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

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

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

kgi

c2k22h2

(fi1

2fi

fi1

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

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

(1r2)fi

kgi

r2(2

i1

f i1

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

x0且0t0.1ut

c2

x1且0t0.1解：第一步，分析并帶入（3）并編寫(xiě)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)；第三步：得出結(jié)果并畫(huà)圖入下1.等值線結(jié)果圖2.三位結(jié)果圖§2.3MATLABpdepe(的調(diào)用格式為：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)，具體程序見(jiàn)附錄得出的結(jié)果為：等值線結(jié)果圖三維結(jié)果圖2.4結(jié)果對(duì)比通過(guò)編寫(xiě)MATLAB的差分方程程序求取結(jié)果和MATLAB自帶函數(shù)求取結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)這兩種方法求得到的結(jié)果是非常相似的，差距不大證明程序編寫(xiě)是成功的?！?橢圓型方程的有限差分解法。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導(dǎo)數(shù)邊界條件Neumann邊界條件確定了u(x,y)邊的法線的方向?qū)?shù)。這里使用零法線導(dǎo)數(shù)條件：Nu(x,y)0對(duì)于熱傳導(dǎo)而言，這表示邊是熱絕緣的而且經(jīng)過(guò)邊的熱通量為零，從而得到：xu(xn,yj)ux(xn,yj)0得出點(diǎn)(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）式這個(gè)近似值，這結(jié)果為：2u u u -4u n1,j n,j+1 n,j-1 n,j

(4)n1,jn,j+1n,j1這個(gè)公式講函數(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）式以如下形式進(jìn)行迭代處理：式中：

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）提高所有余項(xiàng)

ri,j

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

ui,ju

r

i1,j

ui1,j

ui,j4

ui,j

i,j

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

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

0x4第一步：通過(guò)聯(lián)立差分方程及迭代方程編寫(xiě)MATLAB程序如下：functionU=dirich(a,b,h,tol,maxl)%設(shè)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ù)值并計(jì)算a=4;b=4;h=0.5;tol=0.1;maxl=20;執(zhí)行dirich(4,4,0.5,0.1,20)；第三步：得出結(jié)果并畫(huà)圖入下1.等值線結(jié)果圖2.三位結(jié)果圖§3附錄雙曲型及拋物線方程MATLAB語(yǔ)言提供的pdepe()函數(shù)程序。%寫(xiě)主函數(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);%目標(biāo)PDE函數(shù)function[c,f,s]=pdefun(x,t,u,