一維波動方程的有限差分法

1、學生實驗報告實驗課程名稱偏微分方程數(shù)值解開課實驗室 數(shù)統(tǒng)學院學院 數(shù)統(tǒng)年級2013 專業(yè)班信計02班學生姓名 學號開課時間 2015 至 2016 學年第 2 學期總成績教師簽名數(shù)學與統(tǒng)計學院制開課學院、實驗室：數(shù)統(tǒng)學院實驗時間：2016年6月20日實驗項目名 稱一維波動方程的有限差分法實驗項目類型驗證演示綜合設(shè)計其他指導教師曾芳成績是一.實驗目的通過該實驗，要求學生掌握求解一維波動方程的有限差分法，并能通過計算機諦言編程 實現(xiàn)。二.實驗內(nèi)容考慮如下的初值問題：22uu22，X 0,1，t 0,2txu x,0 sin x, u x,00（1）u 0,t u 1,t0,t0,21 .在第三部分

2、寫出問題（1）三層顯格式。2 .根據(jù)你寫出的差分格式，編寫有限差分法程序。將所寫程序放到第四部分。3 .取h 0.1,0.1h ,分別將t 0.5,1.0,1.5,2.0時刻的數(shù)值解畫圖顯不。4,該問題的解析解為u x,t cos tsin x,將四個時刻的數(shù)值解的誤差畫圖顯示， 對數(shù)值結(jié)果進行簡單的討論。三.實驗原理、方法（算法）、步驟1、三層顯格式建立由于題中h 0.1,0.1h, x 0,1 ,t 0,2，取 N 10, M 200,故令網(wǎng)比 r 0.1, hXj j h, j 0,1,2,L10,tk k ,k0,1,L ,200,在xj,tk內(nèi)網(wǎng)個點處,利用二階中心差商得到如下格式:

3、k 1 ujk2ujk 1 ujkkuj 1 2ujkuj 1h2o h2略去誤差項得到：k 12 k2 k 2 k k 1uj r uj 1 2 1 r uj r uj 1 uj其中j 1,2,L 9,k 1,2,L ,199,局部截斷誤差為o 2 h2。對于初始條件u x,0 sin x ,建立差分格式為：u0 sinxjsin jh , j 0,1,2,L 10對于初始條件-u x,00 ,利用中心差商，建立差分格式為:t11uj uj20,即 u：二%1, j 0,1,2,L 10對于邊界條件u 0,tu 1,t0,t0,2 ,建立差分格式為:u； uN 0,k 0,1,L ,200(

4、2)(3)(4)(5)(6)將差分格式延拓使k 0為內(nèi)點，代入（3）得到的式子再與（5）聯(lián)立消去可1后整理得至上 1120.20120uj 2r uj 11 r uj 2r uj 1綜上（3）、（4）、（6）、（7）得到三層顯格式如下：（局部截斷誤差為o 2 h2 ）uk 1r2uk1 21 r2ukr2u：1u：1,j1,2,L 9,k1,2,L ,199(8)u0 sinxj sinjh , j0,1,2,L 1011 2 0201 2 0uj2r uj 11 r uj 2 r uj 1, j 1,2,L 9ukuN 0,k 0,1 ,L ,200其中r 0.1。 h四.實驗環(huán)境（所用軟件

5、、硬件等）及實驗數(shù)據(jù)文件Matlab三層顯格式程序如下：%一維波動方程，三層顯格式求解法h=0.1;tau=0.1*h;r=tau/h;N=1/h;M=2/tau;x=0:h:1;t=0:tau:2;u=sin(pi*x);% 計算t=0時刻的u值 u(1,11)=0;for j=2:Nu(2,j)=0.5*rA2*u(1,j+1)+(1-rA2)*u(1,j)+0.5*rA2*u(1,j-1);end%定義x=0邊界上的數(shù)值for k=1:M+1u(k,1)=0；end%定義x=1邊界上的數(shù)值for k=1:M+1u(k,N+1)=0;end%迭代計算開始，差分格式for k=2:Mfor

