2. 实验内容：
用 -方法导出的隐式格式求如下初边值问题
u u 0 , 0 x 1, 0 t 1 2 t x u ( x , 0 ) sin x , 0 x 1 , u ( 0 , t ) 0 , u (1 , t ) 0 , 0 t 1 ,
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0001 -0.0053 0.3102
0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0009 -0.0706
U
n j 1
2U
n j 2
U
n j 1
( x)

j
( x j ), U
2
0, U
n J
0
以r 为
t /( x )
表示网比。由于本是实验中
( x ) sin x
, 所以上述差分格式简化
rU U
0 j
n 1 j 1
(1 2 r )U
运行结果如下： cita=0 M=10 J=10 x = 1.0e+010 * Columns 1 through 9 0 0.0099 Column 10 0 -2.4999 3.6685 -3.2575 2.0465 -0.9383 0.3102 -0.0706
u = 1.0e+010 * Columns 1 through 9
5
0.0000 0.0000 0 0.0000 0 0.0000 0 0.0000 0 0.0000 0 0.0000 0 -0.0000 0 0.0000 0 -0.0000 0 -0.0001 0 0.0099
0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0001 -0.0023 0.0754 -2.4999
淮
海
工
学
院
实 验 报 告 书
课程名称： 实验名称： 班 姓 日 级 名： 期：
微分方程数值解法 抛物型方程隐格式的程序实现 信息 101
学号： 地点 成绩：
数学实验室
指导教师：
数 理 科 学 系
1
1. 实验目的：
按照实验内容对 -方法导出的几种隐式差分格式编程实现求出数值解； 估计 数值解的误差界；大致比较几种不同格式的优劣；结合格式的相容性、稳定性和 收敛性条件简单分析计算结果。
间步长为 x 线：
x x
j
1/ J
和时间步长为 t
，其中 J ，M 为正整数，用两族平行直
j x , j 1 , 2 , ..., J
平行于 t 轴作竖直线和 平行于 x 轴作水平线,
t t n n t , n 1 , 2 , ..., M
将矩形区域 G 分割成矩形网格，得网格剖分图。 第二步: 差分法的目的是：求初边值问题（2）的解 u ( x , t ) 在节点 ( x j , t n ) 的近 似值 U
0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0012 -0.0505 2.0465
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0004 0.0200 -0.9383
Columns 10 through 11 0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
h = 1.0e+009 * 0.0000 0.0007 -0.7543 1.0667 -0.8874 0.5047 -0.2001 0.0534 -0.0087
0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0031 -0.1067 3.6685
0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0001 -0.0024 0.0887 -3.2575
j
( u ) [(
1 2
) tu
xxt
1 8
1 12
( x) u
2
n xxxx
1 2
]
1 2
j
[
1 24
( t ) u ttt
2
( t) u
2
n xxtt
]
j
O (( t )
3
( x) )
4
利用 Fourier 方法判断差分格式的稳定性条件是： 当0
1 2
，当且仅当 r

1 2
(1 2 )
1
时稳定，当
1 2
1 ，对所有的 r
格式
均稳定。 利用最大值原理可以证明差分格式的收敛条件是： r
1 2 (1 )
1
。
第三步: 有限差分方程的解法。根据问题的规模和计算机的容量和速度，选 取适当的解法。这是数值计算的关键一步，有限差分法解方程的主要计算都集中 在这里。 由于本实验问题的有限差分方程为 -方法导出的隐式差分格式，可以用 Thomas 方法求解。对每个时间层 n 1 需要求解的方程为：
n 0
0, U
n J
0
。
第四步: 根据前面的分析，编制程序，上机计算，求出数值解，必要时画出 解的几何图形，分析观察解的性质。 为了方便比较可以用傅立叶级数方法求出定解问题的解析解为：
u ( x, t) e
t
2
sin x
。
另外，上述差分程组为三对角方程组，可用 Thomas 法（追赶法）求解。 若有三对角方程组：
n j
( j 1, 2 , , J 1; n 1, 2 ,..., M ) 。为此，需要构造逼近微分方程定解问
题的差分格式。采用 -方法，构造如下差分格式
U
n 1 j
U
n jt U0 j NhomakorabeaU
n 1 j 1
2U
n 1 j 2
U
n 1 j 1
( x)
n 0
(1 )
6
>>
cita=0.5 M=10 J=10 x= 1.0e-003 * Columns 1 through 9 0 -0.0037 Column 10 0 -0.9215 0.2009 0.2907 0.0655 -0.0568 -0.0596 -0.0252
7
u= Columns 1 through 9 0.3090 0.3090 0 0.1041 0 0.0305 0 0.0063 0 0.0017 0 0.0006 0 -0.0000 0 0.0001 0 -0.0000 0 0.0000 0 -0.0000 0.5878 0.1335 0.0163 0.0139 -0.0029 0.0037 -0.0023 0.0019 -0.0014 0.0011 -0.0009 0.8090 0.2338 0.0425 0.0174 0.0012 0.0017 -0.0002 0.0001 0.0001 -0.0002 0.0002 0.9511 0.2981 0.0664 0.0183 0.0046 0.0005 0.0009 -0.0005 0.0005 -0.0004 0.0003 1.0000 0.3250 0.0813 0.0185 0.0061 0.0004 0.0009 -0.0003 0.0002 -0.0001 0.0001 0.9511 0.3148 0.0846 0.0178 0.0059 0.0008 0.0005 0.0000 0.0000 0.0001 -0.0001 0.8090 0.2704 0.0760 0.0156 0.0049 0.0010 0.0002 0.0002 -0.0001 0.0001 -0.0001 0.5878 0.1975 0.0571 0.0117 0.0034 0.0009 0.0000 0.0002 -0.0001 0.0001 -0.0000
，
3
xn fn x i f i e i x i 1 , i n 1 , n 2 , ,1
其中 x 1 ,
x 2 , ..., x n
为所求三对角方程组的解。
四． 实验数据记录及分析（或程序及运行结果） ：
程序如下：
clear cita=input('cita='); M=input('M='); J=input('J='); u=zeros(J+1,M+1); d=zeros(1,J); e=zeros(1,J-1); f=zeros(1,J); x=zeros(1,J); uu=zeros(1,J); for n=1:M u(n,1)=0; u(n,M)=0; end for j=1:J u(1,j)=sin(j*pi/J); end tt=1/M; xx=1/J; r=tt/(xx*xx); for n=1:M for j=2:J-1; d(1,j)=(1-cita)*r*u(n,j+1)+(1-2*(1-cita)*r)*u(n,j)+(1-cita)*r*u(n,j-1) ; a=cita*r; b=1+2*cita*r; c=cita*r; end e(1,2)=c/b; f(1,2)=d(1,2)/b; for i=3:J-2 e(1,i)=c/(b-e(1,i-1)*a); end for i=3:J-1