线性代数知识
向量的范数
设x和y为n维实向量,则x和y的任意范数具有以下基本性质:
正定性:
∣
∣
x
∣
∣
≥
0
||x||≥0
∣∣x∣∣≥0,只有x为零向量时||x||=0
齐次性:∀k∈R,
∣
∣
k
x
∣
∣
=
∣
k
∣
⋅
∣
∣
x
∣
∣
||kx||=|k|·||x||
∣∣kx∣∣=∣k∣⋅∣∣x∣∣
三角不等式:
∣
∣
x
+
y
∣
∣
≤
∣
∣
x
∣
∣
+
∣
∣
y
∣
∣
||x+y||≤||x||+||y||
∣∣x+y∣∣≤∣∣x∣∣+∣∣y∣∣
且只要与向量x对应的实值函数||x||满足上述四条性质,它就被定义为x的一种向量范数。
常用范数的计算:
1-范数:
∣
∣
x
∣
∣
1
=
∑
∣
x
i
∣
1
=
∑
∣
x
i
∣
||x||_1=\sqrt[1]{∑|x_i|}=∑|x_i|
∣∣x∣∣1=1∑∣xi∣=∑∣xi∣,史称“街区距离”、“曼哈顿距离”
2-范数:
∣
∣
x
∣
∣
2
=
∑
∣
x
i
∣
2
=
∑
x
i
2
||x||_2=\sqrt{∑|x_i|^2}=\sqrt{∑x_i^2}
∣∣x∣∣2=∑∣xi∣2=∑xi2,史称“向量的长度”、“向量的模”、“欧氏距离”,最常用
p-范数:
∣
∣
x
∣
∣
p
=
∑
∣
x
i
∣
p
p
||x||_p=\sqrt[p]{∑|x_i|^p}
∣∣x∣∣p=p∑∣xi∣p,囊括了其它三种范数
∞-范数:
∣
∣
x
∣
∣
∞
=
lim
p
→
+
∞
∑
∣
x
i
∣
p
p
=
max
(
x
1
,
x
2
,
⋯
,
x
n
)
||x||_∞=\lim\limits_{p→+∞}\sqrt[p]{∑|x_i|^p}=\max(x_1,x_2,⋯,x_n)
∣∣x∣∣∞=p→+∞limp∑∣xi∣p=max(x1,x2,⋯,xn),史称“棋盘距离”、“切比雪夫距离”
范数的性质:
连续性: 对于实向量x,其范数||x||是对于x的各个分量的n元连续函数
等价性: 对于Rn上的任意两种范数||x||a和||x||b,∃m,n∈R,使∀x∈Rn,都有
∣
∣
x
∣
∣
b
∈
[
m
∣
∣
x
∣
∣
a
,
n
∣
∣
x
∣
∣
a
]
||x||_b∈[m||x||_a,n||x||_a]
∣∣x∣∣b∈[m∣∣x∣∣a,n∣∣x∣∣a]
按范数收敛性: 向量序列{x1,x2,⋯,xn}收敛于精确值x*⇔对于某种范数有
lim
k
→
∞
∣
∣
x
∗
−
x
k
∣
∣
=
0
\lim\limits_{k→∞}||x^*-x_k||=0
k→∞lim∣∣x∗−xk∣∣=0⇔对于任意范数有
lim
k
→
∞
∣
∣
x
∗
−
x
k
∣
∣
=
0
\lim\limits_{k→∞}||x^*-x_k||=0
k→∞lim∣∣x∗−xk∣∣=0
矩阵的范数
设A和B为n阶实方阵,则A和B的任意范数具有以下基本性质:
正定性:
∣
∣
A
∣
∣
≥
0
||A||≥0
∣∣A∣∣≥0,只有A=O时||A||=0
齐次性:∀k∈R,
∣
∣
k
A
∣
∣
=
∣
k
∣
⋅
∣
∣
A
∣
∣
||kA||=|k|·||A||
∣∣kA∣∣=∣k∣⋅∣∣A∣∣
三角不等式:
∣
∣
A
+
B
∣
∣
≤
∣
∣
A
∣
∣
+
∣
∣
B
∣
∣
||A+B||≤||A||+||B||
∣∣A+B∣∣≤∣∣A∣∣+∣∣B∣∣
相容性:
∣
∣
A
B
∣
∣
≤
∣
∣
A
∣
∣
⋅
∣
∣
B
∣
∣
||AB||≤||A||·||B||
∣∣AB∣∣≤∣∣A∣∣⋅∣∣B∣∣
且只要与矩阵A对应的实值函数||A||满足上述四条性质,它就被定义为A的一种矩阵范数。
若实方阵A和实向量x满足
∣
∣
A
x
∣
∣
≤
∣
∣
A
∣
∣
⋅
∣
∣
x
∣
∣
||Ax||≤||A||·||x||
∣∣Ax∣∣≤∣∣A∣∣⋅∣∣x∣∣,则||A||是与向量范数||x||相容的矩阵范数。
对于算子范数有||E||=1,对于任意范数有||E||≥ρ(E)=1。
常用范数的计算:
行范数:
∣
∣
A
∣
∣
∞
=
max
(
∑
j
=
1
n
∣
a
1
j
∣
,
∑
j
=
1
n
∣
a
2
j
∣
,
⋯
,
∑
j
=
1
n
∣
a
n
j
∣
)
||A||_∞=\max(∑_{j=1}^n|a_{1j}|,∑_{j=1}^n|a_{2j}|,⋯,∑_{j=1}^n|a_{nj}|)
∣∣A∣∣∞=max(∑j=1n∣a1j∣,∑j=1n∣a2j∣,⋯,∑j=1n∣anj∣)
列范数:
∣
∣
A
∣
∣
1
=
max
(
∑
i
=
1
n
∣
a
i
1
∣
,
∑
i
=
1
n
∣
a
i
2
∣
,
⋯
,
∑
i
=
1
n
∣
a
i
n
∣
)
||A||_1=\max(∑_{i=1}^n|a_{i1}|,∑_{i=1}^n|a_{i2}|,⋯,∑_{i=1}^n|a_{in}|)
∣∣A∣∣1=max(∑i=1n∣ai1∣,∑i=1n∣ai2∣,⋯,∑i=1n∣ain∣)
2-范数:
∣
∣
A
∣
∣
2
=
ρ
(
A
T
A
)
||A||_2=\sqrt{ρ(A^TA)}
∣∣A∣∣2=ρ(ATA)。A为对称矩阵时,
∣
∣
A
∣
∣
2
=
ρ
(
A
)
||A||_2=ρ(A)
∣∣A∣∣2=ρ(A)
F-范数:
∣
∣
A
∣
∣
F
=
∑
a
i
j
2
||A||_F=\sqrt{∑ a_{ij}^2}
∣∣A∣∣F=∑aij2
矩阵的谱半径
谱半径的定义为 ρ ( A ) = max ( ∣ λ 1 ∣ , ∣ λ 2 ∣ , ⋯ , ∣ λ n ∣ ) ρ(A)=\max(|λ_1|, |λ_2|,⋯,|λ_n|) ρ(A)=max(∣λ1∣,∣λ2∣,⋯,∣λn∣)。
设A为n阶实方阵,则A的谱半径具有以下基本性质:
(1)对任意一种算子范数||A||,均有
ρ
(
A
)
≤
∣
∣
A
∣
∣
ρ(A)≤||A||
ρ(A)≤∣∣A∣∣
(2)∀ε>0,必存在一种或多种算子范数|A||,使
ρ
(
A
)
+
ε
≥
∣
∣
A
∣
∣
ρ(A)+ε≥||A||
ρ(A)+ε≥∣∣A∣∣
直接法解方程组
顺序Gauss消元法
最普通、最没有技术含量的解法。
先将A写成增广矩阵[A,b],然后从左到右一列列地消元成上三角阵,再从右往左一列列地回代成对角阵。
也可以从左到右一列列地消元成下三角阵,再从右往左一列列地前推成对角阵。
列主元Gauss消元法
列主元高斯消元法中,每次消元时把该列绝对值最大的一行换到对角线上。
在顺序高斯消元法中,系数矩阵A的顺序主子式不为0时,才能保证消元过程中的对角元素不会为0,而使用列主元高斯消元法时只需|A|≠0即可。
三角分解法
Doolittle分解法:
设L为单位下三角阵,U为上三角阵,则可唯一地将A分解为LU。于是,Ax=b即变为LUx=b,解x需依次解Ly=b和Ux=y。
其分解形式为
A
=
[
1
×
1
×
×
1
]
[
×
×
×
×
×
×
]
A=\begin{bmatrix}1\\×&1\\×&×&1\end{bmatrix}\begin{bmatrix}×&×&×\\&×&×\\&&×\end{bmatrix}
A=
1××1×1
××××××
。
Crout分解法:
设L’为下三角阵,U为单位上三角阵,则可唯一地将A分解为L’U’。于是,Ax=b即变为L’U’x=b,解x需依次解L’y=b和U’x=y。
其分解形式为
A
=
[
×
×
×
×
×
×
]
[
1
×
×
1
×
1
]
A=\begin{bmatrix}×\\×&×\\×&×&×\end{bmatrix}\begin{bmatrix}1&×&×\\&1&×\\&&1\end{bmatrix}
A=
××××××
1×1××1
。
LDU分解法:
设L为单位下三角阵,D为对角矩阵,U’为单位上三角阵,则可唯一地将A分解为LDU’。于是,Ax=b即变为LDU’x=b,解x需依次解Lz=b、Dy=z和U’x=y。
其分解形式为
A
=
[
1
×
1
×
×
1
]
[
×
×
×
]
[
1
×
×
1
×
1
]
A=\begin{bmatrix}1\\×&1\\×&×&1\end{bmatrix}\begin{bmatrix}×\\&×\\&&×\end{bmatrix}\begin{bmatrix}1&×&×\\&1&×\\&&1\end{bmatrix}
A=
1××1×1
×××
1×1××1
。
追赶法与平方根法
三对角方程组的系数矩阵的形式:
A
=
[
b
1
c
1
a
2
b
2
c
2
a
3
b
3
c
3
a
4
b
4
]
A=\begin{bmatrix}b_1&c_1\\a_2&b_2&c_2\\&a_3&b_3&c_3\\&&a_4&b_4\end{bmatrix}
A=
b1a2c1b2a3c2b3a4c3b4
追赶法实际上就是用基于Dooliitle分解的顺序高斯消元法解三对角方程组。
平方根法
当方程组的系数矩阵A对称正定时,根据LDU分解可将其分解为A=LDU’=AT=UT’DLT=LDLT(矩阵的LDU分解结果唯一)=L√D·√DLT=(L√D)(L√D)T。
因此,对称正定矩阵A可用Cholesky分解法分解成两个对称的三角阵,其形式为
A
=
L
1
L
1
T
=
[
a
11
a
21
a
22
a
31
a
32
a
33
]
[
a
11
a
21
a
31
a
22
a
32
a
33
]
A=L_1L_1^T=\begin{bmatrix}a_{11}\\a_{21}&a_{22}\\a_{31}&a_{32}&a_{33}\end{bmatrix}\begin{bmatrix}a_{11}&a_{21}&a_{31}\\&a_{22}&a_{32}\\&&a_{33}\end{bmatrix}
A=L1L1T=
a11a21a31a22a32a33
a11a21a22a31a32a33
。
在作Cholesky分解时,对角元素开根号后取正值。反之,若对称阵能被分解成上述形式,则其正定。
直接法的误差与条件数
用直接法解方程Ax+b的误差估计式为
∣
∣
δ
x
∣
∣
∣
∣
x
∣
∣
≤
∣
∣
A
−
1
∣
∣
⋅
∣
∣
A
∣
∣
1
−
∣
∣
A
−
1
∣
∣
⋅
∣
∣
δ
A
∣
∣
(
∣
∣
δ
A
∣
∣
∣
∣
A
∣
∣
+
∣
∣
δ
b
∣
∣
∣
∣
b
∣
∣
)
\frac{||δx||}{||x||}≤\frac{||A^{-1}||·||A||}{1-||A^{-1}||·||δA||}(\frac{||δA||}{||A||}+\frac{||δb||}{||b||})
∣∣x∣∣∣∣δx∣∣≤1−∣∣A−1∣∣⋅∣∣δA∣∣∣∣A−1∣∣⋅∣∣A∣∣(∣∣A∣∣∣∣δA∣∣+∣∣b∣∣∣∣δb∣∣),其中
∣
∣
δ
x
∣
∣
∣
∣
x
∣
∣
,
∣
∣
δ
A
∣
∣
∣
∣
A
∣
∣
,
∣
∣
δ
b
∣
∣
∣
∣
b
∣
∣
\frac{||δx||}{||x||},\frac{||δA||}{||A||},\frac{||δb||}{||b||}
∣∣x∣∣∣∣δx∣∣,∣∣A∣∣∣∣δA∣∣,∣∣b∣∣∣∣δb∣∣为x、A、b的相对误差。
若A精确,b有扰动,则上式变为
∣
∣
δ
x
∣
∣
∣
∣
x
∣
∣
≤
∣
∣
A
−
1
∣
∣
⋅
∣
∣
A
∣
∣
⋅
∣
∣
δ
b
∣
∣
∣
∣
b
∣
∣
\frac{||δx||}{||x||}≤||A^{-1}||·||A||·\frac{||δb||}{||b||}
∣∣x∣∣∣∣δx∣∣≤∣∣A−1∣∣⋅∣∣A∣∣⋅∣∣b∣∣∣∣δb∣∣。
若A有扰动,b精确,则上式变为
∣
∣
δ
x
∣
∣
∣
∣
x
∣
∣
≤
∣
∣
A
−
1
∣
∣
⋅
∣
∣
A
∣
∣
1
−
∣
∣
A
−
1
∣
∣
⋅
∣
∣
δ
A
∣
∣
⋅
∣
∣
δ
A
∣
∣
∣
∣
A
∣
∣
\frac{||δx||}{||x||}≤\frac{||A^{-1}||·||A||}{1-||A^{-1}||·||δA||}·\frac{||δA||}{||A||}
∣∣x∣∣∣∣δx∣∣≤1−∣∣A−1∣∣⋅∣∣δA∣∣∣∣A−1∣∣⋅∣∣A∣∣⋅∣∣A∣∣∣∣δA∣∣。
若||A-1||·||δA||很小的同时b精确,则上式变为
∣
∣
δ
x
∣
∣
∣
∣
x
∣
∣
≈
∣
∣
A
−
1
∣
∣
⋅
∣
∣
A
∣
∣
⋅
∣
∣
δ
A
∣
∣
∣
∣
A
∣
∣
=
∣
∣
A
−
1
∣
∣
⋅
∣
∣
A
∣
∣
\frac{||δx||}{||x||}≈||A^{-1}||·||A||·\frac{||δA||}{||A||}=||A^{-1}||·||A||
∣∣x∣∣∣∣δx∣∣≈∣∣A−1∣∣⋅∣∣A∣∣⋅∣∣A∣∣∣∣δA∣∣=∣∣A−1∣∣⋅∣∣A∣∣。
条件数的定义为
c
o
n
d
(
A
)
=
∣
∣
A
−
1
∣
∣
⋅
∣
∣
A
∣
∣
\mathrm{cond}(A)=||A^{-1}||·||A||
cond(A)=∣∣A−1∣∣⋅∣∣A∣∣,相当于直接法解方程组的误差的“放大系数”。
若条件数很小,则A为良态矩阵,相应的方程组Ax=b则为良态方程组。
若条件数很大,则A为病态矩阵或病态条件,相应的方程组Ax=b则为病态方程组。
无论||A||为哪种算子范数,其对应的条件数均满足
c
o
n
d
(
A
)
≥
1
\mathrm{cond}(A)≥1
cond(A)≥1。
迭代法解方程组
迭代法简介
解线性代数方程组Ax=b时,可以通过构造形如
x
k
+
1
=
B
x
k
+
f
,
k
∈
N
x_{k+1}=Bx_k+f,k∈N
xk+1=Bxk+f,k∈N的基本迭代公式求方程组的近似解。
在基本型迭代公式(也称单步定常线性迭代公式)中,如果
lim
k
→
∞
x
k
=
x
∗
\lim\limits_{k→∞}x_k=x^*
k→∞limxk=x∗,则该迭代收敛于精确解。
Jacobi迭代法
首先,将A分裂为D、L、U三部分:
A
=
D
−
L
−
U
=
[
a
11
a
22
a
33
]
−
[
0
−
a
21
0
−
a
31
−
a
32
0
]
−
[
0
−
a
12
−
a
13
0
−
a
23
0
]
A=D-L-U=\begin{bmatrix}a_{11}\\&a_{22}\\&&a_{33}\end{bmatrix}-\begin{bmatrix}0\\-a_{21}&0\\-a_{31}&-a_{32}&0\end{bmatrix}-\begin{bmatrix}0&-a_{12}&-a_{13}\\&0&-a_{23}\\&&0\end{bmatrix}
A=D−L−U=
a11a22a33
−
0−a21−a310−a320
−
0−a120−a13−a230
然后,Ax=b ⇒ (D-L-U)x=b ⇒ Dx=(L+U)x+b ⇒ xk+1=D-1(L+U)xk+D-1b。
所以,雅可比迭代法的迭代矩阵为
B
J
=
D
−
1
(
L
+
U
)
=
E
−
D
−
1
A
B_J=D^{-1}(L+U)=E-D^{-1}A
BJ=D−1(L+U)=E−D−1A,常数项为f=D-1b。
Gauss-Seidel迭代法
它是雅可比迭代法的改进版,推导过程为Ax=b ⇒ (D-L-U)x=b ⇒ (D-L)x=Ux+b ⇒xk+1=(D-L)-1Uxk+(D-L)-1b 。
所以,高斯-塞屌迭代法的迭代矩阵为
B
G
S
=
(
D
−
L
)
−
1
U
B_{GS}=(D-L)^{-1}U
BGS=(D−L)−1U,常数项为f=(D-L)-1b。
迭代法收敛性分析
由基本迭代公式得xk=Bxk-1+f,k→∞时有x*=Bx*+f,将两式相减可得xk-x*=B(xk-1-x*)。
将上式循环使用(n-1)次,则可得xk-x*=Bk(x0-x*)。
因此,基本迭代公式收敛等价于
lim
k
→
∞
B
k
=
O
\lim\limits_{k→∞}B^k=O
k→∞limBk=O,而
lim
k
→
∞
B
k
=
O
\lim\limits_{k→∞}B^k=O
k→∞limBk=O等价于ρ(B)<1。
因此,迭代公式收敛的充要条件为其迭代矩阵B满足ρ(B)<1,迭代公式收敛的充分条件为其某种算子范数满足||B||<1。
由ρ(B)≤||B||可知,||B||<1时必有ρ(B)<1,但ρ(B)<1时不一定有||B||<1。
迭代公式的误差估计式:
∣
∣
x
k
−
x
∗
∣
∣
≤
∣
∣
B
∣
∣
1
−
∣
∣
B
∣
∣
∣
∣
x
k
−
x
k
−
1
∣
∣
||x_k-x^*||≤\frac{||B||}{1-||B||}||x_k-x_{k-1}||
∣∣xk−x∗∣∣≤1−∣∣B∣∣∣∣B∣∣∣∣xk−xk−1∣∣或
∣
∣
x
k
−
x
∗
∣
∣
≤
∣
∣
B
∣
∣
k
1
−
∣
∣
B
∣
∣
∣
∣
x
1
−
x
0
∣
∣
||x_k-x^*||≤\frac{||B||^k}{1-||B||}||x_1-x_0||
∣∣xk−x∗∣∣≤1−∣∣B∣∣∣∣B∣∣k∣∣x1−x0∣∣
当某个迭代公式的ρ(B)≥1时称其发散,即“该迭代公式不是对任意初始向量都收敛”。
若方程组的系数矩阵A是严格对角占优矩阵,即
∣
a
i
i
∣
>
∑
j
=
1
,
j
≠
i
n
∣
a
i
j
∣
|a_{ii}|>\sum_{j=1,j≠i}^n|a_{ij}|
∣aii∣>∑j=1,j=in∣aij∣,则解该方程组的J法和GS法均收敛。
若A是不可约弱对角占优矩阵,即
∣
a
i
i
∣
≥
∑
j
=
1
,
j
≠
i
n
∣
a
i
j
∣
|a_{ii}|≥\sum_{j=1,j≠i}^n|a_{ij}|
∣aii∣≥∑j=1,j=in∣aij∣,则解该方程组的J法和GS法均收敛。
若A是对称正定矩阵,则解方程组的GS法收敛。
当基本迭代公式收敛时,由xk-x*=Bk(x0-x*)可得
∣
∣
x
k
−
x
∗
∣
∣
≤
∣
∣
B
∣
∣
k
∣
∣
x
0
−
x
∗
∣
∣
||x_k-x^*||≤||B||^k||x_0-x^*||
∣∣xk−x∗∣∣≤∣∣B∣∣k∣∣x0−x∗∣∣,可见:
(1)迭代初值x0越靠近精确值x*,迭代公式收敛得越快。
(2)在||B||<1的前提下,||B||越小,迭代公式收敛得越快。
(3)由于ρ(B)≤||B||,所以ρ(B)越小,迭代公式收敛得越快。
迭代公式的渐进收敛速度的定义为
R
(
B
)
=
−
ln
ρ
(
B
)
R(B)=-\lnρ(B)
R(B)=−lnρ(B)。
可见,R(B)>0,且R(B)越大,迭代公式收敛得越快。
SOR迭代法
SOR迭代法又叫“逐次超松弛迭代法”,在上述两种方法的基础上引入了松弛因子ω。
在0<ω<1时,这种迭代法为亚松弛算法;ω=1时,它和GS法一致;ω>1时,它成为超松弛算法。
与前文相同,SOR法收敛的充要条件为ρ(Bω)<1,充分条件为||Bω||<1。
不同的是,SOR法收敛的必要条件为0<ω<2。
若A为对称正定矩阵,且0<ω<2,则解对应方程组的SOR法收敛。
若A为对称正定的三对角矩阵,则ρ(BGS)=ρ2(BJ)<1,ρ(Bω)=ω-1,最佳松弛因子为
ω
o
p
t
=
2
1
+
1
−
ρ
(
B
G
S
)
ω_{opt}=\frac2{1+\sqrt{1-ρ(B_{GS})}}
ωopt=1+1−ρ(BGS)2。
常微分方程的数值解法
数值解的概念
初值问题的基本形式为
{
y
′
=
f
(
x
,
y
)
y
(
x
0
)
=
y
0
\left\{\begin{matrix}y'=f(x,y)\\y(x_0)=y_0\end{matrix}\right.
{y′=f(x,y)y(x0)=y0。
上式的解析解y(x)在指定的离散点xi(或指定步长)处的值y(xi)的近似值yi即被称为数值解。
对于一个定解问题,离散点xi或步长由求解者按需确定。在用数值法解初值问题时,通常截取等步长。
微分方程初值问题的数值解法的特点是从已知的y0出发,通过一定的公式依次求出y1、y2、⋯。
Euler方法及其局部截断误差
解初值问题的欧拉公式为
y
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
y_{n+1}=y_n+hf(x_n,y_n)
yn+1=yn+hf(xn,yn),它由一阶向前差分公式
y
n
+
1
−
y
n
h
≈
f
(
x
n
,
y
n
)
\frac{y_{n+1}-y_n}h≈f(x_n,y_n)
hyn+1−yn≈f(xn,yn)演变而来,h为步长。
隐式欧拉公式为
y
n
+
1
=
y
n
+
h
f
(
x
n
+
1
,
y
n
+
1
)
y_{n+1}=y_n+hf(x_{n+1},y_{n+1})
yn+1=yn+hf(xn+1,yn+1),它由一阶向后差分公式
y
n
+
1
−
y
n
h
≈
f
(
x
n
+
1
,
y
n
+
1
)
\frac{y_{n+1}-y_n}h≈f(x_{n+1},y_{n+1})
hyn+1−yn≈f(xn+1,yn+1)演变而来。
梯形公式为
y
n
+
1
=
y
n
+
h
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
n
+
1
)
]
y_{n+1}=y_n+\frac h2[f(x_n,y_n)+f(x_{n+1},y_{n+1})]
yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)],它由上面两式求平均值而得。
改进欧拉公式为
{
y
ˉ
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
y
n
+
1
=
y
n
+
h
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
ˉ
n
+
1
)
]
\left\{\begin{matrix}\bar y_{n+1}=y_n+hf(x_n,y_n)\\y_{n+1}=y_n+\frac h2[f(x_n,y_n)+f(x_{n+1},\bar y_{n+1})]\end{matrix}\right.
{yˉn+1=yn+hf(xn,yn)yn+1=yn+2h[f(xn,yn)+f(xn+1,yˉn+1)],它将隐式的梯形公式用欧拉公式显式化。
单步法的局部截断误差定义为
T
n
+
1
=
y
(
x
n
+
1
)
−
y
n
+
1
=
y
(
x
n
+
1
)
−
y
(
x
n
)
−
h
φ
(
x
n
,
y
(
x
n
)
,
y
(
x
n
+
1
)
,
h
)
T_{n+1}=y(x_{n+1})-y_{n+1}=y(x_{n+1})-y(x_n)-hφ(x_n,y(x_n),y(x_{n+1}),h)
Tn+1=y(xn+1)−yn+1=y(xn+1)−y(xn)−hφ(xn,y(xn),y(xn+1),h)。
由泰勒公式可得:
欧拉公式的局部截断误差为
T
n
+
1
=
1
2
y
′
′
(
x
n
)
h
2
+
O
(
h
3
)
T_{n+1}=\frac12y''(x_n)h^2+O(h^3)
Tn+1=21y′′(xn)h2+O(h3),它具有一阶精度。
隐式欧拉公式的局部截断误差为
T
n
+
1
=
−
1
2
y
′
′
(
x
n
)
h
2
+
O
(
h
3
)
T_{n+1}=-\frac12y''(x_n)h^2+O(h^3)
Tn+1=−21y′′(xn)h2+O(h3),它具有一阶精度。
梯形公式的局部截断误差为
T
n
+
1
=
−
1
12
y
′
′
′
(
x
n
)
h
3
+
O
(
h
4
)
T_{n+1}=-\frac1{12}y'''(x_n)h^3+O(h^4)
Tn+1=−121y′′′(xn)h3+O(h4),它具有二阶精度。
显然,当Tn+1中的第一个非零项含hp+1时,对应公式的截断误差即含有p阶精度,含hp+1的项也被称为局部截断误差主项。
单步法的收敛性与稳定性
解初值问题的显式单步法公式可统一写作
y
n
+
1
=
y
n
+
h
φ
(
x
n
,
y
n
,
h
)
y_{n+1}=y_n+hφ(x_n,y_n,h)
yn+1=yn+hφ(xn,yn,h)。
如果φ中没有xn,且能将上式变成
y
n
+
1
=
f
(
h
)
y
n
y_{n+1}=f(h)y_n
yn+1=f(h)yn,则y的近似数值解为
y
n
=
f
n
(
h
)
y
0
y_n=f^n(h)y_0
yn=fn(h)y0。
对于上式,若对于任一固定的
x
n
=
x
0
+
n
h
x_n=x_0+nh
xn=x0+nh皆有
lim
h
→
0
,
n
→
∞
y
n
=
y
(
x
n
)
\lim\limits_{h→0,n→∞}y_n=y(x_n)
h→0,n→∞limyn=y(xn),则单步法公式收敛于精确解,隐式单步法也一样。
在初值问题中,如果用某一数值方法按某一固定步长h计算时,在节点值yi上产生扰动δi,而且仅由这个扰动引起的yi之后各节点值的变化量的绝对值≤|δi|,则称这个数值方法在取步长h时是绝对稳定的。
研究初值问题的稳定性时,可采用含负实数λ的试验方程/模型方程
y
′
=
λ
y
y'=λy
y′=λy,一般方程总可以线性化为该形式。
将
y
′
=
λ
y
y'=λy
y′=λy代入初值问题欧拉公式
y
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
y_{n+1}=y_n+hf(x_n,y_n)
yn+1=yn+hf(xn,yn)得
y
n
+
1
=
(
λ
h
+
1
)
y
n
y_{n+1}=(λh+1)y_n
yn+1=(λh+1)yn。
根据上式,若yn有一扰动δ,则yn+1的扰动为E(λh)δ=(λh+1)δ。
若初值问题欧拉公式数值稳定,则需误差放大倍数|E(λh)|<1,解得-2<λh<0,故该公式的绝对稳定区间为(-2,0)。
Runge-Kutta方法的原理
解初值问题的Runge-Kutta方法为
{
y
n
+
1
=
y
n
+
h
2
(
k
1
+
k
2
)
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
h
,
y
n
+
h
k
1
)
\left\{\begin{matrix}y_{n+1}=y_n+\frac h2(k_1+k_2)\\k_1=f(x_n,y_n)\\k_2=f(x_n+h,y_n+hk_1)\end{matrix}\right.
⎩
⎨
⎧yn+1=yn+2h(k1+k2)k1=f(xn,yn)k2=f(xn+h,yn+hk1)。
其特点是将
y
′
=
f
(
x
n
,
y
n
)
,
y
′
=
f
(
x
n
+
h
,
y
n
+
h
k
1
)
y'=f(x_n,y_n),y'=f(x_n+h,y_n+hk_1)
y′=f(xn,yn),y′=f(xn+h,yn+hk1)两个斜率加权求和为增量部分。
经典四阶Runge-Kutta公式为
{
y
n
+
1
=
y
n
+
h
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
h
2
,
y
n
+
h
2
k
1
)
k
3
=
f
(
x
n
+
h
2
,
y
n
+
h
2
k
2
)
k
4
=
f
(
x
n
+
h
,
y
n
+
h
k
3
)
\left\{\begin{matrix}y_{n+1}=y_n+\frac h6(k_1+2k_2+2k_3+k_4)\\k_1=f(x_n,y_n)\\k_2=f(x_n+\frac h2,y_n+\frac h2k_1)\\k_3=f(x_n+\frac h2,y_n+\frac h2k_2)\\k_4=f(x_n+h,y_n+hk_3)\end{matrix}\right.
⎩
⎨
⎧yn+1=yn+6h(k1+2k2+2k3+k4)k1=f(xn,yn)k2=f(xn+2h,yn+2hk1)k3=f(xn+2h,yn+2hk2)k4=f(xn+h,yn+hk3),它是一种四级四阶公式。
虽然四阶RK公式精度较高,但计算量比较大,每求解一点需调用函数四次,而欧拉公式和梯形公式只需调用函数一两次。
掌握经典4阶Runge-Kutta公式的特点(性质)
线性多步法的概念
在求解初值问题的多步法中,某一步解的计算不只依赖于前一步解的值,还依赖于前两步(及以上)的解的值。
例如,中心差商公式为
y
(
x
n
+
1
)
−
y
(
x
n
−
1
)
2
h
≈
y
′
(
x
n
)
=
f
(
x
n
,
y
(
x
n
)
)
\frac{y(x_{n+1})-y(x_{n-1})}{2h}≈y'(x_n)=f(x_n,y(x_n))
2hy(xn+1)−y(xn−1)≈y′(xn)=f(xn,y(xn)),变形可得欧拉两步法
y
n
+
1
=
y
n
−
1
+
2
h
f
(
x
n
,
y
n
)
y_{n+1}=y_{n-1}+2hf(x_n,y_n)
yn+1=yn−1+2hf(xn,yn)。
构造线性多步法的例子:
设有两步法公式
y
n
+
1
=
a
0
y
n
+
a
1
y
n
−
1
+
h
(
b
0
f
(
x
n
,
y
n
)
+
b
1
f
(
x
n
−
1
,
y
n
−
1
)
)
y_{n+1}=a_0y_n+a_1y_{n-1}+h(b_0f(x_n,y_n)+b_1f(x_{n-1},y_{n-1}))
yn+1=a0yn+a1yn−1+h(b0f(xn,yn)+b1f(xn−1,yn−1)),为使其精度尽量高,求参数a、b、c、d。
先将其变换成局部截断误差公式
T
n
+
1
=
y
(
x
n
+
1
)
−
a
0
y
(
x
n
)
−
a
1
y
(
x
n
−
1
)
−
b
0
h
y
′
(
x
n
)
−
b
1
h
y
′
(
x
n
−
1
)
T_{n+1}=y(x_{n+1})-a_0y(x_n)-a_1y(x_{n-1})-b_0hy'(x_n)-b_1hy'(x_{n-1})
Tn+1=y(xn+1)−a0y(xn)−a1y(xn−1)−b0hy′(xn)−b1hy′(xn−1)。
然后代入向前和向后差商公式得
T
n
+
1
=
[
1
1
1
2
!
1
3
!
−
a
0
−
a
1
a
1
−
a
1
2
!
a
1
3
!
−
b
0
−
b
1
b
1
−
b
1
2
!
]
[
y
(
x
n
)
h
y
′
(
x
n
)
h
2
y
′
′
(
x
n
)
h
3
y
′
′
′
(
x
n
)
]
+
1
−
a
1
+
4
b
1
4
!
h
4
y
(
4
)
(
x
n
)
+
O
(
h
5
)
T_{n+1}=\begin{bmatrix}1&1&\frac1{2!}&\frac1{3!}\\-a_0\\-a_1&a_1&-\frac {a_1}{2!}&\frac {a_1}{3!}\\&-b_0\\&-b_1&b_1&-\frac{b_1}{2!}\end{bmatrix}\begin{bmatrix}y(x_n)\\hy'(x_n)\\h^2y''(x_n)\\h^3y'''(x_n)\end{bmatrix}+\frac{1-a_1+4b_1}{4!}h^4y^{(4)}(x_n)+O(h^5)
Tn+1=
1−a0−a11a1−b0−b12!1−2!a1b13!13!a1−2!b1
y(xn)hy′(xn)h2y′′(xn)h3y′′′(xn)
+4!1−a1+4b1h4y(4)(xn)+O(h5)。
令h的0~3次项为0,可得方程组
{
1
−
a
0
−
a
1
=
0
1
+
a
1
−
b
0
−
b
1
=
0
1
2
−
a
1
2
+
b
1
=
0
1
6
+
a
1
6
−
b
1
2
=
0
\left\{\begin{matrix}1-a_0-a_1=0\\1+a_1-b_0-b_1=0\\\frac12-\frac{a_1}2+b_1=0\\\frac16+\frac{a_1}6-\frac{b_1}2=0\end{matrix}\right.
⎩
⎨
⎧1−a0−a1=01+a1−b0−b1=021−2a1+b1=061+6a1−2b1=0,解得
a
0
=
−
4
,
a
1
=
5
,
b
0
=
4
,
b
1
=
2
a_0=-4,a_1=5,b_0=4,b_1=2
a0=−4,a1=5,b0=4,b1=2。
所以,所求公式为
y
n
+
1
=
−
4
y
n
+
5
y
n
−
1
+
h
(
4
f
(
x
n
,
y
n
)
+
2
b
1
f
(
x
n
−
1
,
y
n
−
1
)
)
y_{n+1}=-4y_n+5y_{n-1}+h(4f(x_n,y_n)+2b_1f(x_{n-1},y_{n-1}))
yn+1=−4yn+5yn−1+h(4f(xn,yn)+2b1f(xn−1,yn−1)),其局部截断误差为
T
n
+
1
=
1
6
h
4
y
(
4
)
(
x
n
)
+
O
(
h
5
)
T_{n+1}=\frac16h^4y^{(4)}(x_n)+O(h^5)
Tn+1=61h4y(4)(xn)+O(h5),所以该公式具有三阶精度。