热传导方程差分格式

热传导方程的差分格式第2页
一维抛物方程的初边值问题
分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题：
.: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)=ef t s in (二x)进行比较
1差分格式形式
2
设空间步长h =1/ N ,时间步长• • 0, T =M •，网比r = • / h .
(1) 向前差分格式
该问题是第二类初边值问题(混合问题)，我们要求出所需次数的偏微商的函数
Eu c2u
u(x,t)，满足方程—a—^, 0 ：：：x ：：：1,和初始条件u(x,0)= sin x , 0 x ：：： 1抚ex
及边值条件u(0,t) =u(1,t) =0, t 0。

已知sin二x在相应区域光滑，并且在x =0,丨与边值相容，使问题有唯一充分光滑的解。

取空间步长h =1/ N，和时间步长-T /M，其中N,M都是正整数。

用两族平行
直线x= j X =( j h0Hj 1 ，和tlNt k =X(k=0,1,|||,M) 将矩形域
G二｛0 — x — 1,t — 0｝分割成矩形网络，网络格节点为(X j,t k)。

以G h表示网格内点集合，
即位于矩形G的网点集合；G h表示闭矩形G的网格集合；丨h=G h-G h是网格界点的集合。

向前差分格式，即
k 1 k k k k
u, -u, u, 1 -2比■ u, 4
- j二a 亠2亠t ( 1)
£h
k 1 k
U j _Uj
[ T
2
k 1 c k 1
a U j 1 -2U j
h 2
k k k
U j 1 _2U j U j 和
]
f j
(3)
0 U j
=(X j ), U 0 = U N = 0.
f i = f(x),
k k
U j j = (X j ), U o = U N = 0
其中，j =1,2，…，N _1,k =1,2,…，M 一1.以r =a ./h 2表示网比。

(1)式可改写成如下：
u ：+ =ru
+(1 -2r)u ： + ru ：」j
此格式为显格式。

其矩阵表达式如下:
1 -2r r
(j 、
U 1
厂j *、
U 1 r
1 -2r
+
U
2
*
U 2*
I-
r 1 -2r
r j
U
N 」
j * U N J
r
1 -2r
丿
j
<UN 丿
j * <UN 丿
(2) 向后差分格式
向后差分格式，即
此种差分格式被称为隐格式。

其矩阵表达式如下：
勺+2r
—r
\ f j 八
U 1
f j 、
U 1 —r
1 +2r
+ r
U 2
9
j U 2
■
-r 1 +2r
-r
j +
U N4
j
U N4
-r 1 + 2r >
j + <U N 丿
j < U N 丿
(3) 六点对称格式
六点差分格式:
k 1 U
j
T
k
— U k 1
k 1 k 1 U
j 1 -纠
U
j 」 二
a
h
(2)
u 0
其中 j =1,2, , N -1,k =1,2, ,M
一1・(2)
式可改写成 -ru k ； (1 2r)u k 1
k 1 _rU j 」
=Uk f j
热传导方程的差分格式
第4页
将(3)式改写成
r k .1 一 、 k i r k 1 r k _
、 k r k
_2出 1 (1 r)U j ―刁山/ =?山
i
(1_r
)U j
2u j-<
f
j
其矩阵表达式如下： 广1+r -r/2
、
-r/2 1 +r
+
j
出
U 2
■
'1 - r r/2
、
r/2
1-r
5 \ u
2 *
—r/2 1 +r -r/2
u
N A
r/2 1-r
r/2
u
N_i
、、
_r 1 +2r /
u
j + 凹丿
r/2 1-2r /
<u N 」
2利用MATLAB 求解问题的过程
对每种差分格式依次取 N=40.,
=1/1600,
=1/ 3200，.=1/6400，用 MATLAB
求解并图形比较数值解与精确解，用表格列出不同剖分时的 L 2
误差。

向前差分格式：
t 705 :
=1/1600 :
*数值解 ——精确解
-1
0 0,1 0.2 0 3 0.4 0.6 0.6 0.7
D.日
0.9 1
=1/3200：
.=1/6400:
*数值解
---- 精确解
05
0.4
02
n £ I
，
$ * N li n ■ I
中
0 0.1 0.2
0 3
0.4 0.6 0.6 0.7 0,0 0.9
1
t =0.1:
=1/3600：
-
、1
*数值解 -- 精确解
/
■
/
- / / / L / /
- 4 / /
-/
/ #
fl 1 11 li \
\ - \ \ \ \ \
06
05
0.4
0 3
0 2
0 1
°0
0.1 0.2 0.3 Q.4 0.6 G.6 0.7 0.0 0.9 1
06
热传导方程的差分格式
第6页
x 10 2.5 | ---- 2 - 15 -
-1 - -
+ *
-1.5 -
*
-
*
显- *
* -
■
9 片 ____ | _____ i ______ I ________ a _____ | ______ | _____ | ______ | ______ i _______ ''0
0.1
0.2
0 3
0.4
0.5
0.6
0.7
0,0
0.9
1
.=1/3200:
*数值解
精确解-
0.25
0.15
0.05
n £ H * N li n ■ I ]
0 0.1 0.2 0.3 0.4 0.6 0.6 0.7 0,0 0.9
1
*数值解
---- 精确解
0 36
■ =1/6400:
0 J
t =0.2
.=1/1600：
=1/3200：
0.36
0 3
0.25
0 2
0.15
0 05
3
10
-1
'3x 1D 13S
0,08
0 06
0 04
0 02
Q ] I I H i I I I I I
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9
1
.=1/6400:
0.14
*数值解
精确解
0,08
0 06
0 04
0 02
n £ H
， i * N k n i j
0 0.1 0.2
0.3 0.4 0.6 0.6 0.7 0,0 0.9
1
向后差分格式:
D.14
0.12
*数值解 ---- 精确解
0.12
t 二0.05 :
;
U 32
O O
m6400
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9
1
05
0.4
Q ] I I H i I I I I ]
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9
1
t =0.1:
.=1/1600:
*数值解 精确解
0.25
0.15
0.05
Q ] I I H i I I I I ]
06
*数值解 ---- 精确解
0.36
=1/3200
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9 1
n £ I
，
$ * N li n ■ i ]
0 0.1
0.2
0.3 0.4 0.6 0.6 0.7 0,0 0.9
1
0.25
0.15
0.05
Q ] I I H il i I I I I ]
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9
1
.=1/6400
*数值解
精确解
0.25
02
0.15
0 05
t =0.2
= 1/1600:
0.36
*数值解 ---- 精确解
0 36
热传导方程的差分格式 第15页
0 0.1
0.2
0 3 0.4 0.6 0.6 0.7 0,0 0.9
1
0,08
0 06
0 04
0 02
Q ] I I H il i I I I I I
0 0.1 0.2 0 3 0.4 0.5 0.6 0.7 0,0 0.9
1
.=1/3200
0,08
0 06
0 04
0 02
n £ H
， $ * N k n i j
■ =1/6400
0.12
*数值解 ---- 精确解
D.1J
0.12
*数值解 ---- 精确解
D 14
0.12
0 1
0.00
0.06
0 04
0 02
°00.1 0.2 0 3 0.4
0.5 0.6 0.7 O.e 0.9 1
六点差分格式：
t =0.05 :
.=1/1600:
Q ]
I
I
H
i
I
I
I
I
]
0 0.1 0.2 0.3 0.4 0.6 0.6 0.7 0,0 0.9
1
=1/3200
06
05
0.4
*数值解 ---- 精确解
热传导方程的差分格式
第18页
n
『 【I
ii n I
N ]
0 0.1
0.2
0.3
Q.4
0.6
G.6 0.7 0.0 0.9 1
0 7
06
05
0.4
0 3
02
1
0.2 0 3 0.4 0.5 0.6 0.7 0,0 0.9
1
.=1/6400
06
05
0.4
*数值解
---- 精确眸
02
t =0.1:
= 1/1600:
0 J
0.15
0.05
Q ] I I H il i I I I I ]
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9
1
.=1/3200
02
0.15
0.05
n £ 【I
， $ * N li n i ]
0 0.1 0.2 0.3
0.4 0.6 0.6 0.7 0,0 0.9
1
■ =1/6400
0.36
0.25
*数值解 ---- 精确解
0 36
0.25
*数值解 ---- 精确解
热传导方程的差分格式 第20页
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9
1
0.15
0.05
Q ] I I H il i I I I I ]
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9
1
t =0.2
.=1/1600:
0,08
0 06
0 04
0.36
0.25
*数值解 ---- 精确解-
0.12
/W /*
*数值解 ---- 精确解
热传导方程的差分格式第21页
0 02
Q ] I I H il i I I I I I =1/3200
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0,0 0.9 1
热传导方程的差分格式 第22页
0,08
0 06
D.14
0.12
0 1
0.08
0.06
0 04
0 02
0.1
.=1/6400
0.14
0.12
*数值解 ---- 精确解
0.2 0.3 0.4 0.5
0 04
0 02
n £I ，$ * N R n ■I ] 0 0.1 0.2 0.3 0.4 0.6 0.6 0.7 0,0 0.9 1
热传导方程的差分格式第24页
L2误差:
3方法总结及分析
本文向前差分格式，向后差分格式及六点差分格式都是使用三对角系数矩阵，计算简单。

根据matlab作，特别明显的是，t =0.05® =1/1600 :时，图像解析解波动特别大，与真
1
解差距很大。

这是因为r =1 一一差分格式不稳定。

根据误差比较，解本题时向后差分格式
2
误差最小。

衡量一个差分格式是否经济实用，有点多方面因素决定，应考虑不同的情况决定。

附件程序
向前差分格式：
fun cti on [ u , err ] = xq( t , t1 )
% t 是时刻值；
% t1 是时间步长；N = 40;
h = 1 / N; x=[h : h : (1-h)]';
r = t1 / (h A2);
u(:,1)=sin(pi * x); %t=0 时刻
R = zeros(N-1 , N-1);
for i = 2 : N-2
R(i , i) = 1 - 2 * r;
R(i , i+1) = r; R(i, i-1) = r;
end
R(1 , 1) = 1 -2 * r; R(1 , 2) = r;
R(N-1 , N-1) = 1 - 2 * r; R(N-1, N-2) = r; k = t / t1;
for i = 1 : k
u(: , i+1) =R * u(: , i); end u=[0 ; u(: , i+1) ; 0]; for i=1 : k for j=1 : N+1 up(j , i)=exp(-pi*pi*t) * sin(pi*(j-1)*h); end end x=0 : h : 1;
plot(x , u , 'b.' , x , up(: , k) ,
'r--' );
legend( ' 数值解' , ' 精确解' ); err=norm(u - up(: , k));
end
向后差分格式：
function [ u , err ] = xh( t , t1 ) N = 40;
h = 1 / N;
x=[h:h:(1-h)]';
r = t1 / (hA2);
u(:,1)=sin(pi * x); %t=0 时刻
R = zeros(N-1 , N-1);
for i = 2 : N-2
R(i ,i) = 1 + 2 * r;
R(i , i+1) = -r;
R(i, i-1) = -r；
end
R(1 , 1) = 1 + 2 * r;
R(1 , 2) = -r;
R(N-1 , N-1) = 1 + 2 * r;
R(N-1, N-2) = -r;
k = t / t1;
for i = 1 : k
u(:,i+1) =inv(R) * u(:,i);
end
u=[0;u(:,i+1);0];
for i=1:k
for j=1:N+1
热传导方程的差分格式第26页
up(j,i)=exp(-pi*pi*t)*si n(pi*(j-1)*h);
end
end
x=0:h:1;
plot(x,u, b' ,x,up(:,k), 'r--');
legend('数值解’，‘精确解');
err= norm(u-up(:,k));
end
六点差分格式：
fun cti on [ u , err ] = ld( t , t1 )
N = 40;
h = 1 / N;
x=[h:h:(1-h)]';
r = t1 / (h A2);
u(:,1)=sin(pi * x); %t=0 时刻
R = zeros(N-1 , N-1);
for i = 2 : N-2
R(i ,i) = 1 + r;
R(i , i+1) = -r/2;
R(i, i-1) = -r/2;
end
R(1 , 1) = 1 + r;
R(1 , 2) = -r/2;
R(N-1 , N-1) = 1 + 2 * r;
R(N-1, N-2) = -r;
for i = 2 : N-2
R1(i ,i) = 1 - r;
R1(i , i+1) = r/2;
R1(i, i-1) = r/2;
end
R1(1 , 1) = 1 - r;
R1(1 , 2) = r/2;
R1(N-1 , N-1) = 1 -2 * r;
R1(N-1, N-2) = r/2;
k = t / t1;
for i = 1 : k
u(:,i+1) =inv(R) *R1* u(:,i);
end
u=[0;u(:,i+1);0];
for i=1:k
for j=1:N+1
up(j,i)=exp(-pi*pi*t)*sin(pi*(j-1)*h);
end
end
x=0:h:1;
plot(x,u, 'b.' ,x,up(:,k), legend( ' 数值
'r--' ); 解' , ' 精确解' ) err=norm(u-up(:,k));
end。