熱傳導方程的差分格式

1、一維拋物方程的初邊值問題分別用向前差分格式、向后差分格式、六點對稱格式，求解下列問題:dudid2u0<x<l,u(x, 0) = sin 7tx. 0< x< 1w(0,/) = w(l,/) = 0,/ > 0在t = 0.05,0.1和0.2時刻的數(shù)值解，并與解析解u(x, t)=護sin()進行比較。1差分格式形式設空i'可步長h = ( n,時i'可步長r>0, t = mt,網(wǎng)比r = t/h2.(1)向前差分格式該問題是第二類初邊值問題(混合問題)，我們要求出所需次數(shù)的偏微商的函數(shù)如,滿足方礙m器,0<x< 1,和初

2、始條件u(x,0) = sinxr, 0<x< 1,及邊值條件以0,/) = 心)= 0,/ >0 o已知sin/rx在相應區(qū)域光滑，并且在x = 0j與邊值相容，使問題有唯一充分光滑的 解。取空間步長h = in ,和時間步長t = t/m ,其中n,m都是正整數(shù)。用兩族平行 直 線 x = xj = jh(j = 0,，n) 和 t = tk=lcr(k ",，m) 將 矩形域g = 0<x<l,z>0分割成矩形網(wǎng)絡，網(wǎng)絡格節(jié)點為w 以表示網(wǎng)格內(nèi)點集合, 即位于矩形g的網(wǎng)點集合；gh表示閉矩形石的網(wǎng)格集合；rh =gh-gh是網(wǎng)格界點的集 合。

3、向前差分格式，即fi = /u,)，u(j = (pj =(p(xj), wq =u； = 0其中,) = 1,2,，n_1,e = 1,2,，m_1以尸=m/臚表示網(wǎng)比。(1)式可改寫成如下:r+1wy+i = rukj+i + (1 - 2廠)iij + rukh + tfj此格式為顯格式。其矩陣表達式如下：<l-2rrl-2r(j、njuiw,v-i“n-l7出丿<)r l-2r rr l-2r（2）向后差分格式向后差分格式，即ukk+l 。知1 i -/+1uj uj _/屮一弘 + 勺t , r - 一+力,h2wj = (pj =(p(xj),uq =wn = 

4、6;，其中丿=1,2,，n =,m1. （2）式可改寫成一皿;:;+ (1 + 2廠)uj+ -= uj + tfj此種差分格式被稱為隱格式。其矩陣表達式如下：<l + 2r -r、缶田、(j、-r1+2廠2+，-r 1 + 2廠-r 必一】<一廠1 + 2乙k wiv >（3）六點對稱格式六點差分格式:%""t 2喑-2咁+屹;"+廠2與+屹打jh2h2心=<pj =(p(x),u；=說=0.將（3）式改寫成_ 彳 w./+l +（1 + r）w./+1 - f w；-l =彳 0+1 +（1 _ 廠）uj + 彳 0_1 + 戲j其矩陣

5、表達式如下：l + r-r/2<l-rr/2(j y-r/2l + rr/21-r-r/2l + r -r/2r/21-rr/2un-l-r 14- 2廣丿才丿kr/21一2乙k un )2利用matlab求解問題的過程對每種差分格式依次取n = 40，el/1600, r=l/3200, r=l/6400,用 matlab 求解并圖形比較數(shù)值解與精確解，用表格列岀不同剖分時的厶彳誤差。向前差分格式：t = 0.05:r=l/1600：x 10x 102131 i i i i i i i i i數(shù)值解精確解2 - - 0 1 -1 - . - -2 - - 二；iiiiiiiii&quo

6、t;00.10.20.30.40.50.60.70.80.91r = 1/3200：0.70.60.50.40.30.2數(shù)值解精確解0 彳 i1iiii'00.10.20.30.40.50.60.70.80.91r = 1/6400：0.1:r = 1/3600：ii 數(shù)值解精確解0.50.60.70.80.91x 10592.5 |2 -1.5 -10.5 -0-0.5 -.-1.5 -2 - t 2511100.10.20.30.4r = 1/3200:0.30.250.20.0521%ij.2 0.30.40.50.60.70.80.913注t = q.27 = 1/1600：數(shù)

7、值解精確解-0.350.15數(shù)值解精確解0-1-2 j iiiiiiiii'00.10.20.30.40.50.60.70.80.91r = 1/6400：向后差分格式:t 0.05:"1/1600:r = 1/3200r = 0.1:21/1600:0.40.350.30.250.20.150.10.05000.10.20.30.40.50.60.70.80.910.30.250.20.05%0.20.30.40.50.6070.80.910.350.15數(shù)值解精確解t = 0.2"1/1600:r = 1/3200六點差分格式：t = 0.05:7 = 1/16

8、00:r = 1/64000.1:“1/1600:0.30.250.20.05%數(shù)值解精確解0.20.30.40.50.6070.80.90.350.15r = 1/3200r = 0.221/1600:數(shù)值解精確解0.140.120.10.080.060.040.02000.10.20.30.40.50.60.70.80.910.140.12數(shù)值解精確解0.080.060.040.020 iiiiilli00.10.20.30.40.50.60.70.80.910.140.12”一、 /w/ >數(shù)值解精確解0.080.060.040.02q iiiiii11i00.10.20.30.4

9、0.50.60.70.80.91厶2誤差:3 1向前差分格式00 «ia8+2od e1x o - t20.-t向后差分格式2 8 o 2 o 3 5 o 211 o-i2 o-t六點差分格式052勿 3 o 35 o o 4 7 o 75 7o180.8/1xo-t4 oo120.o4 oo60o 44 40.09012 o-t2 9o040.92 8 o 8 3oz- 11 o o 53方法總結及分析本文向前差分格式，向后差分格式及六點差分格式都是使用三對角系數(shù)矩陣，計算簡單。 根據(jù)matlab作，特別明顯的是，/二0.05" = 1 /1600:時，圖像解析解波動特別

10、大，與真 解差距很大。這是因為r = l a丄差分格式不穩(wěn)定。根據(jù)誤差比較，解本題時向后差分格式2誤差最小。衡量一個差分格式是否經(jīng)濟實用，有點多方面因素決定，應考慮不同的情況決定。附件程序向前差分格式：function u , err = xq( t , tl ) % t是時刻值；% tl是時間步長；n = 40;h = 1 / n;x=h : h :(1-h) *;r = tl / (ha2);u (: a 1) =sin (pi * x) ; %t=0時亥9r = zeros (n-l , n-l);for i = 2 : n-2r (i , i) = 1 - 2 * r;r(i , i+

11、1) = r;r(i, i-1) = r;endr(1 ,1) = 1 -2 * r;r(1 z 2) = r;r(n-1 , n-l) = 1 - 2 * r;r(n-1, n-2) = r;k = t / tl;for i = 1 : ku(: , i + 1) =r * u (: , i);endu=0 ; u(: , i + 1) ; 0;for i=l : kfor j=l : n+lup ( j , i)=exp(-pi*pi*t)* sin (pi* ( j-1)*h);endendx=0 : h : 1;plot (x , u , 1b. 1, x , up(: , k) ,

12、1r-1); legend ('數(shù)值解','精確解');err=norm(u 一 up(: , k);end向后差分格式:function u , err = xh( t , tl ) n = 40;h = 1 / n;x=h:h: (1-h) 1;r = tl / (ha2);u (:z1)=sin(pi * x);%t=0時刻r = zeros (n-l , n-l);for i = 2 : n-2r(i ,i) = 1 + 2 * r;r(i , i+1) = -r;r(i/ i-1) = -r；endr(1 z 1) = 1 + 2 * r;r(1 /

13、2) = -r;r(n-1 , n-l) = 1 + 2 * r;r(n-1, n-2) = -r;k = t / tl;for i = 1 : ku(:, i + 1) =inv(r)* u (:,i);endu=0;u(:,i + 1);0;for i=l:kfor j=l:n+lup(jz i)=exp(-pi*pi*t)*sin(pi*(j-1)*h);endendx=0:h:1;plot(x,u, fb. 1,x,up (:,k)/ 1r-1);legend (*數(shù)值解j，精確解j;err=norm(u-up (:,k);end六點差分格式：function u , err = id

14、( t , tl )n = 40;h = 1 / n;x=h:h:(1-h) *;r = tl / (ha2);u (: , 1) =sin (pi * x) ; %t=0時刻r = zeros(n-l , n-l);for i = 2 : n-2r(i ,i) = 1 + r;r(i , i+1) = -r/2;r(i, i-1) = -r/2;endr(1 , 1) = 1 + r;r(1 , 2) = -r/2;r(n-1 , n-l) = 1 + 2 * r;r(n-l, n-2) = -r;for i = 2 : n-2rl(i ,i) = 1 - r;r1(i , i+1) = r/2;r1(i, i-1) = r/2;endrl(l , 1) = 1 - r;r1 d,2)= r/2;r1 (n-l . n-l) = 1 -2 * r;r1(n-l, n-2) = r/2;k = t / tl;