抛物型方程的古典显格式
微分方程数值解
计算科学系 杨韧
抛物型方程的古典显格式
第二章 抛物型方程的差分格式
模型 长度为 L 的绝缘杆上的一维热流
x = 0
u(0,t)=c1 x =L x u(L,t)=c
2
杆
绝缘体
抛物型方程的古典显格式
热传导方程在时间 t 和位置 x 处的温度 u (x,t)
表示为
初始温度分布为
杆端点的边界值为
其中 k 是导热率系数,σ 是热量,ρ 是杆的密度。
?????
?
??
?
? tLx
t
txu
x
txuk 0,0),(),(
2
2
??
.0),()0,( Lxxxu ????
????
????
tctLu
tctu
0,),(
0,),0(
2
1
抛物型方程的古典显格式
§ 2.1 差分格式建立的基础
一、网格
将求解区域 Ω 分割成 M?N个小矩形，长宽分
别为
⊿ x=h
⊿ t=k
h k
0 x0 x1 … xM x
t
tN
┆
t0
抛物型方程的古典显格式
二、差商
一阶偏导数 的向前、向后、中心差商为
二阶偏导数 的中心差商为
n
mx
u ?
?
??
?
?
?
?
h
uu
h
uu
h
uu nmnmnmnmnmnm
2
,,1111 ???? ???
n
mx
u
???
?
???
?
?
?
2
2
2
11 2
h
uuu nmnmnm ?? ??
抛物型方程的古典显格式
用 (xm,tn)处的一阶向前差商近似代替微商 ux,
即
考虑 u(x,t) 的 Taylor 展式
)5.2(1
h
uu
x
u nmnmn
m
??
?
?
??
?
?
?
? ?
????
?
?
???
?
?
??
???
?
???
?
?
??
???
?
???
?
?
???
?
n
m
n
m
n
m
nmnm x
uh
x
uh
x
uhtxutxu
3
33
2
22
1 !3!2!1),(),(
????
?
?
???
?
?
??
???
?
???
?
?
??
???
?
???
?
?
???
?
n
m
n
m
n
m
nmnm x
uh
x
uh
x
uhtxutxu
3
33
2
22
1 !3!2!1),(),(
抛物型方程的古典显格式
截断误差
)(
!2
!3!2
!3!2!1
12
2
3
32
2
2
3
33
2
22
1
??
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
mm
t
x
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
xxx
x
uh
x
uh
x
uh
hu
x
uh
x
uh
x
uh
u
x
u
h
uu
x
u
E
n
?
?
抛物型方程的古典显格式
用 (xm,tn)处的一阶中心差商近似代替微商 ux,即
截断误差为
)7.2(
2
11
h
uu
x
u nmnmn
m
?? ???
?
?
?
?
?
?
?
n
t
x
n
m
n
m
n
m
n
m
n
m
n
m
x
uh
x
uh
x
uh
h
uu
x
u
E
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
3
32
5
54
3
32
11
!3!5!3
2
?
抛物型方程的古典显格式
三、算子
为 x 方向偏导数算子
Tx为 x 方向位移算子
μ x 为 x 方向平均算子
x
D x
?
??
n
m
n
mx
n
m
n
mx uuTuuT 1
1
1,?
?
? ??
? ?nmnmnmx uuu
2
1
2
12
1
?? ???
抛物型方程的古典显格式
x 方向差分算子
前差算子
后差算子
中心差算子
（ I为恒等算子）
nmnmnmxx uuu ?? ? 1,??
nmnmnmxx uuu 1,?????
n
m
n
m
n
mxx uuu 2121,?? ????
n
mx
n
mxx
n
m
n
m
n
m
n
m
n
m
uhDuD
h
D
h
I
x
uh
x
uh
x
uh
uu
)e x p (
!2!1
!3!2!1
2
2
3
33
2
22
1
??
?
?
?
?
?
????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
抛物型方程的古典显格式
§ 2.2 显式差分格式
差分方法：用差商代替微商，在网格节点上
求出微分方程解的近似值的一种方法。
一、一维常系数热传导方程的古典显式格式
考虑一维热传导方程 (2.26)
在网格节点
2
2
x
u
t
u
?
??
?
?
m– 1,n m,n m+1,n
m,n+1
抛物型方程的古典显格式
的差分近似
取 r = k / h2 为步长比，得显式向前差分方程
(2.29)
(2.29)也称为解热传导方程 (2.26)的古典显式格
式。截断误差为
2
11
1 2
h
uuu
k
uu nmnmnmnmnm ??? ????
n
m
n
m
n
m
n
m rUUrrUU 11
1 )21(
??
? ????
).( 2hkO ?
抛物型方程的古典显格式
例 1 用古典显式格式求解抛物型方程
初始条件为
边界条件为
取步长 ⊿x = h = 0.2,⊿t = k = 0.02 。
20.00,0)(),1(
20.00,0)(),0(
2
1
?????
?????
tttu
tttu
10,44)()0,( 2 ?????? xxxxxu
20.0010,2
2
????
?
??
?
? tx
x
u
t
u 且
抛物型方程的古典显格式
解 r = k / h2 = 0.02 / 0.22 = 0.5,古典显式格式为
? ?nmnmnm UUU 111
2
1
??
? ??
…
┇ ┇
┇ ┇
…
0,0 0.2,0 0.4,0 1,0
0,0.2
┇
0,0.02
1,0.2
┇
1,0.02
抛物型方程的古典显格式
x=0.0 x=0.2 x=0.4 x=0.6 x=0.8 x=1.0
U =
t=0.00 0 0.6400 0.9600 0.9600 0.6400 0
t=0.02 0 0.4800 0.8000 0.8000 0.4800 0
t=0.04 0 0.4000 0.6400 0.6400 0.4000 0
t=0.06 0 0.3200 0.5200 0.5200 0.3200 0
t=0.08 0 0.2600 0.4200 0.4200 0.2600 0
t=0.10 0 0.2100 0.3400 0.3400 0.2100 0
t=0.12 0 0.1700 0.2750 0.2750 0.1700 0
t=0.14 0 0.1375 0.2225 0.2225 0.1375 0
t=0.16 0 0.1113 0.1800 0.1800 0.1113 0
t=0.18 0 0.0900 0.1456 0.1456 0.0900 0
t=0.20 0 0.0728 0.1178 0.1178 0.0728 0
抛物型方程的古典显格式
抛物型方程的古典显格式
function u=gu_dian(f,a,b,c1,c2,m,n)
%输入初值和 U
h=a/(m-1);
k=b/(n-1);
r=k/h^2;
U=zeros(n,m);
%赋边界条件
U(2:n,1)=c1;
U(2:n,m)=c2;
function y=fg(x)
y=4.*(x-x.^2);
抛物型方程的古典显格式
%赋初始条件
U(1,1:m)=fg(0:h:h*(m-1));
%计算内点上 u的数值解 U
for i=2:n
for j=2:m-1
U(i,j)=(1-2*r)*U(i-1,j)+r*(U(i-1,j-1)+U(i-1,j+1));
end
end
抛物型方程的古典显格式
% gu_dianl1.m 步长 h=0.20,k=0.02,r = k / h2 = 0.5
a=1;
b=0.20;
c1=0;
c2=0;
m=6;
n=11;
U=gu_dian('fg',a,b,c1,c2,m,n)
x=0:0.2:a;
y=0:0.02:b;
[X,Y]=meshgrid(x,y);
surf(X,Y,U) % 输入 U后再画图
抛物型方程的古典显格式
若取步长 h = 0.2,k = 0.0333,则 r = k / h2
=0.8333,
古典显式格式为
求解区域 Ω={0≤x≤1,0≤t≤0.3333},由
于 r =0.8333 <1/2，计算结果不稳定。
? ?nmnmnmnm UUUU 111 8 3 3 3.06 6 6 7.0 ??? ????
抛物型方程的古典显格式
x=0.0 x=0.20 x=0.40 x=0.60 x=0.80 x=1.00
U =
t=0.0000 0 0.6400 0.9600 0.9600 0.6400 0
t=0.0333 0 0.3734 0.6934 0.6934 0.3734 0
t=0.0667 0 0.3289 0.4267 0.4267 0.3289 0
t=0.1000 0 0.1364 0.3452 0.3452 0.1364 0
t=0.1333 0 0.1968 0.1712 0.1712 0.1968 0
t=0.1667 0 0.0115 0.1925 0.1925 0.0115 0
t=0.2000 0 0.1527 0.0417 0.0417 0.1527 0
t=0.2333 0 -0.0671 0.1342 0.1342 -0.0671 0
t=0.2667 0 0.1565 -0.0335 -0.0335 0.1565 0
t=0.3000 0 -0.1323 0.1249 0.1249 -0.1323 0
t=0.3333 0 0.1922 -0.0894 -0.0894 0.1922 0
抛物型方程的古典显格式
抛物型方程的古典显格式
% gu_dianl2.m 步长 h=0.20,k=0.0333,r = k / h2 = 0.8333
a=1;
b=0.3333;
c1=0;
c2=0;
m=6;
n=11;
U=gu_dian('fg',a,b,c1,c2,m,n)
x=0:0.2:a;
y=0:0.03333:b;
[X,Y]=meshgrid(x,y);
surf(X,Y,U)
抛物型方程的古典显格式
二、一维变系数热传导方程的显式格式
1,(2.36)
差分近似
方程 (2.36)的显式格式 (2.37)
截断误差为
0)()( 2
2
?
?
??
?
? xa
x
uxa
t
u
2
11
1 2
)(
h
uuuxa
k
uu nmnmnm
m
n
m
n
m ??
? ??
??
))(())(21( 111 nmnmmnmmnm UUxraUxraU ??? ????
).( 2hkO ?
抛物型方程的古典显格式
2,
差分近似
)38.2(0)()( ???
?
?
???
?
?
?
?
??
?
? xa
x
uxa
xt
u
2
2
)()()(
x
uxa
x
uxa
x
uxa
x ?
??
?
???
???
?
???
?
?
?
?
?
2
1111
1 2
2 h
uuua
h
uua
k
uu nmnmnmnmnmnmnm ????? ???????
抛物型方程的古典显格式
整理得方程 (2.38)的显式格式 (2.39)
截断误差为
作业 (P.108) 习题二 1
n
m
n
m
n
m
n
m
UaharUahar
UraU
11
1
)
2
1
()
2
1
(
)21(
??
?
??????
??
).( 2hkO ?
微分方程数值解
计算科学系 杨韧
抛物型方程的古典显格式
第二章 抛物型方程的差分格式
模型 长度为 L 的绝缘杆上的一维热流
x = 0
u(0,t)=c1 x =L x u(L,t)=c
2
杆
绝缘体
抛物型方程的古典显格式
热传导方程在时间 t 和位置 x 处的温度 u (x,t)
表示为
初始温度分布为
杆端点的边界值为
其中 k 是导热率系数,σ 是热量,ρ 是杆的密度。
?????
?
??
?
? tLx
t
txu
x
txuk 0,0),(),(
2
2
??
.0),()0,( Lxxxu ????
????
????
tctLu
tctu
0,),(
0,),0(
2
1
抛物型方程的古典显格式
§ 2.1 差分格式建立的基础
一、网格
将求解区域 Ω 分割成 M?N个小矩形，长宽分
别为
⊿ x=h
⊿ t=k
h k
0 x0 x1 … xM x
t
tN
┆
t0
抛物型方程的古典显格式
二、差商
一阶偏导数 的向前、向后、中心差商为
二阶偏导数 的中心差商为
n
mx
u ?
?
??
?
?
?
?
h
uu
h
uu
h
uu nmnmnmnmnmnm
2
,,1111 ???? ???
n
mx
u
???
?
???
?
?
?
2
2
2
11 2
h
uuu nmnmnm ?? ??
抛物型方程的古典显格式
用 (xm,tn)处的一阶向前差商近似代替微商 ux,
即
考虑 u(x,t) 的 Taylor 展式
)5.2(1
h
uu
x
u nmnmn
m
??
?
?
??
?
?
?
? ?
????
?
?
???
?
?
??
???
?
???
?
?
??
???
?
???
?
?
???
?
n
m
n
m
n
m
nmnm x
uh
x
uh
x
uhtxutxu
3
33
2
22
1 !3!2!1),(),(
????
?
?
???
?
?
??
???
?
???
?
?
??
???
?
???
?
?
???
?
n
m
n
m
n
m
nmnm x
uh
x
uh
x
uhtxutxu
3
33
2
22
1 !3!2!1),(),(
抛物型方程的古典显格式
截断误差
)(
!2
!3!2
!3!2!1
12
2
3
32
2
2
3
33
2
22
1
??
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
mm
t
x
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
n
m
xxx
x
uh
x
uh
x
uh
hu
x
uh
x
uh
x
uh
u
x
u
h
uu
x
u
E
n
?
?
抛物型方程的古典显格式
用 (xm,tn)处的一阶中心差商近似代替微商 ux,即
截断误差为
)7.2(
2
11
h
uu
x
u nmnmn
m
?? ???
?
?
?
?
?
?
?
n
t
x
n
m
n
m
n
m
n
m
n
m
n
m
x
uh
x
uh
x
uh
h
uu
x
u
E
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
3
32
5
54
3
32
11
!3!5!3
2
?
抛物型方程的古典显格式
三、算子
为 x 方向偏导数算子
Tx为 x 方向位移算子
μ x 为 x 方向平均算子
x
D x
?
??
n
m
n
mx
n
m
n
mx uuTuuT 1
1
1,?
?
? ??
? ?nmnmnmx uuu
2
1
2
12
1
?? ???
抛物型方程的古典显格式
x 方向差分算子
前差算子
后差算子
中心差算子
（ I为恒等算子）
nmnmnmxx uuu ?? ? 1,??
nmnmnmxx uuu 1,?????
n
m
n
m
n
mxx uuu 2121,?? ????
n
mx
n
mxx
n
m
n
m
n
m
n
m
n
m
uhDuD
h
D
h
I
x
uh
x
uh
x
uh
uu
)e x p (
!2!1
!3!2!1
2
2
3
33
2
22
1
??
?
?
?
?
?
????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
抛物型方程的古典显格式
§ 2.2 显式差分格式
差分方法：用差商代替微商，在网格节点上
求出微分方程解的近似值的一种方法。
一、一维常系数热传导方程的古典显式格式
考虑一维热传导方程 (2.26)
在网格节点
2
2
x
u
t
u
?
??
?
?
m– 1,n m,n m+1,n
m,n+1
抛物型方程的古典显格式
的差分近似
取 r = k / h2 为步长比，得显式向前差分方程
(2.29)
(2.29)也称为解热传导方程 (2.26)的古典显式格
式。截断误差为
2
11
1 2
h
uuu
k
uu nmnmnmnmnm ??? ????
n
m
n
m
n
m
n
m rUUrrUU 11
1 )21(
??
? ????
).( 2hkO ?
抛物型方程的古典显格式
例 1 用古典显式格式求解抛物型方程
初始条件为
边界条件为
取步长 ⊿x = h = 0.2,⊿t = k = 0.02 。
20.00,0)(),1(
20.00,0)(),0(
2
1
?????
?????
tttu
tttu
10,44)()0,( 2 ?????? xxxxxu
20.0010,2
2
????
?
??
?
? tx
x
u
t
u 且
抛物型方程的古典显格式
解 r = k / h2 = 0.02 / 0.22 = 0.5,古典显式格式为
? ?nmnmnm UUU 111
2
1
??
? ??
…
┇ ┇
┇ ┇
…
0,0 0.2,0 0.4,0 1,0
0,0.2
┇
0,0.02
1,0.2
┇
1,0.02
抛物型方程的古典显格式
x=0.0 x=0.2 x=0.4 x=0.6 x=0.8 x=1.0
U =
t=0.00 0 0.6400 0.9600 0.9600 0.6400 0
t=0.02 0 0.4800 0.8000 0.8000 0.4800 0
t=0.04 0 0.4000 0.6400 0.6400 0.4000 0
t=0.06 0 0.3200 0.5200 0.5200 0.3200 0
t=0.08 0 0.2600 0.4200 0.4200 0.2600 0
t=0.10 0 0.2100 0.3400 0.3400 0.2100 0
t=0.12 0 0.1700 0.2750 0.2750 0.1700 0
t=0.14 0 0.1375 0.2225 0.2225 0.1375 0
t=0.16 0 0.1113 0.1800 0.1800 0.1113 0
t=0.18 0 0.0900 0.1456 0.1456 0.0900 0
t=0.20 0 0.0728 0.1178 0.1178 0.0728 0
抛物型方程的古典显格式
抛物型方程的古典显格式
function u=gu_dian(f,a,b,c1,c2,m,n)
%输入初值和 U
h=a/(m-1);
k=b/(n-1);
r=k/h^2;
U=zeros(n,m);
%赋边界条件
U(2:n,1)=c1;
U(2:n,m)=c2;
function y=fg(x)
y=4.*(x-x.^2);
抛物型方程的古典显格式
%赋初始条件
U(1,1:m)=fg(0:h:h*(m-1));
%计算内点上 u的数值解 U
for i=2:n
for j=2:m-1
U(i,j)=(1-2*r)*U(i-1,j)+r*(U(i-1,j-1)+U(i-1,j+1));
end
end
抛物型方程的古典显格式
% gu_dianl1.m 步长 h=0.20,k=0.02,r = k / h2 = 0.5
a=1;
b=0.20;
c1=0;
c2=0;
m=6;
n=11;
U=gu_dian('fg',a,b,c1,c2,m,n)
x=0:0.2:a;
y=0:0.02:b;
[X,Y]=meshgrid(x,y);
surf(X,Y,U) % 输入 U后再画图
抛物型方程的古典显格式
若取步长 h = 0.2,k = 0.0333,则 r = k / h2
=0.8333,
古典显式格式为
求解区域 Ω={0≤x≤1,0≤t≤0.3333},由
于 r =0.8333 <1/2，计算结果不稳定。
? ?nmnmnmnm UUUU 111 8 3 3 3.06 6 6 7.0 ??? ????
抛物型方程的古典显格式
x=0.0 x=0.20 x=0.40 x=0.60 x=0.80 x=1.00
U =
t=0.0000 0 0.6400 0.9600 0.9600 0.6400 0
t=0.0333 0 0.3734 0.6934 0.6934 0.3734 0
t=0.0667 0 0.3289 0.4267 0.4267 0.3289 0
t=0.1000 0 0.1364 0.3452 0.3452 0.1364 0
t=0.1333 0 0.1968 0.1712 0.1712 0.1968 0
t=0.1667 0 0.0115 0.1925 0.1925 0.0115 0
t=0.2000 0 0.1527 0.0417 0.0417 0.1527 0
t=0.2333 0 -0.0671 0.1342 0.1342 -0.0671 0
t=0.2667 0 0.1565 -0.0335 -0.0335 0.1565 0
t=0.3000 0 -0.1323 0.1249 0.1249 -0.1323 0
t=0.3333 0 0.1922 -0.0894 -0.0894 0.1922 0
抛物型方程的古典显格式
抛物型方程的古典显格式
% gu_dianl2.m 步长 h=0.20,k=0.0333,r = k / h2 = 0.8333
a=1;
b=0.3333;
c1=0;
c2=0;
m=6;
n=11;
U=gu_dian('fg',a,b,c1,c2,m,n)
x=0:0.2:a;
y=0:0.03333:b;
[X,Y]=meshgrid(x,y);
surf(X,Y,U)
抛物型方程的古典显格式
二、一维变系数热传导方程的显式格式
1,(2.36)
差分近似
方程 (2.36)的显式格式 (2.37)
截断误差为
0)()( 2
2
?
?
??
?
? xa
x
uxa
t
u
2
11
1 2
)(
h
uuuxa
k
uu nmnmnm
m
n
m
n
m ??
? ??
??
))(())(21( 111 nmnmmnmmnm UUxraUxraU ??? ????
).( 2hkO ?
抛物型方程的古典显格式
2,
差分近似
)38.2(0)()( ???
?
?
???
?
?
?
?
??
?
? xa
x
uxa
xt
u
2
2
)()()(
x
uxa
x
uxa
x
uxa
x ?
??
?
???
???
?
???
?
?
?
?
?
2
1111
1 2
2 h
uuua
h
uua
k
uu nmnmnmnmnmnmnm ????? ???????
抛物型方程的古典显格式
整理得方程 (2.38)的显式格式 (2.39)
截断误差为
作业 (P.108) 习题二 1
n
m
n
m
n
m
n
m
UaharUahar
UraU
11
1
)
2
1
()
2
1
(
)21(
??
?
??????
??
).( 2hkO ?
