抛物型方程的差分格式

u t

2u x2
(2.26)
在网格节点
m , n+1
m–1, n m , n m+1 , n
抛物型方程的古典显格式
的差分近似
umn1 umn umn 1 2umn umn 1
k
h2
取 r = k / h2 为步长比，得显式向前差分方程
U n1 m

rU
n m1
n m
的中心差商为
umn 1 2umn umn 1
h2
抛物型方程的古典显格式
用(xm , tn)处的一阶向前差商近似代替微商ux，
即

u
n


umn 1 umn
(8.5)
x m
h
考虑 u(x , t) 的 Taylor 展式
u(
xm1,
tn
)

u ( xm
,
tn
U
n1 m

1 2
U
n m1

U
n m1
0 , 0.2
┇
┇
0,0.02
பைடு நூலகம்
┇
0 , 0 0.2 , 0
…
┇ ┇
…
0.4 , 0
1 , 0.2 ┇
1, 0.02
1,0
抛物型方程的古典显格式
x=0.0 x=0.2 U= t=0.00 0 0.6400 t=0.02 0 0.4800 t=0.04 0 0.4000 t=0.06 0 0.3200 t=0.08 0 0.2600 t=0.10 0 0.2100 t=0.12 0 0.1700 t=0.14 0 0.1375 t=0.16 0 0.1113 t=0.18 0 0.0900 t=0.20 0 0.0728
抛物型方程的古典显格式
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
n
m


I

h 1!
Dx

h2 2!
Dx2
umn

exp( hDx )umn
（I为恒等算子）
抛物型方程的古典显格式
8.2 显式差分格式
差分方法：用差商代替微商，在网格节点上 求出微分方程解的近似值的一种方法。
一、一维常系数热传导方程的古典显式格式
考虑一维热传导方程
x=0.8
0.6400 0.4800 0.4000 0.3200 0.2600 0.2100 0.1700 0.1375 0.1113 0.0900 0.0728
x=1.0
0 0 0 0 0 0 0 0 0 0 0
抛物型方程的古典显格式
抛物型方程的古典显格式
function u=gu_dian(f,a,b,c1,c2,m,n)

(1
2r
)U
n m

rU
n m1
(2.29)
(2.29)也称为解热传导方程(2.26)的古典显式格
式。截断误差为 O(k h2 ).
抛物型方程的古典显格式
例1 用古典显式格式求解抛物型方程
u 2u t x2 ,
0 x 1且0 t 0.20
初始条件为 u(x,0) (x) 4x 4x2, 0 x 1
u n

umn 1 umn 1
x m
2h
(8.7)
截断误差为
Emn


n
u x

m

umn 1 umn 1 2h

h2 3!

3u x3
n m

h4 5!

5u x5
n

m

h2 3!

3u x3
前差算子
x,
xumn

un m1
umn
后差算子
x,
xumn

umn

un m1
中心差算子
x , xumn
un
m

1 2
un
m

1 2
un m1

umn

h 1!

u x
n m

h2 2!

2u x2
n
m

h3 3!

3u x3
n m



umn

h


h 2!


2u x2
n
m

h2 3!


3u x3
n
m




h 2!


2u x2
tn x
(xm x xm1)
抛物型方程的古典显格式
用(xm , tn)处的一阶中心差商近似代替微商ux，即
x=0.4
0.9600 0.8000 0.6400 0.5200 0.4200 0.3400 0.2750 0.2225 0.1800 0.1456 0.1178
x=0.6
0.9600 0.8000 0.6400 0.5200 0.4200 0.3400 0.2750 0.2225 0.1800 0.1456 0.1178
%输入初值和U
h=a/(m-1);
k=b/(n-1);
r=k/h^2;
U=zeros(n,m); %赋边界条件 U(2:n,1)=c1;
function y=fg(x) y=4.*(x-x.^2);
U(2:n,m)=c2;
抛物型方程的古典显格式
%赋初始条件
U(1,1:m)=fg(0:h:h*(m-1));
a umn 1 umn 1 2h

a
umn 1

2umn h2
umn 1
抛物型方程的古典显格式
整理得方程(2.38)的显式格式(2.39)
U n1 m

(1
2ra)U
n m

r
(a

1 2
ha)U
n m1

r
(a

1 2
ha)U
n m1
截断误差为 O(k h2 ).
抛物型方程的古典显格式
若取步长 h = 0.2 , k = 0.0333 , 则 r = k / h2 =0.8333，
古典显式格式为
U n1 m

0.6667U
n m

0.8333
Un m1

U
n m1
求解区域 Ω ={0≤x≤1 ,0≤t≤0.3333},由
于r =0.8333 <1/2，计算结果不稳定。
边界条件为 u(0,t) 1(t) 0, 0 t 0.20 u(1,t) 2(t) 0, 0 t 0.20
取步长⊿x = h = 0.2 , ⊿t = k = 0.02 。
抛物型方程的古典显格式
解 r = k / h2 = 0.02 / 0.22 = 0.5, 古典显式格式为
抛物型方程的古典显格式
抛物型方程的古典显格式
% 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)
第8章 抛物型方程的差分格式
模型 长度为 L 的绝缘杆上的一维热流
绝缘体
x=0
u(0 , t)=c1
杆
x =L u(L , t)=c2
x
抛物型方程的古典显格式
热传导方程在时间 t 和位置 x 处的温度 u (x , t)
表示为
k 2u(x,t) u(x,t) 0 x L,0 t

