热传导方程的差分格式汇总

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