二維波動(dòng)方程的有限差分法

1、學(xué)生實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)課程名稱偏微分方程數(shù)值解開(kāi)課實(shí)驗(yàn)室 數(shù)統(tǒng)學(xué)院學(xué)院 數(shù)統(tǒng)年級(jí) 20院專業(yè)班 信計(jì)02班學(xué)生姓名 學(xué) 號(hào)開(kāi)課時(shí) 間 2015至2016學(xué)年第 2 學(xué)期總成績(jī)教師簽名數(shù)學(xué)與統(tǒng)計(jì)學(xué)院制x,y,t=meshgrid(0:h:1,0:h:1,0:tau:1.4);%空間網(wǎng)格剖分開(kāi)課學(xué)院、實(shí)驗(yàn)室： 數(shù)統(tǒng)學(xué)院實(shí)驗(yàn)時(shí)間：2016年6月20日實(shí)驗(yàn)項(xiàng)目名 稱二維波動(dòng)方程的有限差分法實(shí)驗(yàn)項(xiàng)目類型驗(yàn)證演示綜合設(shè)計(jì)其他指導(dǎo)教師目力成 績(jī)是一.實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)該實(shí)驗(yàn)，要求學(xué)生掌握求解二維波動(dòng)方程的有限差分法， 并能通過(guò)計(jì)算機(jī)語(yǔ)言編程 實(shí)現(xiàn)。.實(shí)驗(yàn)內(nèi)容考慮如下的初值問(wèn)題:-2_2_2鬟=劈十粵(x,y)wc =(

2、0,1)2,t10,1.4):t二 x: y22«u (x, y,0 產(chǎn)sinxsin y y,一u (x, y,0 )= 0,(x, y 產(chǎn)(0,1)x,y,t =0, x,y r,t10,1.41在第三部分寫(xiě)出問(wèn)題（1）三層顯格式。2.3.4.根據(jù)你寫(xiě)出的差分格式，編寫(xiě)有限差分法程序。將所寫(xiě)程序放到第四部分。取h =0.1,7 =0.1h ,分別將t =0.5,1.0,1.4時(shí)刻的數(shù)值解畫(huà)圖顯示。該問(wèn)題的解析解為u（x,y,t ） = cosJ2ntsinnxsinny ,將四個(gè)時(shí)刻的數(shù)值解的誤差畫(huà)圖顯示，對(duì)數(shù)值結(jié)果進(jìn)行簡(jiǎn)單的討論。.實(shí)驗(yàn)原理、方法（算法）、步驟網(wǎng)格劃分h =0.

3、1,7 =0.1h，故N =1.4htk =kT , k =0,1,111,140。在內(nèi)網(wǎng)點(diǎn)(K, yj,tk )= 10, M = 140, Xi =ih,yj = jh, i, j =0,1川，10 ,利用二階中心差商，對(duì)（1）建立差分格式:k 1 八 kUi,j-2ui,j2k -1kk kkk k' ui,j =ui1,j-2ui,j ' ui 4,j . ui,j 1 - 2ui,j ' ui,j 4h2h2(2)整理得到:k 12 kkui,j =r ui u ' IUM Uikj42-4r2 u：j y7(3)其中，i, j =12|,9,k=12

4、|,139 ,網(wǎng)比 r=?=0.1,局部截?cái)嗾`差為 o（d + h2）。 h考慮邊界條件u x,y,t =0, x,y）三1"t . 0,1.41,差分格式為：_ _k _ _ k_ _ k_ _ ku0,0 = u0,N = uN,0 = uN ,N = 0,k =0,1,H|,140(4)考慮初始條件u（x,y,0）= sinnxsinny ,差分格式為:0ui,j =sin 二 xi sin r： yj =sin 二 ih sin 二 jh ,i, j = 0,1,l|l ,10(5)2考慮初始條件ut（x,y,0 ）=0,（x,y產(chǎn)（0,1）,利用二階差冏近似:1 jUi j

5、 - Uj j,j,j =0,i, j =0,1,111,102.設(shè)k=0時(shí)刻的點(diǎn)為內(nèi)點(diǎn)，則滿足差分格式(2),代入上式得到:u1,j =r2(u*j +u0,j +Wj+u0j)+(24r2 u0j uj(6)將（6）得到的結(jié)果u：j =5：代入（7）中，整理得到:1ui,j綜上（2）、(4)、k 1uui,j12000020=2rUi 1,jujui,j 1 Ui,j1-2rui,j(5)、(8)得到三層顯格式的差分格式為：2kkkk2k krui 1,j uj "j 1 U,j2-4r uJ-u5(8)0 u ui,ju1,ji,j =1,2,|l|,9,k =1,2,小，13

6、9k k k ku0,0 = u0,N =uN,0 =uN,N - 0, k - 0,1,| I I ,140二sin - xi sin 二 yj l =sin 二ih sin 二 jh ,i, j =0,1,|l|,10= 2r2(u0+j +u0,j +u0j + +ui0j)+(12r2 )u0j,i, j =0,1川,10(9)其中=三=0.1,局部截?cái)嗾`差為o（T2+h2）。 h四.實(shí)驗(yàn)環(huán)境（所用軟件、硬件等）及實(shí)驗(yàn)數(shù)據(jù)文件Matlab%二維波動(dòng)方程數(shù)值計(jì)算（關(guān)鍵：怎么運(yùn)用i,j,k 三個(gè)指標(biāo)建立循環(huán)）clc;%可以將代碼換成函數(shù)m文件h=0.1;tau=0.1*h;% 定義步長(zhǎng) r

7、=tau/h;% 網(wǎng)比uu=cos(sqrt(2)*pi*t).*sin(pi*x).*sin(pi*y);%精確解計(jì)算%第一層網(wǎng)點(diǎn)計(jì)算u=sin(pi*x).*sin(pi*y);% 初始條件u1=u(:,:,1);%因?yàn)榇藭r(shí)得到的u為11x11x141,故只取第一層%第二層網(wǎng)點(diǎn)計(jì)算fo門(mén)=2:10for j=2:10u(i,j,2)=0.5*rA2*(u(i+1,j,1)+u(i-1,j,1)+u(i,j+1,1)+u(i,j-1,1)+(1-2*rA2)*u(i,j, 1)；u(11,:,2)=0;u(:,11,2)=0;endendu2=u(:,:,2);%S 3-141層網(wǎng)點(diǎn)計(jì)算fo

8、r k=2:140for i=2:10for j=2:10u(i,j,k+1)=rA2*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k)+(2-4*rA2)*u(i,j,k)-u(i,j,k-1)；u(11,:,k+1)=0;u(:,11,k+1)=0;endendend%wucha=abs(u-uu);% 求絕對(duì)誤差矩陣 11x11x141wucha1=wucha(:,:,11);% 計(jì)算t=0.1時(shí)刻的絕對(duì)誤差矩陣11x11wucha2=wucha(:,:,51);% 計(jì)算t=0.5時(shí)刻的絕對(duì)誤差矩陣11x11wucha3=wucha(:,:,101

9、);% 計(jì)算t=1.0時(shí)刻的絕對(duì)誤差矩陣11x11wucha4=wucha(:,:,141);% 計(jì)算 t=1.4 時(shí)亥1I的絕對(duì)誤差矢口陣 11x11x0=0:h:1;y0=0:h:1;% %作t=0.1時(shí)刻的絕對(duì)誤差圖 subplot（2,2,1）;mesh（x0,y0,wucha1）;title（'t=0.1時(shí)刻的絕對(duì)誤差'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（' 絕對(duì)誤差值'）； %作t=0.5時(shí)刻的絕對(duì)誤差圖 subplot（2,2,2）;mesh（x0,y0,wucha2）;ti