ra(
xm
)(U
n m1

U
n m1
)
截断误差为 O(k h2 ).
抛物型方程的古典显格式
2、u t

x
a(x)
u x

a(x) 0
(2.38)
x

a(x)
u x


a(x)
u x

a(x)

2u x2
差分近似
umn1 umn k
3u x3
n m

抛物型方程的古典显格式
截断误差
Emn


u x
n
m

un m1

umn
h


u x
n
m

umn

h 1!

u x
nm

h2 2!

2u x2
n m

h3 3!

3u x3
抛物型方程的古典显格式
二、一维变系数热传导方程的显式格式
1、
u t

a(x)
2u x2
a(x) 0
差分近似
(2.36)
umn1 umn k

a(
xm
)
umn 1

2umn h2
umn 1
方程(2.36)的显式格式(2.37)
U n1 m

(1
2ra(
xm
))U
n m
tn x
抛物型方程的古典显格式
三、算子
Dx

x
为 x 方向偏导数算子
Tx为 x 方向位移算子
Txumn umn 1, Tx1umn umn 1
μ x 为 x 方向平均算子

xu
n m

1 2
un
m

1 2
un
m

1 2
抛物型方程的古典显格式
x 方向差分算子
%计算内点上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后再画图
别为
⊿x=h ⊿t=k
t tN
┆
t0
k
h
0 x0 x1 …
xM x
抛物型方程的古典显格式
二、差商
一阶偏导数
u
n

的向前、向后、中心差商为
x m
umn 1 umn , umn umn 1 , umn 1 umn 1
h
h
2h
二阶偏导数

2u x 2

x2
t
初始温度分布为 u(x,0) (x), 0 x L.
杆端点的边界值为
u(0,t) c1, 0 t
u(L,t) c2, 0 t
其中 k 是导热率系数，σ 是热量，ρ 是杆的密度。
抛物型方程的古典显格式
8.1 差分格式建立的基础
一、网格
将求解区域Ω 分割成 M∙N个小矩形，长宽分
抛物型方程的古典显格式
)

h 1!

u x
n
m

h2 2!

2u x2

n m

h3 3!

3u x3

n m


u(xm1, tn
)

u ( xm
,tn
)

h 1!

u x
n
m

h2 2!

2u x2
n m

h3 3!