误差分析
误差的基本概念
截断误差是指在构造数值计算方法时,用有限的过程代替无限的过程或用简易的方法代替复杂的方法造成的计算结果的误差。
截断误差是由计算的方法造成的,所以也叫方法误差。
舍入误差是计算机的计算精度限制所造成的误差,也叫计算误差;此外,实验的测量精度所造成的误差也归为此类。
绝对误差的定义:e=x-x*,x*在本章为估计值,在后面代表精确值。
绝对误差限的定义:|e|=|x-x*|≤ε,或x=x*±ε
相对误差的定义:
e
r
=
x
−
x
∗
x
e_r=\frac{x-x^*}x
er=xx−x∗,或
e
r
=
x
−
x
∗
x
∗
e_r=\frac{x-x^*}{x^*}
er=x∗x−x∗
相对误差限的定义:
∣
e
r
∣
=
∣
x
−
x
∗
x
∣
≤
ε
r
|e_r|=|\frac{x-x^*}x|≤ε_r
∣er∣=∣xx−x∗∣≤εr,或
ε
r
=
ε
∣
x
∗
∣
ε_r=\fracε{|x^*|}
εr=∣x∗∣ε
函数的误差估计:
e
(
y
)
=
y
′
(
x
∗
)
e
(
x
)
,
e
r
(
y
)
=
e
(
y
)
y
∗
=
y
′
(
x
∗
)
e
(
x
)
y
(
x
∗
)
e(y)=y'(x^*)e(x),e_r(y)=\frac{e(y)}{y^*}=\frac{y'(x^*)e(x)}{y(x^*)}
e(y)=y′(x∗)e(x),er(y)=y∗e(y)=y(x∗)y′(x∗)e(x)
近似数的有效数字的半位判断法:
(1)用3.14近似表示π,因为|3.14-π|≈0.0015926∈[5×10-4,5×10-3],所以0.14对π具有3位有效数字。
(2)用
22
7
\frac{22}7
722近似表示π,因为
∣
22
7
−
π
∣
|\frac{22}7-π|
∣722−π∣≈0.0013∈[5×10-4,5×10-3],所以
22
7
\frac{22}7
722对π具有3位有效数字。
(3)用3.1415近似表示π,因为|3.1415-π|≈0.0000926∈[5×10-5,5×10-4],所以
22
7
\frac{22}7
722对π具有4位有效数字。
(4)用3.1416近似表示π,因为|3.1415-π|≈0.000007∈[5×10-6,5×10-5],所以
22
7
\frac{22}7
722对π具有5位有效数字。
算法的数值稳定性:在迭代时如果误差不会被积累或放大,它就是数值稳定的,当然误差在迭代后变得越小就越好。
做题时一般可得到迭代公式xk+1=f(xk),据此得出误差估计式e(xk+1)=Le(xk),若|L|<1则迭代算法稳定,否则不稳定。线性代数中类似地有ρ(B)<1或||B||<1。
数值算法设计若干原则
(1)尽量简化计算步骤,减少运算次数,尤其是减少乘除法运算次数
多项式y=anxn+⋯+a1x+a0需要
n
(
n
+
1
)
2
\frac{n(n+1)}2
2n(n+1)次乘法运算。
如果将其改写为秦九韶算法的格式,形如y=((a3x+a2)x+a1)x+a0,则只需n次乘法运算。
(2)尽量减少舍入误差的影响,避免有效数字的损失,保证算法的数值稳定性
如通过对公式进行变换避免两相近数相减、大数吃掉小数或分母极小等情形。
函数的插值拟合法
插值问题的提法
理解插值问题的基本概念、插值多项式的存在唯一性。
Lagrange插值法
若知道n个点x0、x1、⋯、xn与对应的函数值f0、f1、⋯、fn,就可以据此作出一个唯一的n次拉格朗日插值多项式。
设lk(x)为拉格朗日插值基函数,则:
n=1时,拉格朗日线性插值公式为
L
1
(
x
)
=
l
0
f
0
+
l
1
f
1
=
x
−
x
1
x
0
−
x
1
f
0
+
x
−
x
0
x
1
−
x
0
f
1
L_1(x)=l_0f_0+l_1f_1=\frac{x-x_1}{x_0-x_1}f_0+\frac{x-x_0}{x_1-x_0}f_1
L1(x)=l0f0+l1f1=x0−x1x−x1f0+x1−x0x−x0f1。
n=2时,拉格朗日抛物线插值公式为
L
2
(
x
)
=
∑
l
i
f
i
=
(
x
−
x
1
)
(
x
−
x
2
)
(
x
0
−
x
1
)
(
x
0
−
x
2
)
f
0
+
(
x
−
x
0
)
(
x
−
x
2
)
(
x
1
−
x
0
)
(
x
1
−
x
2
)
f
1
+
(
x
−
x
0
)
(
x
−
x
1
)
(
x
2
−
x
0
)
(
x
2
−
x
1
)
f
2
L_2(x)=∑l_if_i=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f_0+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f_1+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f_2
L2(x)=∑lifi=(x0−x1)(x0−x2)(x−x1)(x−x2)f0+(x1−x0)(x1−x2)(x−x0)(x−x2)f1+(x2−x0)(x2−x1)(x−x0)(x−x1)f2。
高阶插值公式同理。此外,在拉格朗日插值中,有
[
l
0
(
x
0
)
l
0
(
x
1
)
l
1
(
x
2
)
l
1
(
x
0
)
l
1
(
x
1
)
l
1
(
x
2
)
l
2
(
x
0
)
l
2
(
x
1
)
l
2
(
x
2
)
]
=
[
1
1
1
]
=
E
\begin{bmatrix}l_0(x_0)&l_0(x_1)&l_1(x_2)\\l_1(x_0)&l_1(x_1)&l_1(x_2)\\l_2(x_0)&l_2(x_1)&l_2(x_2)\end{bmatrix}=\begin{bmatrix}1\\&1\\&&1\end{bmatrix}=E
⎣⎡l0(x0)l1(x0)l2(x0)l0(x1)l1(x1)l2(x1)l1(x2)l1(x2)l2(x2)⎦⎤=⎣⎡111⎦⎤=E。
原函数在[a,b]上可以表示成f(x)=Ln(x)+Rn(x),即拉格朗日插值多项式及其余项之和。
拉格朗日插值的余项公式:
R
n
(
x
)
=
f
(
n
+
1
)
(
ξ
)
(
n
+
1
)
!
ω
(
n
+
1
)
(
x
)
=
f
(
n
+
1
)
(
ξ
)
(
n
+
1
)
!
∏
i
=
0
n
(
x
−
x
i
)
R_n(x)=\frac{f^{(n+1)}(ξ)}{(n+1)!}ω_{(n+1)}(x)=\frac{f^{(n+1)}(ξ)}{(n+1)!}∏_{i=0}^n(x-x_i)
Rn(x)=(n+1)!f(n+1)(ξ)ω(n+1)(x)=(n+1)!f(n+1)(ξ)∏i=0n(x−xi),其中ξ∈(a,b)且与x有关。
Hermite插值法
若已知两点x0、x1,对应的函数值f0、f1和导数值f’0、f’1,则相应的两点三次Hermite插值公式为:
α
0
(
x
)
=
(
1
+
2
x
−
x
0
x
1
−
x
0
)
(
x
−
x
1
x
0
−
x
1
)
2
α_0(x)=(1+2\frac{x-x_0}{x_1-x_0})(\frac{x-x_1}{x_0-x_1})^2
α0(x)=(1+2x1−x0x−x0)(x0−x1x−x1)2
α
1
(
x
)
=
(
1
+
2
x
−
x
1
x
0
−
x
1
)
(
x
−
x
0
x
1
−
x
0
)
2
α_1(x)=(1+2\frac{x-x_1}{x_0-x_1})(\frac{x-x_0}{x_1-x_0})^2
α1(x)=(1+2x0−x1x−x1)(x1−x0x−x0)2
β
0
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
x
0
−
x
1
)
2
β_0(x)=(x-x_0)(\frac{x-x_1}{x_0-x_1})^2
β0(x)=(x−x0)(x0−x1x−x1)2
β
1
(
x
)
=
(
x
−
x
1
)
(
x
−
x
0
x
1
−
x
0
)
2
β_1(x)=(x-x_1)(\frac{x-x_0}{x_1-x_0})^2
β1(x)=(x−x1)(x1−x0x−x0)2
H
3
(
x
)
=
α
0
f
0
+
α
1
f
1
+
β
0
f
0
′
+
β
1
f
1
′
H_3(x)=α_0f_0+α_1f_1+β_0f'_0+β_1f'_1
H3(x)=α0f0+α1f1+β0f0′+β1f1′
其余项表达式为
R
3
(
x
)
=
f
(
x
)
−
H
3
(
x
)
=
f
(
4
)
4
!
(
x
−
x
0
)
2
(
x
−
x
1
)
2
R_3(x)=f(x)-H_3(x)=\frac{f^{(4)}}{4!}(x-x_0)^2(x-x_1)^2
R3(x)=f(x)−H3(x)=4!f(4)(x−x0)2(x−x1)2。
若是知道三点的f0、f1、f2和f’1,可先作三点的二次拉格朗日插值函数L2,再设H3=L2+k(x-x0)(x-x1)(x-x2),解出k即可得H3。
也可直接设H3(x)=a3x3+a2x2+a1x+a0,然后暴力代入求解。
最后,其余项为
R
3
(
x
)
=
f
(
4
)
4
!
(
x
−
x
0
)
(
x
−
x
1
)
2
(
x
−
x
2
)
R_3(x)=\frac{f^{(4)}}{4!}(x-x_0)(x-x_1)^2(x-x_2)
R3(x)=4!f(4)(x−x0)(x−x1)2(x−x2)。
差商
一阶差商公式:
f
[
x
0
,
x
1
]
=
f
(
x
0
)
−
f
(
x
1
)
x
0
−
x
1
f[x_0,x_1]=\frac{f(x_0)-f(x_1)}{x_0-x_1}
f[x0,x1]=x0−x1f(x0)−f(x1)。若设零阶差商为函数值,分子也可看做两个零阶差商的差。
二阶差商公式:
f
[
x
0
,
x
1
,
x
2
]
=
f
[
x
0
,
x
1
]
−
f
[
x
1
,
x
2
]
x
0
−
x
2
f[x_0,x_1,x_2]=\frac{f[x_0,x_1]-f[x_1,x_2]}{x_0-x_2}
f[x0,x1,x2]=x0−x2f[x0,x1]−f[x1,x2]
n阶差商公式:
f
[
x
0
,
x
1
,
⋯
,
x
n
]
=
f
[
x
0
,
x
1
,
⋯
,
x
n
−
1
]
−
f
[
x
1
,
x
2
,
⋯
,
x
n
]
x
1
−
x
n
f[x_0,x_1,⋯,x_n]=\frac{f[x_0,x_1,⋯,x_{n-1}]-f[x_1,x_2,⋯,x_n]}{x_1-x_n}
f[x0,x1,⋯,xn]=x1−xnf[x0,x1,⋯,xn−1]−f[x1,x2,⋯,xn]
求差商时需要做阶梯型差商表。
差商的性质:在[a,b]中求差商时,有
f
[
x
0
,
x
1
,
⋯
,
x
n
]
=
f
(
n
)
(
ξ
)
n
!
f[x_0,x_1,⋯,x_n]=\frac{f^{(n)}(ξ)}{n!}
f[x0,x1,⋯,xn]=n!f(n)(ξ),其中ξ∈[a,b]。
Newton插值法
牛顿插值多项式是利用数据点的各阶差商拟合的多项式。
一阶牛顿插值多项式:
N
1
=
f
(
x
0
)
+
f
[
x
0
,
x
1
]
(
x
−
x
0
)
N_1=f(x_0)+f[x_0,x_1](x-x_0)
N1=f(x0)+f[x0,x1](x−x0)
二阶牛顿插值多项式:
N
2
=
f
(
x
0
)
+
f
[
x
0
,
x
1
]
(
x
−
x
0
)
+
f
[
x
0
,
x
1
,
x
2
]
(
x
−
x
0
)
(
x
−
x
1
)
N_2=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)
N2=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)
n阶牛顿插值多项式:
N
n
=
N
n
−
1
+
f
[
x
0
,
x
1
,
⋯
,
x
n
]
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
n
−
1
)
N_n=N_{n-1}+f[x_0,x_1,⋯,x_n](x-x_0)(x-x_1)⋯(x-x_{n-1})
Nn=Nn−1+f[x0,x1,⋯,xn](x−x0)(x−x1)⋯(x−xn−1)
牛顿插值多项式的余项估计:
∣
R
n
(
x
)
∣
≤
∣
f
[
x
0
,
x
1
,
⋯
,
x
n
,
x
n
+
1
]
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
n
)
∣
|R_n(x)|≤|f[x_0,x_1,⋯,x_n,x_{n+1}](x-x_0)(x-x_1)⋯(x-x_n)|
∣Rn(x)∣≤∣f[x0,x1,⋯,xn,xn+1](x−x0)(x−x1)⋯(x−xn)∣
牛顿插值法比拉格朗日插值法节省3~4倍的工作量。
分段线性插值法
Runge现象即用三次以上的高次多项式拟合数据时产生的过拟合现象。
为了避免数据过多时高次多项式产生的过拟合现象,可采用分段低次插值法拟合函数。
分段线性插值法即分段使用一次拉格朗日插值公式,虽然表达式写作φ(x)=∑lifi,但li只在前后区间中起作用,而在其它区间内为0。
三次样条函数与三次样条插值法
所谓三次样条曲线就是挠曲线,即光滑的分段三次多项式。
挠曲线中的挠度、转角和弯矩连续,因此三次样条函数的f(x)、f’(x)、f"(x)也连续。剪力不一定连续,所以f"’(x)不一定连续。
特别地,自然边界条件指挠曲线的最左端和最右端处有f"(x)=0,即弯矩为0。
如果将挠曲线用于拟合周期函数,则需保证曲线头尾的挠度、转角和弯矩相等。
函数的最小二乘拟合法
曲线拟合的最小二乘法
对离散实验数据做最小二乘拟合的两个主要步骤是选定拟合模型的形式和求模型的最小二乘解。
曲线拟合的最小二乘法是指用偏导数求“曲线上各点的拟合值与近似值的误差平方和最小”时的曲线参数。
这也就是机器学习里经典的用伪逆求超定方程组
A
x
=
b
Ax=b
Ax=b的近似解
x
=
(
A
T
A
)
−
1
A
T
b
x=(A^TA)^{-1}A^Tb
x=(ATA)−1ATb的过程。
直接用最小二乘法的定义求,即列出误差函数e=∑(yi-fi)2后,对拟合函数的各参数列
∂
e
∂
k
i
=
0
\frac{∂ e}{∂ k_i} =0
∂ki∂e=0,计算参数比法方程方便。
在上述解方程过程中所使用的方程
A
T
A
x
=
A
T
b
A^TAx=A^Tb
ATAx=ATb被称为法方程。
若把A看做列向量集[φ0,φ1,⋯,φn],那么法方程的系数矩阵G=ATA的元素aij即为对应列向量的内积(φi,φj)。
正交多项式
在使用三次及以上的多项式拟合函数时,法方程的条件数会变得很大,使之变得病态。
因此,需要用“基于正交多项式的曲线拟合方法”改善这一病态现象。
函数f(x)、g(x)在区间[a,b]上的权重为ρ的函数内积被定义为
(
f
,
g
)
=
∫
a
b
ρ
(
x
)
f
(
x
)
g
(
x
)
d
x
(f,g)=∫_a^bρ(x)f(x)g(x)\mathrm dx
(f,g)=∫abρ(x)f(x)g(x)dx。
在区间[a,b]内,如果函数f(x)、g(x)和权函数ρ(x)满足(f,g)=0,则称函数f与g在[a,b]上带权ρ正交。
假设有一函数序列{φ0,φ1,⋯,φn},其中的每个元素都是一个函数。
在某一区间内,若该序列中任意两个不同的多项式带权ρ正交,则该函数序列带权ρ正交。
在此基础上,若序列中每个元素与自己的(ρ=1的)内积都为1,则该序列为规范正交序列。
假设有多项式序列{φ0,φ1,⋯,φn},其中φi为i次多项式。
在某一区间内,若该序列中任意两个不同的多项式带权ρ正交,则该多项式序列带权ρ正交。
线性代数中将向量组A变为正交向量组B的施密特正交化公式为
{
β
1
=
α
1
β
2
=
α
2
−
(
α
2
,
β
1
)
(
β
1
,
β
1
)
⋅
β
1
β
3
=
α
3
−
(
α
3
,
β
1
)
(
β
1
⋅
β
1
)
⋅
β
1
−
(
α
3
,
β
2
)
(
β
2
,
β
2
)
⋅
β
2
⋯
\left\{\begin{matrix}β_1=α_1\\β_2=α_2-\frac{(α_2,β_1)}{(β_1,β_1)}·β_1\\β_3=α_3-\frac{(α_3,β_1)}{(β_1·β_1)}·β_1-\frac{(α_3,β_2)}{(β_2,β_2)}·β_2\\⋯\end{matrix}\right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧β1=α1β2=α2−(β1,β1)(α2,β1)⋅β1β3=α3−(β1⋅β1)(α3,β1)⋅β1−(β2,β2)(α3,β2)⋅β2⋯
同理可得Gram-Schmidt正交化方法的公式为
{
φ
0
=
1
φ
1
=
x
−
(
x
,
φ
0
)
(
φ
0
,
φ
0
)
⋅
φ
0
φ
2
=
x
2
−
(
x
2
,
φ
0
)
(
φ
0
,
φ
0
)
⋅
φ
0
−
(
x
2
,
φ
1
)
(
φ
1
,
φ
1
)
⋅
φ
1
⋯
\left\{\begin{matrix}φ_0=1\\φ_1=x-\frac{(x,φ_0)}{(φ_0,φ_0)}·φ_0\\φ_2=x^2-\frac{(x^2,φ_0)}{(φ_0,φ_0)}·φ_0-\frac{(x^2,φ_1)}{(φ_1,φ_1)}·φ_1\\⋯\end{matrix}\right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧φ0=1φ1=x−(φ0,φ0)(x,φ0)⋅φ0φ2=x2−(φ0,φ0)(x2,φ0)⋅φ0−(φ1,φ1)(x2,φ1)⋅φ1⋯
Legendre多项式的定义为
P
0
(
x
)
=
1
,
P
n
(
x
)
=
1
2
n
!
d
n
d
x
n
[
(
x
2
−
1
)
n
]
,
x
∈
[
−
1
,
1
]
P_0(x)=1,P_n(x)=\frac1{2n!}\frac{\mathrm d^n}{\mathrm dx^n}[(x^2-1)^n],x∈[-1,1]
P0(x)=1,Pn(x)=2n!1dxndn[(x2−1)n],x∈[−1,1],其性质为:
(1)正交性:两个不同次数的勒让德多项式(带权ρ=1)正交。
(2)递推公式:
P
0
(
x
)
=
1
,
P
1
(
x
)
=
x
,
P
n
+
1
(
x
)
=
2
n
+
1
n
+
1
x
P
n
(
x
)
−
n
n
+
1
P
n
−
1
(
x
)
P_0(x)=1,P_1(x)=x,P_{n+1}(x)=\frac{2n+1}{n+1}xP_n(x)-\frac n{n+1}P_{n-1}(x)
P0(x)=1,P1(x)=x,Pn+1(x)=n+12n+1xPn(x)−n+1nPn−1(x)
(3)奇偶性:当勒让德多项式的最高次数为奇数时,多项式为奇函数,最高次数为偶数时为偶函数
(4)勒让德多项式的所有根都是单根,且与原点对称。
连续函数的最佳平方逼近
Span代表线性空间,Φ=Span(1,x,x2,⋯ ,xn)代表一个线性多项式空间。比如,p(x)∈P=Span(1,x2)等价于p(x)=a0+a1x2。
为了求区间[a,b]上对f(x)的最佳平方逼近函数g(x),需让区间平方误差
e
=
∫
a
b
[
f
(
x
)
−
g
(
x
)
]
2
d
x
e=∫_a^b[f(x)-g(x)]^2\mathrm dx
e=∫ab[f(x)−g(x)]2dx最小。
为了让e最小,自然要让g内各参数满足
∂
e
∂
k
i
=
0
\frac{∂ e}{∂ k_i}=0
∂ki∂e=0,从而求出各参数。
一般考试中,要求的是积分比较好算的最佳平方逼近多项式。
数值法求定积分
数值求积的基本思想、插值型求积公式与代数精度
定积分的定义式为
∫
a
b
f
(
x
)
d
x
=
lim
λ
→
∞
∑
Δ
x
k
f
(
ξ
k
)
∫_a^bf(x)\mathrm dx=\lim\limits_{λ→∞}∑Δx_kf(ξ_k)
∫abf(x)dx=λ→∞lim∑Δxkf(ξk)。
若积分的每个小区间很小,但不是无限小,则上式可改写为
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
A
k
f
(
x
k
)
+
R
n
[
f
]
∫_a^bf(x)\mathrm dx=∑_{k=0}^nA_kf(x_k)+R_n[f]
∫abf(x)dx=∑k=0nAkf(xk)+Rn[f]。
上式被称为函数f在[a,b]上的数值求积公式及其余项,xk为求积节点,Ak为与被积函数无关的求积系数。
拉格朗日插值公式为
f
(
x
)
=
∑
i
=
0
n
l
i
(
x
)
f
i
+
f
(
n
+
1
)
(
ξ
)
(
n
+
1
)
!
ω
(
n
+
1
)
(
x
)
f(x)=∑_{i=0}^nl_i(x)f_i+\frac{f^{(n+1)}(ξ)}{(n+1)!}ω_{(n+1)}(x)
f(x)=∑i=0nli(x)fi+(n+1)!f(n+1)(ξ)ω(n+1)(x)
将上式两边积分,即可得插值型求积公式
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
f
k
∫
a
b
l
k
(
x
)
d
x
+
1
(
n
+
1
)
!
∫
a
b
f
(
n
+
1
)
(
ξ
(
x
)
)
ω
(
n
+
1
)
(
x
)
d
x
∫_a^bf(x)\mathrm dx=∑_{k=0}^nf_k∫_a^bl_k(x)\mathrm dx+\frac1{(n+1)!}∫_a^bf^{(n+1)}(ξ(x))ω_{(n+1)}(x)\mathrm dx
∫abf(x)dx=∑k=0nfk∫ablk(x)dx+(n+1)!1∫abf(n+1)(ξ(x))ω(n+1)(x)dx。
因此,插值型求积公式中的求积系数为
∫
a
b
l
k
(
x
)
d
x
∫_a^bl_k(x)\mathrm dx
∫ablk(x)dx,余项为
1
(
n
+
1
)
!
∫
a
b
f
(
n
+
1
)
(
ξ
(
x
)
)
ω
(
n
+
1
)
(
x
)
d
x
\frac1{(n+1)!}∫_a^bf^{(n+1)}(ξ(x))ω_{(n+1)}(x)\mathrm dx
(n+1)!1∫abf(n+1)(ξ(x))ω(n+1)(x)dx。
若插值型求积公式中的被积函数为次数≤n的多项式,因为其余项在求(n+1)次导后会变为0,所以插值型求积公式会变为准确求积公式。(多项式不是直接硬积就完事了嘛)
对于数值求积公式
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
A
k
f
(
x
k
)
+
R
n
[
f
]
∫_a^bf(x)\mathrm dx=∑_{k=0}^nA_kf(x_k)+R_n[f]
∫abf(x)dx=∑k=0nAkf(xk)+Rn[f],若f(x)为m次及以下的多项式时公式的余项为0,f(x)为(m+1)次多项式时余项不为0,则该数值求积公式具有m次代数精度。
用(n+1)个结点求得的插值型求积公式至少具有n次代数精度,在做题时可以将f(x)依次代为
1
,
x
,
x
2
,
⋯
,
x
n
1,x,x^2,⋯,x^n
1,x,x2,⋯,xn求积分系数,并代入更高次式验算所得公式是否有更高的代数精度。
Newton-Cotes型求积公式(等距节点插值型求积公式)
Newton-Cotes型求积公式的通式为
∫
a
b
f
(
x
)
d
x
=
(
b
−
a
)
∑
k
=
0
n
C
k
(
n
)
f
(
x
k
)
+
R
[
f
]
∫_a^bf(x)\mathrm dx=(b-a)∑_{k=0}^nC_k^{(n)}f(x_k)+R[f]
∫abf(x)dx=(b−a)∑k=0nCk(n)f(xk)+R[f]。
其中,
C
k
(
n
)
=
(
−
1
)
n
−
k
n
k
!
(
n
−
k
)
!
∫
0
n
t
(
t
−
1
)
(
t
−
2
)
⋯
(
t
−
n
)
d
t
C_k^{(n)}=\frac{(-1)^{n-k}}{nk!(n-k)!}∫_0^nt(t-1)(t-2)⋯(t-n)\mathrm dt
Ck(n)=nk!(n−k)!(−1)n−k∫0nt(t−1)(t−2)⋯(t−n)dt被称为Cotes系数。
-
用两个端点求积时,设Newton-Cotes求积公式为 ∫ − 1 1 f ( x ) d x ≈ A f ( − 1 ) + B f ( 1 ) ∫_{-1}^1f(x)\mathrm dx≈Af(-1)+Bf(1) ∫−11f(x)dx≈Af(−1)+Bf(1)。
为了使求积公式精度最高,列出 { ∫ − 1 1 d x = A + B = 2 ∫ − 1 1 x d x = − A + B = 0 \left\{\begin{matrix}∫_{-1}^1\mathrm dx=A+B=2\\∫_{-1}^1x\mathrm dx=-A+B=0\end{matrix}\right. {∫−11dx=A+B=2∫−11xdx=−A+B=0,解得 A = 1 , B = 1 A=1,B=1 A=1,B=1。
将上述系数除以区间长度,即可得一阶Cotes系数为 1 2 , 1 2 \frac12,\frac12 21,21。
因此,梯形公式(一阶Newton-Cotes求积公式)为 ∫ a b f ( x ) d x = ( b − a ) [ 1 2 f ( a ) + 1 2 f ( b ) ] − ( b − a ) 3 12 f ′ ′ ( ξ ) , ξ ∈ [ a , b ] ∫_a^bf(x)\mathrm dx=(b-a)[\frac12f(a)+\frac12f(b)]-\frac{(b-a)^3}{12}f''(ξ),ξ∈[a,b] ∫abf(x)dx=(b−a)[21f(a)+21f(b)]−12(b−a)3f′′(ξ),ξ∈[a,b]。
显然,将f(x)代为x2时上述公式左右两边不等,所以一阶Newton-Cotes求积公式具有一次代数精度。 -
用三个等距节点求积时,设Newton-Cotes求积公式为 ∫ − 1 1 f ( x ) d x ≈ A f ( − 1 ) + B f ( 0 ) + C f ( 1 ) ∫_{-1}^1f(x)\mathrm dx≈Af(-1)+Bf(0)+Cf(1) ∫−11f(x)dx≈Af(−1)+Bf(0)+Cf(1)。
为了使求积公式精度最高,列出 { ∫ − 1 1 d x = A + B + C = 2 ∫ − 1 1 x d x = − A + C = 0 ∫ − 1 1 x 2 d x = A + C = 2 3 \left\{\begin{matrix}∫_{-1}^1\mathrm dx=A+B+C=2\\∫_{-1}^1x\mathrm dx=-A+C=0\\∫_{-1}^1x^2\mathrm dx=A+C=\frac23\end{matrix}\right. ⎩⎪⎨⎪⎧∫−11dx=A+B+C=2∫−11xdx=−A+C=0∫−11x2dx=A+C=32,解得 A = 1 3 , B = 4 3 , C = 1 3 A=\frac13,B=\frac43,C=\frac13 A=31,B=34,C=31。
将上述系数除以区间长度,即可得二阶Cotes系数为 1 6 , 4 6 , 1 6 \frac16,\frac46,\frac16 61,64,61。
因此,Simpson公式(二阶Newton-Cotes求积公式)为 ∫ a b f ( x ) d x = ( b − a ) [ 1 6 f ( a ) + 4 6 f ( a + b 2 ) + 1 6 f ( b ) ] − ( b − a ) 5 2880 f ( 4 ) ( ξ ) , ξ ∈ [ a , b ] ∫_a^bf(x)\mathrm dx=(b-a)[\frac16f(a)+\frac46f(\frac{a+b}2)+\frac16f(b)]-\frac{(b-a)^5}{2880}f^{(4)}(ξ),ξ∈[a,b] ∫abf(x)dx=(b−a)[61f(a)+64f(2a+b)+61f(b)]−2880(b−a)5f(4)(ξ),ξ∈[a,b]。
显然,将f(x)代为x3时上述公式左右两边也相等,代入x4时则不等,所以一阶Newton-Cotes求积公式具有三次代数精度。
以此类推,奇数阶Newton-Cotes型求积公式的代数精度=阶数,偶数阶公式的代数精度=阶数+1。
求余项系数的方法:
例如
∫
0
1
f
(
x
)
d
x
=
3
4
f
(
1
3
)
+
1
4
f
(
1
)
+
C
f
′
′
′
(
ξ
)
∫_0^1f(x)\mathrm dx=\frac34f(\frac13)+\frac14f(1)+Cf'''(ξ)
∫01f(x)dx=43f(31)+41f(1)+Cf′′′(ξ)二阶收敛,将f(x)代入为x3,原式变为
1
4
=
1
36
+
1
4
+
6
C
\frac14=\frac1{36}+\frac14+6C
41=361+41+6C,故
C
=
−
1
216
C=-\frac1{216}
C=−2161。
令插值型求积公式
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
A
k
f
(
x
k
)
+
R
n
[
f
]
∫_a^bf(x)\mathrm dx=∑_{k=0}^nA_kf(x_k)+R_n[f]
∫abf(x)dx=∑k=0nAkf(xk)+Rn[f]中的f(x)=1,可得系数之和
∑
k
=
0
n
A
k
=
b
−
a
∑_{k=0}^nA_k=b-a
∑k=0nAk=b−a。
显然,同阶Cotes系数不仅是对称的,且它们的和为1。阶数够高时,Cotes系数会为负。
若插值型求积公式
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
A
k
f
(
x
k
)
+
R
n
[
f
]
∫_a^bf(x)\mathrm dx=∑_{k=0}^nA_kf(x_k)+R_n[f]
∫abf(x)dx=∑k=0nAkf(xk)+Rn[f]中的求积系数都>0,则该求积公式是数值稳定的。
与分段低次插值类似,用分段低次插值求积法可以提高积分精度,这也被称为复化求积方法。
将[a,b]等分为长度为h的n个区间,用(n+1)个区间点进行插值积分,则:
复化梯形公式及其余项:
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
−
1
h
[
1
2
f
(
x
k
)
+
1
2
f
(
x
k
+
1
)
]
−
(
b
−
a
)
h
2
12
f
′
′
(
ξ
)
,
ξ
∈
(
a
,
b
)
∫_a^bf(x)\mathrm dx=∑_{k=0}^{n-1}h[\frac12f(x_k)+\frac12f(x_{k+1})]-\frac{(b-a)h^2}{12}f''(ξ),ξ∈(a,b)
∫abf(x)dx=∑k=0n−1h[21f(xk)+21f(xk+1)]−12(b−a)h2f′′(ξ),ξ∈(a,b)
复化Simpson公式及其余项:
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
−
1
h
[
1
6
f
(
x
k
)
+
4
6
f
(
x
k
+
x
k
+
1
2
)
+
1
6
f
(
x
k
+
1
)
]
−
(
b
−
a
)
h
4
2880
f
(
4
)
(
ξ
)
,
ξ
∈
(
a
,
b
)
∫_a^bf(x)\mathrm dx=∑_{k=0}^{n-1}h[\frac16f(x_k)+\frac46f(\frac{x_k+x_{k+1}}2)+\frac16f(x_{k+1})]-\frac{(b-a)h^4}{2880}f^{(4)}(ξ),ξ∈(a,b)
∫abf(x)dx=∑k=0n−1h[61f(xk)+64f(2xk+xk+1)+61f(xk+1)]−2880(b−a)h4f(4)(ξ),ξ∈(a,b)
Gauss型求积公式
Gauss型求积公式是指形如
∫
a
b
ρ
(
x
)
f
(
x
)
d
x
=
∑
k
=
0
n
A
k
f
(
x
k
)
+
R
n
[
f
]
∫_a^bρ(x)f(x)\mathrm dx=∑_{k=0}^nA_kf(x_k)+R_n[f]
∫abρ(x)f(x)dx=∑k=0nAkf(xk)+Rn[f]的带权ρ的、具有(2n+1)次代数精度的插值求积公式。
高斯型求积公式的插值节点被称为高斯点。
高斯点等价于在[a,b]上使ωn+1(x)和任何次数≤n的多项式p(x)带权ρ正交的点,写成公式就是
∫
a
b
ρ
(
x
)
ω
n
+
1
(
x
)
p
(
x
)
d
x
=
0
∫_a^bρ(x)ω_{n+1}(x)p(x)\mathrm dx=0
∫abρ(x)ωn+1(x)p(x)dx=0。
当x∈[-1,1]时,令权 ρ = 1 ρ=1 ρ=1,即得Gauss-Legendre求积公式 ∫ − 1 1 f ( x ) d x = ∑ k = 0 n A k f ( x k ) + R n [ f ] ∫_{-1}^1f(x)\mathrm dx=∑_{k=0}^nA_kf(x_k)+R_n[f] ∫−11f(x)dx=∑k=0nAkf(xk)+Rn[f]。
-
用两个节点求积时,设Gauss-Legendre求积公式为 ∫ − 1 1 f ( x ) d x ≈ A f ( − α ) + B f ( α ) ∫_{-1}^1f(x)\mathrm dx≈Af(-α)+Bf(α) ∫−11f(x)dx≈Af(−α)+Bf(α)。
为了使求积公式精度最高,列出 { ∫ − 1 1 d x = A + B = 2 ∫ − 1 1 x d x = − A α + B α = 0 ∫ − 1 1 x 2 d x = A α 2 + B α 2 = 2 3 \left\{\begin{matrix}∫_{-1}^1\mathrm dx=A+B=2\\∫_{-1}^1x\mathrm dx=-Aα+Bα=0\\∫_{-1}^1x^2\mathrm dx=Aα^2+Bα^2=\frac23\end{matrix}\right. ⎩⎪⎨⎪⎧∫−11dx=A+B=2∫−11xdx=−Aα+Bα=0∫−11x2dx=Aα2+Bα2=32,解得 A = 1 , B = 1 , α = 3 3 A=1,B=1,α=\frac{\sqrt3}3 A=1,B=1,α=33。
因此,两点Gauss-Legendre求积公式为 ∫ − 1 1 f ( x ) d x ≈ f ( − 3 3 ) + f ( 3 3 ) ∫_{-1}^1f(x)\mathrm dx≈f(-\frac{\sqrt3}3)+f(\frac{\sqrt3}3) ∫−11f(x)dx≈f(−33)+f(33),两个求积结点即为高斯点。
显然,将f(x)代为x3时上述公式左右两边也相等,代入x4时则不等,所以一阶Gauss-Legendre求积公式具有三次代数精度。 -
用三个节点求积时,设Gauss-Legendre求积公式为 ∫ − 1 1 f ( x ) d x ≈ A f ( − α ) + B f ( 0 ) + C f ( α ) ∫_{-1}^1f(x)\mathrm dx≈Af(-α)+Bf(0)+Cf(α) ∫−11f(x)dx≈Af(−α)+Bf(0)+Cf(α)。
为了使求积公式精度最高,列出 { ∫ − 1 1 d x = A + B + C = 2 ∫ − 1 1 x d x = − A α + C α = 0 ∫ − 1 1 x 2 d x = A α 2 + C α 2 = 2 3 ∫ − 1 1 x 3 d x = − A α 3 + C α 3 = 0 ∫ − 1 1 x 4 d x = A α 4 + C α 4 = 2 5 \left\{\begin{matrix}∫_{-1}^1\mathrm dx=A+B+C=2\\∫_{-1}^1x\mathrm dx=-Aα+Cα=0\\∫_{-1}^1x^2\mathrm dx=Aα^2+Cα^2=\frac23\\∫_{-1}^1x^3\mathrm dx=-Aα^3+Cα^3=0\\∫_{-1}^1x^4\mathrm dx=Aα^4+Cα^4=\frac25\end{matrix}\right. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧∫−11dx=A+B+C=2∫−11xdx=−Aα+Cα=0∫−11x2dx=Aα2+Cα2=32∫−11x3dx=−Aα3+Cα3=0∫−11x4dx=Aα4+Cα4=52,解得 A = 5 9 , B = 8 9 , C = 5 9 , α = 15 5 A=\frac59,B=\frac89,C=\frac59,α=\frac{\sqrt{15}}5 A=95,B=98,C=95,α=515。
因此,三点Gauss-Legendre求积公式为 ∫ − 1 1 f ( x ) d x ≈ 5 9 f ( − 15 5 ) + 8 9 f ( 0 ) + 5 9 f ( 15 5 ) ∫_{-1}^1f(x)\mathrm dx≈\frac59f(-\frac{\sqrt{15}}5)+\frac89f(0)+\frac59f(\frac{\sqrt{15}}5) ∫−11f(x)dx≈95f(−515)+98f(0)+95f(515),三个求积结点即为高斯点。
显然,将f(x)代为x5时上述公式左右两边也相等,代为x6时则不等,所以二阶Gauss-Legendre求积公式具有五次代数精度。
以此类推,n阶Gauss-Legendre求积公式具有(2n+1)次代数精度,其实高斯型求积公式都具有(2n+1)次代数精度。
当x∈[-1,1]时,令权
ρ
=
1
1
−
x
2
ρ=\frac1{\sqrt{1-x^2}}
ρ=1−x21,即得Gauss-Chebyshev求积公式
∫
−
1
1
f
(
x
)
1
−
x
2
d
x
=
∑
k
=
0
n
A
k
f
(
x
k
)
+
R
n
[
f
]
∫_{-1}^1\frac{f(x)}{\sqrt{1-x^2}}\mathrm dx=∑_{k=0}^nA_kf(x_k)+R_n[f]
∫−111−x2f(x)dx=∑k=0nAkf(xk)+Rn[f]。
其系数为
A
k
=
π
n
+
1
A_k=\fracπ{n+1}
Ak=n+1π,故求积公式为
∫
−
1
1
f
(
x
)
1
−
x
2
d
x
≈
π
n
+
1
∑
k
=
0
n
f
(
x
k
)
∫_{-1}^1\frac{f(x)}{\sqrt{1-x^2}}\mathrm dx≈\fracπ{n+1}∑_{k=0}^nf(x_k)
∫−111−x2f(x)dx≈n+1π∑k=0nf(xk)。
令高斯型求积公式的
f
(
x
)
=
l
k
2
(
x
)
f(x)=l_k^2(x)
f(x)=lk2(x),即拉格朗日插值基函数的平方,它是2n次多项式,故
∫
a
b
ρ
(
x
)
l
k
2
(
x
)
d
x
=
∑
i
=
0
n
A
i
l
k
2
(
x
i
)
∫_a^bρ(x)l_k^2(x)\mathrm dx=∑_{i=0}^nA_il_k^2(x_i)
∫abρ(x)lk2(x)dx=∑i=0nAilk2(xi)。
因为
∫
a
b
ρ
(
x
)
l
k
2
(
x
)
d
x
>
0
∫_a^bρ(x)l_k^2(x)\mathrm dx>0
∫abρ(x)lk2(x)dx>0,
∑
i
=
0
n
A
i
l
k
2
(
x
i
)
=
A
k
∑_{i=0}^nA_il_k^2(x_i)=A_k
∑i=0nAilk2(xi)=Ak,所以高斯求积系数Ak都为正,所以高斯型求积公式都是数值稳定的。
基于Taylor公式的数值微分公式
向前泰勒公式:
f
(
x
+
h
)
=
f
(
x
)
+
h
f
′
(
x
)
+
h
2
2
!
f
′
′
(
x
)
+
h
3
3
!
f
′
′
′
(
ξ
)
f(x+h)=f(x)+hf'(x)+\frac{h^2}{2!}f''(x)+\frac{h^3}{3!}f'''(ξ)
f(x+h)=f(x)+hf′(x)+2!h2f′′(x)+3!h3f′′′(ξ)
向后泰勒公式:
f
(
x
−
h
)
=
f
(
x
)
−
h
f
′
(
x
)
+
h
2
2
!
f
′
′
(
x
)
−
h
3
3
!
f
′
′
′
(
ξ
)
f(x-h)=f(x)-hf'(x)+\frac{h^2}{2!}f''(x)-\frac{h^3}{3!}f'''(ξ)
f(x−h)=f(x)−hf′(x)+2!h2f′′(x)−3!h3f′′′(ξ)
一阶向前差商公式:
f
′
(
x
)
≈
f
(
x
+
h
)
−
f
(
x
)
h
−
h
2
f
′
′
(
ξ
)
f'(x)≈\frac{f(x+h)-f(x)}h-\frac h2f''(ξ)
f′(x)≈hf(x+h)−f(x)−2hf′′(ξ)
一阶向后差商公式:
f
′
(
x
)
≈
f
(
x
)
−
f
(
x
−
h
)
h
+
h
2
f
′
′
(
ξ
)
f'(x)≈\frac{f(x)-f(x-h)}h+\frac h2f''(ξ)
f′(x)≈hf(x)−f(x−h)+2hf′′(ξ)
中心差商公式:
f
(
x
+
h
)
−
f
(
x
−
h
)
2
h
=
f
′
(
x
)
+
h
2
3
!
f
′
′
′
(
ξ
)
+
f
′
′
′
(
η
)
2
\frac{f(x+h)-f(x-h)}{2h}=f'(x)+\frac{h^2}{3!}\frac{f'''(ξ)+f'''(η)}2
2hf(x+h)−f(x−h)=f′(x)+3!h22f′′′(ξ)+f′′′(η)
单个方程的数值解
二分法
二分法可用于快速缩小方程f(x)=0的有根区间。
首先,设置一个有根区间。由介值定理可知,当f(a0)f(b0)<0时,[a0,b0]为f(x)=0的有根区间。
在第n次搜索时,若
f
(
a
n
−
1
)
f
(
a
n
−
1
+
b
n
−
1
2
)
<
0
f(a_{n-1})f(\frac{a_{n-1}+b_{n-1}}2)<0
f(an−1)f(2an−1+bn−1)<0,则将搜索范围收缩至
[
a
n
−
1
,
a
n
−
1
+
b
n
−
1
2
]
[a_{n-1},\frac{a_{n-1}+b_{n-1}}2]
[an−1,2an−1+bn−1],反之则收缩至
[
a
n
−
1
+
b
n
−
1
2
,
b
n
−
1
]
[\frac{a_{n-1}+b_{n-1}}2,b_{n-1}]
[2an−1+bn−1,bn−1]。
新的搜索区间则成为[an,bn],进行n次搜索后,二分法的误差估计式为
∣
x
n
−
x
∗
∣
≤
b
0
−
a
0
2
n
|x_n-x^*|≤\frac{b_0-a_0}{2^n}
∣xn−x∗∣≤2nb0−a0。
若要求|xn-x*|≤ε,则搜索次数满足
n
≥
log
2
b
0
−
a
0
ε
n≥\log_2\frac{b_0-a_0}ε
n≥log2εb0−a0,但至少要搜索
f
l
o
o
r
(
log
2
b
0
−
a
0
ε
)
\mathrm{floor}(\log_2\frac{b_0-a_0}ε)
floor(log2εb0−a0)次可使误差≤ε。但答案一般是(n+1)。
美中不足的是,二分法不能求偶数重根或复根。
不动点迭代法
不动点迭代法的构造方式是将f(x)改写成x=φ(x)的形式,从而得到迭代公式xk+1=φ(x),φ(x)也被称为迭代函数。
在迭代公式中,若
lim
k
→
∞
x
k
=
x
∗
\lim\limits_{k→∞}x_k=x^*
k→∞limxk=x∗,则迭代公式收敛,且原方程的根x*被称为φ(x)的不动点。
收敛性基本定理:若上述迭代公式中的迭代函数φ∈C[a,b]满足以下性质:
映内性:∀x∈[a,b],φ(x)∈[a,b]
压缩性:∃压缩系数L∈(0,1),使∀xm,xn∈[a,b],|φ(xm)-φ(xn)|≤L|xm-xn|
则迭代公式具有以下性质:
(1)迭代函数φ(x)在[a,b]上存在唯一的不动点x*
(2)在[a,b]上无论选取什么初值,迭代公式最后都能收敛到x*
(3)迭代公式的误差估计式为
∣
x
k
−
x
∗
∣
≤
L
1
−
L
∣
x
k
−
x
k
−
1
∣
|x_k-x^*|≤\frac L{1-L}|x_k-x_{k-1}|
∣xk−x∗∣≤1−LL∣xk−xk−1∣或
∣
x
k
−
x
∗
∣
≤
L
k
1
−
L
∣
x
1
−
x
0
∣
|x_k-x^*|≤\frac {L^k}{1-L}|x_1-x_0|
∣xk−x∗∣≤1−LLk∣x1−x0∣。
局部收敛性定理:设x*为φ的不动点,φ’(x)在x*的某个邻域内存在、连续且满足|φ’(x*)|<1,则迭代公式局部收敛。
收敛阶的判定
若函数f(x)满足f(x*)=f’(x*)=f’’(x*)=⋯=f(n-1)(x*)=0,而f(n)(x*)≠0,则称x*是f(x)的n重根、n阶零点或f(x)在x*处n阶收敛。
若∃C≠0,p≥1,使
lim
k
→
+
∞
∣
x
k
+
1
−
x
∗
(
x
k
−
x
∗
)
p
∣
=
C
\lim\limits_{k→+∞}|\frac{x_{k+1}-x^*}{(x_k-x^*)^p}|=C
k→+∞lim∣(xk−x∗)pxk+1−x∗∣=C,则定义该迭代公式p阶收敛,C为渐进误差常数。
在上式中,若p=1,则C<1时迭代公式才收敛。
对于迭代公式xk+1=φ(xk),若φ’(x*)≠0则其一阶/线性收敛;若φ’(x*)=0、φ"(x*)≠0则其二阶收敛;若φ’(x*)=φ"(x*)=0、φ"’(x*)≠0则其三阶收敛⋯
由泰勒展开又可得f(x*)的误差估计式:x→x*时,
f
(
x
)
≈
f
(
n
)
(
ξ
)
n
!
(
x
−
x
∗
)
n
=
g
(
x
)
(
x
−
x
∗
)
n
f(x)≈\frac{f^{(n)}(ξ)}{n!}(x-x^*)^n=g(x)(x-x^*)^n
f(x)≈n!f(n)(ξ)(x−x∗)n=g(x)(x−x∗)n。
其中,ξ在x与x*之间,g(x*)≠0。
Newton迭代法
牛顿迭代法的迭代公式:
x
k
+
1
=
x
k
−
f
(
x
k
)
f
′
(
x
k
)
,
x
∈
N
x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)},x∈N
xk+1=xk−f′(xk)f(xk),x∈N。
牛顿迭代公式在单根附近至少能保证2阶局部收敛,对于重根的收敛性则没那么好。
在重根附近,如果知道根的重数n,则可将牛顿迭代法改进为
x
k
+
1
=
x
k
−
n
⋅
f
(
x
k
)
f
′
(
x
k
)
,
x
∈
N
x_{k+1}=x_k-n·\frac{f(x_k)}{f'(x_k)},x∈N
xk+1=xk−n⋅f′(xk)f(xk),x∈N,使其保证2阶局部收敛。
但在实际应用中,往往不知道重根的重数。