10、tle（'t=0.5時(shí)刻的絕對(duì)誤差'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（' 絕對(duì)誤差值'）； %乍t=1.0時(shí)刻的絕對(duì)誤差圖subplot（2,2,3）;mesh（x0,y0,wucha3）;title（'t=1.0時(shí)刻的絕對(duì)誤差'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（' 絕對(duì)誤差值'）； %作t=i.4時(shí)刻的絕對(duì)誤差圖 subplot（2,2,4）;mesh（x0,y0,wucha4）

11、;title（'t=1.4時(shí)亥1J的絕對(duì)誤差'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（' 絕對(duì)誤差值'）；%|解 % %t=0.1、0.5時(shí)刻的數(shù)值解與精確解 subplot（2,2,1）;mesh（x0,y0,u（:,:,11）;% 作 t=0.1 時(shí)刻的數(shù)值解title（'t=0.1時(shí)刻的數(shù)值解'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（'u 值'）；subplot（2,2,2）;mesh（

12、x0,y0,uu（:,:,11）;% 作 t=0.1 時(shí)刻的精確解title（'t=0.1 時(shí)刻的精確解）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（'u 值'）；%t=0.5時(shí)刻的數(shù)值解與精確解subplot（2,2,3）;mesh（x0,y0,u（:,:,51）;% 作 t=0.5 時(shí)刻的數(shù)值解title（'t=0.5時(shí)刻的數(shù)值解'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（'u 值'）；subplot（2

13、,2,4）;mesh（x0,y0,uu（:,:,51）;% 作 t=0.5 時(shí)刻的精確解title（'t=0.5 時(shí)刻的精確解'）;%xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（'u 值'）；%?%t=1.0、1.4時(shí)刻的數(shù)值解與精確解subplot（2,2,1）;mesh（x0,y0,u（:,:,101）;% 作 t=1.0 時(shí)刻的數(shù)值解title（'t=1.0 時(shí)刻的數(shù)值解'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel

