熱傳導方程的差分格式匯總

1、熱傳導方程的差分格式 上機實習報告 二零一四年五月熱傳導方程的差分格式 第 1 頁 一維拋物方程的初邊值問題 分別用向前差分格式、向后差分格式、六點對稱格式，求解下列問題: .:u ；：2u a 2， 0 ： x ： 1, .:t ；x u(x,0) =sin 二 x, 0 ： x : 1 u(O,t) =u(1,t) =0, t 0 2 在 t =0.05,0.1 和 0.2 時刻的數值解，并與解析解 u(x,t)=efts in (二 x)進行比較。 1 差分格式形式 2 設空間步長h =1/ N , ,時間步長0, T = M ，網比r z/h . . (1) 向前差分格式 向前差分格式

2、，即 k: H k k k k Uj _Uj u j 1 _2Uj Uj j L 乂 U J J ( 1 1) T h fi = f (xi )， uj = j = (Xj),u0 二 uN = 0 其中,j =1,2, , N -1, k =1,2, ,M -1.以r /h2表示網比。(1 1)式可改寫成如下: k 1 k k k uj ru j 1 (1 -2r)uj ruj 4 fj 此格式為顯格式。 其矩陣表達式如下： 勺-2r r (j 、 u1 j + u1 r 1 -2r + j u2 j + u2 9 r 1 -2r r j uN -4 j + uN 4 r 1 -2r j I

3、 uN丿 j + u1 j卅 u2 R- r j 、 u1 j u2 * -r 1 +2r -r j卅 uN d. uN -r 1 +2r ? j卅 j I uN丿 (3) 六點對稱格式 六點差分格式: uj = j = (Xj), u0 = uN 二 0- 將(3 3)式改寫成 r k 1 k 1 r k 1 r k k r k uj 1 (1 r)uj ujj U. 1 (1 r)uj U. fj 2 2 2 2 其矩陣表達式如下： 1 +r r/2 f j八 u1 *1 - r r/2 f j 、 u1 -r/2 1 + r + j + u2 a r/2 1 -r + u2 -r/2 1

4、+r -r/2 j出 uN J r/2 1 -r r/2 j uN J -r 1 +2r j j + luN丿 r/2 1-2r j TLS 11 through 20 h Q015 -L0556 I. OHB -L 1310 1.1524 -1,1633 1.1539 -1.1549 Cluns 2 through 30 1. 0775 -1.0376 0. 9922 -0. 9422 0.8B35 -0. 8321 0. 7736 -0.7139 CDIUJUIS 31 through 40 0.5337 -0. 4749 0.41：3 -0. 3B0d 0.3000 -0. 2S25 0

5、. 2002 -0.1491 CDIULUI 0. B613 L 1370 Q. 6537 0, 0000 -0.936 -I, 110 -Q. 593 -0. 0鵡 熱傳導方程的差分格式 第 8 頁 (2) 向后差分格式 back (0.025, 1/6400, 0. 05 Cuxriit plot held ans 二 1.0e-fi3 * ClunuiiE 1 through S 0 -0.0303 -0. 0605 -0.0903 -0. 1195 -0.L4S0 -0.1756 -0. 2021 -0, 2273 匚 clunuis 10 through. L3 -0.2512 -0

6、.2735 -a 2941 -0.312B -ft. 32M -0.3446 -山 3573 -0+3678 -0.3761 Colums 19 through 27 -0.3820 -0.3856 -0.3863 -CL 3656 -CL 3820 -0.3761 二山 3678 二血 3673 二 6 3446 Columns 28 through 36 -0.32S3 -0.3129 -0.2941 -0.2735 -0.2S12 -0.2273 -0.2021 仇 1756 仇 1480 Colums 37 throBih 41 -0. 1155 -0.0903 -0. 06G5 -0

7、.0303 0.0000 熱傳導方程的差分格式 第 9 頁 back (0.025, 1/3200.0. ! Current plot held ans = L 0?-03 * Columns ) through 9 0 -0. 0593 -D. 1182 -Q.1763 -0. 2334 -0.2891 -0. 3429 -0. 3947 -0.4440 Columns 10 throush 18 -0.4906 -0.5341 -0.5744 Or 61 JI -0. 6440 -Q. 6730 7 8S78 -0. 7184 -(J-. 7345 Coltnms 19 th工口隅h 27

8、 -0.7460 - th 7530 -0. 7553 -M 7530 -0. 7460 -0. 7345 -0* 71S4 一D 6978 -0.6730 Columns 28 through 36 -0 . 6440 -0. 6111 -0.6744 -0.5341 -0. 4906 T 4440 -0. 3947 P 3429 -0.28&1 Columns 37 through 41 -0.2334 -0.1763 -0.1182 -0. 0S93 0. 0000 熱傳導方程的差分格式 第 10 頁 0000 bKkOS 1/1600,02) Current plot hld

9、ans = L 0e-03 * Coluara 1 through & 0 -0-0773 -a. 1542 -0. 2301 -0 Columns 10 through 16 &40 2 -0.6570 -0.7496 -0. 7975 -c CDI ULBTIS 19 through. 27 -0. 9736 -0. 9827 -0, A958 -0. 9327 -0 CDIUBTLS 2S through 36 7 強05 -0.775 7 7496 E0.6B70 -0 Coluans 37 through 4) -0. 3046 -a. 2301 -0.1542 -0

10、. 0773 03O4 0.1) Current plot held arts = L 0e*03 * 匚olumns 1 through 10 0 00148 -0. 0295 -0. 0441 -0.0583 -0.0723 CL 0857 -0. 0987 -0.1110 -0a 1226 11 thx口蟲 h 20 -0.1335 -0. 1436 -0.15Z7 -0, 1610 -0.1682 -0.1744 - (L 1796 -0. 1836 -0.1865 0.1382 匚alumns 21 through 30 -Q.18B8 7 1882 -OH im -0. 1836

11、7. 1796 -Q.1744 -0+ 1682 -0. 1510 -0.1527 0el436 Columns 31 through 40 -0.1335 -0 1226 -0.1110 -0.G907 -0.0857 -0. 0723 -Ox 0583 -0. 0441 -0.0295 -0.O14B 熱傳導方程的差分格式 第 13 頁 0.0000 Column 41熱傳導方程的差分格式 第 14 頁 0.0000 31x(0.025v L/1B00, 0.2) Current plot held I ans三 L 0e-03 * Columns 1 through 10 0 -0,01

12、10 -0” 021& -0.0327 -0. 0433 CL 0S36 -0.0636 -0.0732 -0.0S24 -C.0S1 Colunns 11 throush 20 -0.0991 -0.1060 -0, 1134 -0, 1195 -ft. 1249 -i 1ZB5 -0. 1333 -0. 1353 -0. 1384 -山 139 Colunns 21 through 30 -OH 1401 -0+ 139? -0+1384 -QU363 -0 1333 -也 129E -0. 1249 -0. UB5 -0. 1134 -山 1C6 Colunns 31 thruu

13、sh 40 -0. 0991 -0.0910 -0+ 0324 -0,0732 -0. Q63S -Q. 0&36 -0, M33 -Q. 0327 -0. 021B -0.01 J 熱傳導方程的差分格式 第 15 頁 0.0000 Column A1熱傳導方程的差分格式 第 16 頁 3 方法總結及分析 六點對稱方法克服了向前差分格式與向后差分格式方法誤差偏大的缺點， 更接近 準確數值 本文向前差分格式，向后差分格式及六點差分格式都是使用三對角系 數矩陣， 計算簡單。 根據 matlab 作， 特別明顯的是，t = o.O5：. =1/1600： 時，圖像解析解波動特別大，與真解差距

14、很大。 附件程序 (1) 向前差分格式 源程序： fun cti on e=forward(h,t,T) N=1/h; x=zeros(1,N+1); for i=1:N x(i+1)=i*h; end u=s in( pi.*x); r=t/(h*h); for j=1:T/t; u=for_climb(u,r,t); end plot(x,u, o) hold u_xt=exp(-pi*pi*T)*si n(pi.*x); plot(x,u_xt, r) e=u_xt-u; fun cti on b=for_climb(a,r,t) L=le ngth(a); b=zeros(1,L);

15、for i=1:L-2 b(i+1)=r*a(i+2)+(1-2*r)*a(i+1)+r*a(i); end 熱傳導方程的差分格式 第 17 頁 (2) 向后差分格式 源程序： fun cti on e=back(dx,dt,T) M=1/dx; N=T/dt; u1=zeros(1,M+1); x=1:M-1*dx; u1(2:M)=si n( pi.*x); r=dt/dx/dx; r2=1+2*r; for i=1:M-1 A(i,i)=r2; if i1 A(i-1,i)=-r; A(i,i-1)=-r; end end u=zeros(N+1,M+1); u(N+1,:)=u1; f

16、or k=1:N b=u(N+2-k,2:M); u(N+1-k,2:M )=in v(A)*b; end uT=u(1,:); x1=0,x,1; plot(x1,uT, o) hold u_xt=exp(-pi*pi*T)*si n( pi.*x1); plot(x1,u_xt, r) e=u_xt-uT; (3) 六點對稱格式: 源程序： fun cti on e=six(dx,dt,T) M=1/dx; N=T/dt; u1=zeros(1,M+1); x=1:M-1*dx; u1(2:M)=si n( pi*x); r=dt/dx/dx;r2=2+2*r;r3=2-2*r; for i=1:M-1 A(i,i)=r2; if i1 A(i-1,i)=-r; A(i,i-1)=-r; end end for i=1:M-1 B(i