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

二維波動(dòng)方程的差分士Documentserialnumber【UU89WT-UU98YT-UU8CB-UUUT-UUT108】學(xué)生實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)課程名稱 偏微分方程數(shù)值解開課實(shí)驗(yàn)室 數(shù)統(tǒng)學(xué)院學(xué)院數(shù)統(tǒng)年級(jí)2013專業(yè)班信計(jì)02班學(xué)生姓名學(xué)號(hào)開課時(shí)間2015至2016學(xué)年第2學(xué)期總成績教師簽名數(shù)學(xué)與統(tǒng)計(jì)學(xué)院制實(shí)驗(yàn)項(xiàng)目名稱二維波動(dòng)方程的有限差分法實(shí)驗(yàn)項(xiàng)目類型驗(yàn)證演示綜合設(shè)計(jì)其他指導(dǎo)教師曾芳成績是開課學(xué)院、實(shí)驗(yàn)室：數(shù)統(tǒng)學(xué)院實(shí)驗(yàn)時(shí)間：2016年6月20日一.實(shí)驗(yàn)?zāi)康耐ㄟ^該實(shí)驗(yàn)，要求學(xué)生掌握求解二維波動(dòng)方程的有限差分法，并能通過計(jì)算機(jī)語言編程實(shí)現(xiàn)。二.實(shí)驗(yàn)內(nèi)容考慮如下的初值問題:d2u d2ud2u 二 1 d12 d.xd2u d2ud2u 二 1 d12 d.x2ey2(x,y,0)=sin兀xsin兀y,—udt(x,y,0(x,y,0)=sin兀xsin兀y,—udt(x,y,t)=0,(x,y)￡SQ,tg[0,1.4].在第三部分寫出問題（1）三層顯格式。.根據(jù)你寫出的差分格式，編寫有限差分法程序。將所寫程序放到第四部分。.取h=0.1j=0.1h，分別將t=0.5,1.0,1.4時(shí)刻的數(shù)值解畫圖顯示。.該問題的解析解為u（x,y,t）=cos<2ktsin兀xsin兀y，將四個(gè)時(shí)刻的數(shù)值解的誤差畫圖顯示，對數(shù)值結(jié)果進(jìn)行簡單的討論。三.實(shí)驗(yàn)原理、方法（算法）、步網(wǎng)格劃分h=0.1/=0.1h,故N=：二10,M=1.4二140,x.二ih,y.二jh,i,j=0,1,,10,t=比，k=0,1,,140。在內(nèi)網(wǎng)點(diǎn)(x,y,t)，利用二階中心差商，對(1)建立差分格式： 1」'uk+1——2uk+uk_1uk——2uk+uk uk ——2uk+uk—W ij——W—=i+1,j W i——1j+——W+1 W W——tT2 h2 h2（2）整理得到:

UkUk一Uk-1uk+1=r25k+Uk+Uk+Uki,j i+1,j i-1,j i,j+1 i,j-1其中，i,j=1,2,,9,k=1,2,,139，網(wǎng)比r二=0.1,局部截?cái)嗾`差為oQ2+h2）。h考慮邊界條件u（x,y,t）=0,（x,y）wdQ,t￡【0,1.4],差分格式為：…Uk=.Uk=Uk=Uk=0,k=0,1,,1400,0 0,N N,0 N,N考慮初始條件u（x考慮初始條件u（x,y,0）=sin兀xsin兀y，差分格式為:u0=sin（兀x）sinGy）=sin（兀ih）sin（兀jh）,i,j,0,1,,10（5）i，j考慮初始條件ut（x,y,0）=0,（x,y）式0,1?，利用二階差商近似:（7）將（6）得到的結(jié)果u1二（7）將（6）得到的結(jié)果u1二u-1代入i,j/i,j1（u1二—r2W0+u0i,j2 i+1,j i-1,j（7）中，整理得到:+u0i,j+1+u0）+G-2r2）u0i,j-1 i,jU二"一：=0,i,j=0,1,,10

2t（6）設(shè)k二0時(shí)刻的點(diǎn)為內(nèi)點(diǎn)，則滿足差分格式（2）），代入上式得到:

u1=r2（u0+u0+u0+u0）+（2-4r2）u0-u-1i,j i+1,j i-1,j i,j+1 i,j-1（8）（8）綜上（2）、（4）、（5）、（8）得到三層顯格式的差分格式為:

uk+1=r2（uk+Uk+Uk+Uk）+（2-4r2）uk-Uk-1i,j i+1,j i-1,j i,j+1 i,j-1 i,j i,ji,j=1,2,,9,k:1,2,,139Uk—Uk—Uk—Uk—0,k—0,1,,140U0廣sin（兀x$nCj—sin（兀ih）sin（兀jh）,i,j—0,1,,10u1——r2（u0+u0+u0+u0）+G-2r2）U0,i,j—0,1,,10、i,j2 i+1,j i-1,j i,j+1 i,j-1 i,j（9）其中r90」，局部截?cái)嗾`差為oQ2+h2）。四.實(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二；tau=*h;%定義步長r=tau/h;%網(wǎng)比[x,y,t]=meshgrid(0:h:1,0:h:1,0:tau:;%空間網(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ì)算fori=2:10forj=2:10u(i,j,2)=*r^2*(u(i+1,j,1)+u(iT,j,1)+u(i,j+1,1)+u(i,jT,1))+(1-2*r^2)*u(i,j,1);u(11,:,2)=0;u(:,11,2)=0;endendu2=u(:,:,2);%第3-141層網(wǎng)點(diǎn)計(jì)算fork=2:140fori=2:10forj=2:10u(i,j,k+1)=r^2*(u(i+1,j,k)+u(iT,j,k)+u(i,j+1,k)+u(i,j-1,k))+(2-4*r^2)*u(i,j,k)-u(i,j,k-1);u(11,:,k+1)=0;u(:,11,k+1)=0;endendend%%%%%%%%%%%%%%%%%%%%%%%%%%結(jié)果分析與作圖％%%%%%%%%%%%%%%%%%%%%%%%%%%%wucha=abs(u-uu);%求絕對誤差矩陣11x11x141wucha1二wucha(:,:,11);%計(jì)算t二時(shí)刻的絕對誤差矩陣11x11wucha2二wucha(:,:,51);%計(jì)算t二時(shí)刻的絕對誤差矩陣11x11wucha3二wucha(:,:,101);%計(jì)算t二時(shí)刻的絕對誤差矩陣11x11wucha4二wucha（:,:,141）;%計(jì)算t二時(shí)刻的絕對誤差矩陣11x11x0=0:h:1;y0=0:h:1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%誤差分析％%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%作1二時(shí)刻的絕對誤差圖subplot（2,2,1）;mesh（x0,y0,wucha1）;title（'t二時(shí)刻的絕對誤差'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'絕對誤差值'）；%作1二時(shí)刻的絕對誤差圖subplot（2,2,2）;mesh（x0,y0,wucha2）;title（'t二時(shí)刻的絕對誤差'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'絕對誤差值'）；%作1二時(shí)刻的絕對誤差圖subplot（2,2,3）;mesh（x0,y0,wucha3）;title（'t二時(shí)刻的絕對誤差'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'絕對誤差值'）；%作1二時(shí)刻的絕對誤差圖subplot（2,2,4）;mesh（x0,y0,wucha4）;title（'t二時(shí)刻的絕對誤差'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'絕對誤差值'）；%%%%%%%%%%%%%%%%%%%%%%%%%%四個(gè)時(shí)刻數(shù)值解、精確解％%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%作弋二、時(shí)刻的數(shù)值解與精確解subplot（2,2,1）;mesh（x0,y0,u（:,:,11））;%作t=時(shí)刻的數(shù)值解title（'t=時(shí)刻的數(shù)值解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）；subplot（2,2,2）;mesh（x0,y0,uu（:,:,11））;%作t=時(shí)刻的精確解title（'t二時(shí)刻的精確解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）;%%作t二時(shí)刻的數(shù)值解與精確解subplot（2,2,3）;mesh（x0,y0,u（:,:,51））;%作t=時(shí)刻的數(shù)值解title（'t=時(shí)刻的數(shù)值解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）;subplot（2,2,4）;mesh（x0,y0,uu（:,:,51））;%作t=時(shí)刻的精確解title（'t二時(shí)刻的精確解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%分別復(fù)制粘貼運(yùn)行％%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%作弋二、時(shí)刻的數(shù)值解與精確解subplot（2,2,1）;mesh（x0,y0,u（:,:,101））;%作t=時(shí)刻的數(shù)值解title（'t=時(shí)刻的數(shù)值解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）;subplot（2,2,2）;mesh（x0,y0,uu（:,:,101））;%作t=時(shí)刻的精確解title（'t二時(shí)刻的精確解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）;%%作t二時(shí)刻的數(shù)值解與精確解subplot（2,2,3）;mesh（x0,y0,u（:,:,141））;%作t=時(shí)刻的數(shù)值解title（'t=時(shí)刻的數(shù)值解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）;subplot（2,2,4）;mesh（x0,y0,uu（:,:,141））;%作t=時(shí)刻的精確解title（'t二時(shí)刻的精確解'）；xlabel（'x變量'）；ylabel（'y變量'）；zlabel（'u值'）;五.實(shí)驗(yàn)結(jié)果及實(shí)例分析圖2t二、時(shí)刻的數(shù)值解、精確解1、圖2t二、時(shí)刻的數(shù)值解、精確解1、t=0.1、0.5、1.0、1.4時(shí)刻的數(shù)值解與精確解圖圖1t二、時(shí)刻的數(shù)值解、精確解注：上兩圖為四個(gè)時(shí)刻的數(shù)值解與精確解，r=%r