14、（'u 值'）；subplot（2,2,2）;mesh（x0,y0,uu（:,:,101）;% 作 t=1.0 時(shí)刻的精確解title（'t=1.0時(shí)刻的精確解'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（'u 值'）；%t=1.4時(shí)刻的數(shù)值解與精確解subplot（2,2,3）;mesh（x0,y0,u（:,:,141）;% 作 t=1.4 時(shí)刻的數(shù)值解title（'t=1.4 時(shí)刻的數(shù)值解'）;xlabel（'x 變量'）;ylabel（'

15、;y 變量'）;zlabel（'u 值'）；subplot（2,2,4）;mesh（x0,y0,uu（:,:,141）;% 作 t=1.4 時(shí)刻的精確解title（'t=1.4時(shí)刻的精確解'）;xlabel（'x 變量'）;ylabel（'y 變量'）;zlabel（'u 值'）；五.實(shí)驗(yàn)結(jié)果及實(shí)例分析1、t =0.1、051.01.4時(shí)刻的數(shù)值解與精確解圖圖1 t=0.1、0.5時(shí)刻的數(shù)值解、精確解圖2 t=1.0、1.4時(shí)刻的數(shù)值解、精確解注：上兩圖為四個(gè)時(shí)刻的數(shù)值解與精確解,1r=0.1<j=（

16、p代表維數(shù)）,本文p = 2,三層顯格式達(dá)二階收斂，不難看出，收斂效果很好，符合理論。下圖是四個(gè)時(shí)刻的絕對(duì)誤差圖像，從圖中看出，絕對(duì)誤差較小，且經(jīng)過(guò)計(jì)算得到，收斂階近似于2,正好符合理論值2、t=0.1、051.R 1.4時(shí)刻的絕對(duì)誤差圖t=1。劃的拾"一蝴弱的施好總NQ1時(shí)用的也甘心及M 5m的箔燈我把n W0.02作小序.0.1?10006 .0 .-3,1口圖3四個(gè)時(shí)刻的絕對(duì)誤差3、四個(gè)時(shí)刻（t=0.1、0.5、1.0、1.4）的絕對(duì)誤差表t=0.1時(shí)刻的絕對(duì)誤差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00

17、000.00000.00000.00010.00010.00020.00020.00020.00020.00020.00010.00010.00000.00000.00010.00030.00040.00040.00050.00040.00040.00030.00010.00000.00000.00020.00040.00050.00060.00060.00060.00050.00040.00020.00000.00000.00020.00040.00060.00070.00070.00070.00060.00040.00020.00000.00000.00020.00050.00060.00

18、070.00080.00070.00060.00050.00020.00000.00000.00020.00040.00060.00070.00070.00070.00060.00040.00020.00000.00000.00020.00040.00050.00060.00060.00060.00050.00040.00020.00000.00000.00010.00030.00040.00040.00050.00040.00040.00030.00010.00000.00000.00010.00010.00020.00020.00020.00020.00020.00010.00010.00

19、000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000t=0.5時(shí)刻的絕對(duì)誤差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00070.00130.00180.00210.00220.00210.00180.00130.00070.00000.00000.00130.00250.00340.00400.00420.00400.00340.00250.00130.00000.00000.00180.00340.00

20、470.00550.00580.00550.00470.00340.00180.00000.00000.00210.00400.00550.00650.00680.00650.00550.00400.00210.00000.00000.00220.00420.00580.00680.00710.00680.00580.00420.00220.00000.00000.00210.00400.00550.00650.00680.00650.00550.00400.00210.00000.00000.00180.00340.00470.00550.00580.00550.00470.00340.00

21、180.00000.00000.00130.00250.00340.00400.00420.00400.00340.00250.00130.00000.00000.00070.00130.00180.00210.00220.00210.00180.00130.00070.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000t=1.0時(shí)刻的絕對(duì)誤差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00160.00

22、310.00430.00510.00530.00510.00430.00310.00160.00000.00000.00310.00590.00820.00960.01010.00960.00820.00590.00310.00000.00000.00430.00820.01130.01320.01390.01320.01130.00820.00430.00000.00000.00510.00960.01320.01560.01640.01560.01320.00960.00510.00000.00000.00530.01010.01390.01640.01720.01640.01390.01

23、010.00530.00000.00000.00510.00960.01320.01560.01640.01560.01320.00960.00510.00000.00000.00430.00820.01130.01320.01390.01320.01130.00820.00430.00000.00000.00310.00590.00820.00960.01010.00960.00820.00590.00310.00000.00000.00160.00310.00430.00510.00530.00510.00430.00310.00160.00000.00000.00000.00000.00000.00000.00000.00