偏微分方程数值解法(抛物型方程差分法)2

系数矩阵可逆
k1 [ A(k ) ]1 B(k ) k
记 H (k ) [ A(k ) ]1 B(k ) 称之为过渡矩阵
k1 H (k ) k 常系数差分格式 k1 H k
H 的谱半径:

(
H
)

max
1 jn
|

j
(
H
)
|
定理: 双层差分格式稳定的必要条件是,存在与 无
1 2ra 2 sin2 ( j / 2(n 1)) 1 2ra 2 sin2 ( j / 2(n 1))
13/17
过渡矩阵的谱半径

(
H
)

max
1 jn
|

j
(
H
)
|
1 2ra 2 sin2 ( j / 2(n 1))

max
1 jn
|
1

2ra 2
sin2 (
x=0:h:1;N=length(x);
t=0;uk=sin(pi*x);
II=2:N-1;k=0;
while t<T
t=t+ta;
un=s*uk(II)+r*(uk(II-1)+uk(II+1));
uk=[0,un,0];k=k+1;
end
ux=exp(-pi*pi*t).*sin(pi*x);
error=max(abs(ux-uk))
显格式 隐格式
1

[ukj
1

u
k j
]

a2 h2

2 x
u
k j
1

[ukj
1

u
k j
]

a2 h2
u2 k1 xj
C-N格式
1

[ukj
1

u
k j
]

a2 2h2

2 x
(
ukj
1

ukj )
15/17
数值计算实验 显格式: input T:=1 error = 7.9443e-006 k = 200
nn
特征值

j
(1
1
2ra 2 )
2ra 2[1
2ra 2 cos
cos( j )]
( j
n
1
)
n1
1 4ra 2 sin2 j
2(n 1)
9/17
过渡矩阵的谱半径
(H)
max

|1
max
1 jn
4ra
|
2
j(H
sin2
)
| j
|
1 jn
n
n
u (k ) k1 jm m

u (k ) k
jm m
f
k j
m 1
m1
记矩阵
A( k )

(
(k) jm
)nn
B( k )

( ) (k) jm nn
双层格式的矩阵形式 A(k )uk1 B(k )uk f k
双层差分格式初值稳定概念:
A(k ) k1 B(k ) k 任意解都满足 || k || M || k0 ||
关的常数 M1 ,使得
(H ) 1 M1
7/17
定理 若 H = A-1B 为正规矩阵,即 HH* = H*H,
则条件 (H ) 1 M1
是双层差分格式按欧氏范数稳定的充分条件 注：欧氏范数(或离散L2范数)
n
|| uk ||0 [h (ukj )2 ]1/ 2 j1

(1
2ra2 )ukj

ra 2 (ukj 1

uk j 1
)


f
k j
在实际应用时,取逐层计算形式.当初始层数据有误 差时,误差会逐层传播,影响以后各层的解.
记
ukj
的误差为

k j
,
设
f
k j
无误差,
则有

k 1 j

(1
2ra2 )
k j

ra 2 (
k j 1


) k
n1
显式差分格式稳定充分条件. h2 / 2a2
10/17
简单隐式差分格式矩阵形式
[(1 2ra 2 )I ra 2C ]uhk1 uhk
f
k h
1
过渡矩阵 H [(1 2ra 2 )I ra 2C ]1
特征值
1 2ra ra
1

其中 M 与 无关. k > k0
5/17
简单显式差分格式
uk1 j

(1
2ra2 )ukj

ra 2 (ukj 1

uk j 1
)
u0j ( x j )
u0k

uk n1

0
( j 1,2, , n)
稳定性分析,设 ra 2 1 / 2
|
uk1 j
|
(1
2ra2 ) |
C—N格式 input T:=1 error = 2.6227e-006 k = 50
x 10-5 6
4
2
0
0
0.2
0.4
0.6
0.8
1
x 10-5 6
4
2
0
0
0.2
0.4
0.6
0.8
1
16/17
T=input('input T:='); h=1/10;ta=1/200;
%显格式计算程序
r=ta/(h*h);s=1-2*r;
plot(x,ux,x,uk,':or')
17/17
《偏微分方程数值解法》 7
——抛物型方程差分法2
差分格式稳定性概念 显、隐格式稳定性分析 稳定性分析的矩阵方法
1/17
抛物型方程
u a2 t
2u x 2

f (x,t)
t
简单显式差分格式
x
uk 1 j

ukj


a2 h2

u2 k
xj

f
k j
r / h2

uk1 j
(
H
)

max
1 j N 1
|

j
(H
)
|

max
1 jn
|
1
4ra
1
sin2 ( j
/
2(n

1))
|
1
1 4ra sin2( / 2(n 1)) 1
隐式差分格式无条件稳定.
12/17
C-N 格式矩阵形式
[(1 ra 2 )I ra 2 C ]uk1
2 [(1 ra 2 )I
ra
1 2ra ra




ra
1

2ra


j
(H)
[1


[(1 2ra 2 )
2ra(1 cos
2ra 2 cos j )]1
j
n
]1 1
[1 4ra 2 sin2 nj 1 ]1 1
2(n 1)
11/17
过渡矩阵的谱半径

(1
2ra2 )
k j

ra 2 (
k j 1

) k
j 1


k 1 j
源自文库


k j

(
k j 1

) k
j 1
设初始层上,仅有
0 j0

,其它点处无误差
在各计算层上,误差传播没有得到控制
4/17
无穷大范数定义 ||
uk
||
max |
1 jn
ukj
|
双层差分格式
j 1
2/17
取 ra 2 1 / 2

k 1 j

(1
2ra2 )
k j

ra 2 (
k j 1

) k
j 1


k 1 j

1 2
(
k j1

k j 1
)
设初始层上,仅有
0 j0

,其它点处无误差
在各计算层上,误差传播得到控制
3/17
取
ra 2 1

k 1 j

ra 2
C ]uk

(fk

f k1 )
2
2
H [(1 ra 2 )I ra 2 C ]1[(1 ra 2 )I ra 2 C ]
2
2
特征值
j

(1 (1
ra 2 ) ra 2 )
ra 2 ra 2
cos( cos(
j j
/(n 1)) /(n 1))
8/17
简单显式差分格式矩阵形式
uk1 [(1 2ra2 )I ra2C]uk f k
过渡矩阵 H [(1 2ra 2 )I ra 2C ]
1 2ra 2 ra 2




ra 2
1 2ra 2 ra 2


ra 2
1

2ra
2

2(n 1)
| (ra 2 ) |
(ra 2 )
| n |
| 1 |
ra 2
极值点满足
ra 2 1 2
1 4ra 2 sin2 4ra 2 sin2 n 1
2(n 1)
2(n 1)
(ra 2 ) 1 2sin2 cos 1
2(n 1)
ukj
|
ra 2 (|
uk j 1
|

|
uk j 1
|)
(1 2ra2 ) || uk || 2ra2 || uk ||

|
uk1 j
|||
uk
||
|| uk1 ||C || uk || 此时差分格式稳定
6/17
设齐次方程 A(k ) k1 B(k ) k
j
/
2(n

1))
|
(H) 1
C-N 格式是无条件稳定的.
14/17
数值实验题 用三种差分格式求 初边值问题数值解
ut uxx , 0 x 1, t 0
u( x,0) sinx, 0 x 1
u(0, t) u(1, t) 0, t 0
并与准确解比较 u( x, t) exp( 2t)sin x