基于MATLAB的偏微分方程差分解法_第1頁
基于MATLAB的偏微分方程差分解法_第2頁
基于MATLAB的偏微分方程差分解法_第3頁
已閱讀5頁,還剩7頁未讀, 繼續(xù)免費閱讀

下載本文檔

版權說明:本文檔由用戶提供并上傳,收益歸屬內容提供方,若內容存在侵權,請進行舉報或認領

文檔簡介

基于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,

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
  • 4. 未經權益所有人同意不得將文件中的內容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內容本身不做任何修改或編輯,并不能對任何下載內容負責。
  • 6. 下載文件中如有侵權或不適當內容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論