6、j=2:Nu(k+1,j尸rA2*u(k,j+1)+2*(1-rA2)*u(k,j)+rA2*u(k,j-1)-u(k-1,j);endendu(201,:)=zeros(1,11);%計算k=201行的數(shù)值解 u2(201,11)=0;for j=2:Nu2(201,j)=rA2*u(200,j+1)+2*(1-rA2)*u(200,j)+rA2*u(200,j-1)-u(199,j);endu=u+u2;u=rot90(u,2);% 將矩陣u旋轉(zhuǎn)180度賦值于u%作出圖像x,t=meshgrid(0:0.1:1,0:0.01:2);%劃分網(wǎng)格%作出數(shù)值解的函數(shù)圖像subplot(2,2,1

7、);mesh(x,t,u);title('u(x,t)數(shù)值解的函數(shù)圖像');xlabel('x 變量')；ylabel('t 變量')；zlabel('u 值')；%作出精確解的函數(shù)圖像subplot(2,2,2);u1=cos(pi*t).*sin(pi*x);mesh(x,t,u1);title('u(x,t)精確解的函數(shù)圖像');xlabel('x 變量')；ylabel('t 變量')；zlabel('u 值')；%作出t=0.5 ,1.0, 1.5, 2.0

8、時刻的絕對誤差圖像subplot(2,2,3);wucha=abs(u-u1);x=0:h:1;plot(x,wucha(51,:),'g*-');hold ongrid onplot(x,wucha(101,:),'ro-');hold onplot(x,wucha(151,:),'ks-');hold onplot(x,wucha(201,:),'mp-');title('t=0.5 ,1.0, 1.5, 2.0時刻的絕對誤差函數(shù)圖像');xlabel('x 變量');ylabel('絕

9、對誤差值');legend('t=0.5','t=1.0','t=1.5','t=2.0');%作出t=0.5 ,1.0, 1.5, 2.0時刻的數(shù)值解函數(shù)圖像subplot(2,2,4);x=0:h:1;plot(x,u(51,:),'g*-');hold ongrid onplot(x,u(101,:),'ro-');hold onplot(x,u(151,:),'ks-');hold onplot(x,u(201,:),'mp-');title('

10、;t=0.5 , 1.0, 1.5, 2.0時刻的數(shù)值解函數(shù)圖像');xlabel('x 變量');ylabel('u 值');legend('t=0.5','t=1.0','t=1.5','t=2.0');%當然也可以作出u(x,t)絕對誤差的函數(shù)圖像%mesh(x,t,wucha);%title('u(x,t)絕對誤差的函數(shù)圖像');%xlabel('x 變量')；%ylabel('t 變量')；%zlabel('絕對誤差值'

11、;);五.實驗結(jié)果及實例分析1、u(x,t)在 t=0.5,1.0,1.5,2.0時刻的數(shù)值解、精確解以及絕對誤差表1 u(x,t)在t=0.5,1.0,1.5,2.0 時刻的數(shù)值解時刻tt=0.5,1.0,1.5,2.0時刻的數(shù)值解0-0.005-0.011-0.015-0.018-0.019-0.018-0.015-0.011-0.0050t=0.59352225390-0.309-0.587-0.809-0.951-0.999-0.951-0.809-0.587-0.3090t=1.000700900700t=1.500.00200.00380.00520.00610.00640.006

12、10.00520.00380.00200t=2.00.30900.58780.80900.95111.00000.95110.80900.58780.3090表2 u(x,t)在t=0.5,1.0,1.5,2.0 時刻的精確解時刻tt=0.5,1.0,1.5,2.0時刻的精確解00t=0.50.00000.00000.00000.00000.00000.00000.00000.00000.00000-0.309-0.587-0.809-0.951-1.000-0.951-0.809-0.587-0.3090t=1.000801010800t=1.50.00000.00000.00000.000

13、00.00000.00000.00000.00000.0000t=2.000.30900.58780.80900.95111.00000.95110.80900.58780.30900o表3 u(x,t)在t=0.5,1.0,1.5,2.0時刻的絕對誤差時刻tt=0.5,1.0,1.5,2.0時刻的絕對誤差t=0.5t=1.0t=1.5t=2.00.00590.01130.01550.01820.01920.01820.01550.01130.00590.00000.00000.00010.00010.00010.00010.00010.00000.00000.00200.00380.00520.00610.00640.00610.00520.00380.00200.00000.00000.00000.00000.00000.00000.00000.00000.0000說明：在t=0.5時刻的絕對誤差最大，t=1.5時刻次之，t=1與t=2時刻的絕對誤差均較小，由于r - 0.1 1,該格式穩(wěn)定，由數(shù)值計算得到的矩陣不難看出，數(shù)值