函数的最佳平方逼近与数据的最小二乘拟合
预备知识
- 赋范线性空间: ( V , ∥ ⋅ ∥ ) (V,\parallel \cdot \parallel) (V,∥⋅∥)。 V V V是线性空间, ∥ ⋅ ∥ \parallel \cdot \parallel ∥⋅∥是 V V V上的范数
- 定义区间 [ a , b ] [a,b] [a,b]上连续函数的集合 C [ a , b ] C[a,b] C[a,b]上的范数
∥ f ∥ 2 = [ ∫ a b f 2 ( x ) d x ] 1 2 \parallel f \parallel_2 = \bigg[ \int_a^bf^2(x) dx \bigg]^{\frac12} ∥f∥2=[∫abf2(x)dx]21
-
两个函数的距离: d ( f , g ) = ∥ f − g ∥ d(f,g) = \parallel f-g \parallel d(f,g)=∥f−g∥
-
最佳逼近:
f ( x ) ∈ C [ a , b ] , ϕ ∈ C [ a , b ] f(x) \in C[a,b],\phi \in C[a,b] f(x)∈C[a,b],ϕ∈C[a,b]是简单子集合。找 φ ∗ ( x ) ∈ ϕ \varphi^*(x) \in \phi φ∗(x)∈ϕ,使误差
∥ f ( x ) − φ ∗ ( x ) ∥ = min φ ∈ ϕ ∥ f ( x ) − φ ( x ) ∥ \parallel f(x) - \varphi^*(x) \parallel = \min_{\varphi \in \phi} \parallel f(x) - \varphi(x) \parallel ∥f(x)−φ∗(x)∥=φ∈ϕmin∥f(x)−φ(x)∥
称 φ ∗ ( x ) \varphi^*(x) φ∗(x)为 ϕ \phi ϕ中的对 f ( x ) f(x) f(x)的最佳逼近函数 -
内积空间:设 V V V是线性空间,定义在 V V V上的"二元"实值函数 ( ⋅ , ⋅ ) (\cdot,\cdot) (⋅,⋅),如果满足
- 正定性: ( f , f ) ≥ 0 , & ( f , f ) = 0 ⟺ f = 0 (f,f) \geq 0,\ \& (f,f)=0 \iff f=0 (f,f)≥0, &(f,f)=0⟺f=0
- 对称性: ( f , g ) = ( g , f ) (f,g) = (g,f) (f,g)=(g,f)
- 线性性: ( r f , g ) = r ( f , g ) , ( f + g , h ) = ( f , h ) + ( g , h ) (rf,g) = r(f,g),(f+g,h) = (f,h)+(g,h) (rf,g)=r(f,g),(f+g,h)=(f,h)+(g,h)
称 ( ⋅ , ⋅ ) (\cdot,\cdot) (⋅,⋅)是线性空间 V V V上的内积,称线性空间 V V V为内积空间
-
权函数
如果函数 ρ ( x ) \rho(x) ρ(x)满足
- ρ ( x ) ≥ 0 , x ∈ [ a , b ] \rho(x) \geq 0,x \in [a,b] ρ(x)≥0,x∈[a,b]
- ∫ a b ρ ( x ) d x > 0 \int_a^b \rho(x) dx > 0 ∫abρ(x)dx>0
- ∫ a b x k ρ ( x ) d x \int_a^b x^k\rho(x) dx ∫abxkρ(x)dx存在
则称 ρ ( x ) \rho(x) ρ(x)是 [ a , b ] [a,b] [a,b]上的权函数,称 ( f , g ) p = ∫ a b ρ ( x ) f ( x ) g ( x ) d x (f,g)_p=\int_a^b \rho(x)f(x)g(x)dx (f,g)p=∫abρ(x)f(x)g(x)dx为加权内积
-
Gram \text{Gram} Gram矩阵
设 V V V是内积空间, φ 0 , φ 1 , ⋯ , φ n ∈ V \varphi_0,\varphi_1,\cdots,\varphi_n \in V φ0,φ1,⋯,φn∈V,称矩阵
G = [ ( φ 0 , φ 0 ) ( φ 0 , φ 1 ) ⋯ ( φ 0 , φ n ) ( φ 1 , φ 0 ) ( φ 1 , φ 1 ) ⋯ ( φ 1 , φ n ) ⋯ ⋯ ⋯ ⋯ ( φ n , φ 0 ) ( φ n , φ 1 ) ⋯ ( φ n , φ n ) ] G = \begin{bmatrix} (\varphi_0,\varphi_0) & (\varphi_0,\varphi_1) &\cdots &(\varphi_0,\varphi_n) \\ (\varphi_1,\varphi_0) & (\varphi_1,\varphi_1) &\cdots &(\varphi_1,\varphi_n) \\ \cdots & \cdots &\cdots &\cdots \\ (\varphi_n,\varphi_0) & (\varphi_n,\varphi_1) &\cdots &(\varphi_n,\varphi_n) \\ \end{bmatrix} G= (φ0,φ0)(φ1,φ0)⋯(φn,φ0)(φ0,φ1)(φ1,φ1)⋯(φn,φ1)⋯⋯⋯⋯(φ0,φn)(φ1,φn)⋯(φn,φn)
为关于元素 φ 0 , φ 1 , ⋯ , φ n \varphi_0,\varphi_1,\cdots,\varphi_n φ0,φ1,⋯,φn的 Gram \text{Gram} Gram矩阵,该矩阵对称正定 -
设 V V V是内积空间,则 φ 0 , φ 1 , ⋯ , φ n ∈ V \varphi_0,\varphi_1,\cdots,\varphi_n \in V φ0,φ1,⋯,φn∈V线性无关 ⟺ \iff ⟺相应的 Gram \text{Gram} Gram矩阵 G G G非奇异
正交多项式系
设 V V V是内积空间, f , g ∈ V f,g \in V f,g∈V,如果内积 ( f , g ) = 0 (f,g)=0 (f,g)=0,则称 f , g f,g f,g关于此内积正交。若进一步 ( f , f ) = ( g , g ) = 1 (f,f) = (g,g) = 1 (f,f)=(g,g)=1,则称 f , g f,g f,g单位正交。
零空间和 V V V中任意元素正交
-
正交系:如果线性空间 V V V中的非零元素系两两正交,则称该非零元素系为正交系
-
正交多项式系:如果正交系中的 { f i } i = 0 ∞ \{f_i\}_{i=0}^\infty {fi}i=0∞为多项式,则为正交多项式系
-
正交多项式系的性质
-
正交多项式系为线性无关系
-
正交多项式和次数低的任意多项式是正交的
-
n n n次正交多项式在正交区间内有 n n n个互异的单实零点
-
正交多项式系相邻三项之间有关系
f k + 1 ( x ) = ( a k x + b k ) f k ( x ) + c k − 1 f k − 1 ( x ) f_{k+1}(x) = (a_kx+b_k)f_k(x) + c_{k-1}f_{k-1}(x) fk+1(x)=(akx+bk)fk(x)+ck−1fk−1(x)
-
-
一般正交多项式系的构造
g 0 = 1 , g 1 = x − a , g 2 = x 2 − b x − c g_0=1,g_1 = x-a,g_2 = x^2-bx-c g0=1,g1=x−a,g2=x2−bx−c
常用的正交多项式
多项式 | 区域 | 权函数 | 定义 |
---|---|---|---|
Legendre \text{Legendre} Legendre | [ − 1 , 1 ] [-1,1] [−1,1] | ρ ( x ) = 1 \rho(x) = 1 ρ(x)=1 | { p 0 ( x ) = 1 p n ( x ) = 1 2 n n ! d n d x n ( x 2 − 1 ) n \begin{cases}p_0(x) = 1 \\p_n(x) = \cfrac{1}{2^n n!}\cfrac{d^n}{dx^n}(x^2-1)^n\end{cases} ⎩ ⎨ ⎧p0(x)=1pn(x)=2nn!1dxndn(x2−1)n |
Chebyshev \text{Chebyshev} Chebyshev | [ − 1 , 1 ] [-1,1] [−1,1] | ρ ( x ) = 1 1 − x 2 \rho(x) = \cfrac{1}{\sqrt{1-x^2}} ρ(x)=1−x21 | { T n ( x ) = cos ( n arccos x ) T n + 1 = 2 x T n − T n − 1 \begin{cases}T_n(x) = \cos (n \arccos x) \\ T_{n+1} = 2xT_n-T_{n-1} \end{cases} {Tn(x)=cos(narccosx)Tn+1=2xTn−Tn−1 |
Laguerre \text{Laguerre} Laguerre | [ 0 , + ∞ ) [0,+\infty) [0,+∞) | ρ = e − x \rho=e^{-x} ρ=e−x | L n ( x ) = e x d n d x n ( x n e − x ) L_n(x) = e^x \frac{d^n}{dx^n}(x^ne^{-x}) Ln(x)=exdxndn(xne−x) |
Hermite \text{Hermite} Hermite | ( − ∞ , + ∞ ) (-\infty,+\infty) (−∞,+∞) | ρ = e − x 2 \rho=e^{-x^2} ρ=e−x2 | H n ( x ) = ( − 1 ) n e x 2 d n d x n ( e − x 2 ) H_n(x) = (-1)^ne^{x^2}\frac{d^n}{dx^n}(e^{-x^2}) Hn(x)=(−1)nex2dxndn(e−x2) |
连续函数的最佳平方逼近
f
(
x
)
∈
C
[
a
,
b
]
,
ϕ
∈
C
[
a
,
b
]
f(x) \in C[a,b],\phi \in C[a,b]
f(x)∈C[a,b],ϕ∈C[a,b]是简单子集合。找
φ
∗
(
x
)
∈
ϕ
\varphi^*(x) \in \phi
φ∗(x)∈ϕ,使误差
∥
f
(
x
)
−
φ
∗
(
x
)
∥
2
=
min
φ
∈
ϕ
∥
f
(
x
)
−
φ
(
x
)
∥
2
\parallel f(x) - \varphi^*(x) \parallel_2 = \min_{\varphi \in \phi} \parallel f(x) - \varphi(x) \parallel_2
∥f(x)−φ∗(x)∥2=φ∈ϕmin∥f(x)−φ(x)∥2
称
φ
∗
(
x
)
\varphi^*(x)
φ∗(x)为
ϕ
\phi
ϕ中的对
f
(
x
)
f(x)
f(x)的最佳平方逼近函数
- 解法
已知基或正交基(使用正交基可简化运算)
φ
i
\varphi_i
φi,求解以下方程组求出
c
∗
c^*
c∗,则
φ
∗
(
x
)
=
∑
i
=
0
n
c
i
∗
φ
i
(
x
)
\varphi^*(x) = \sum_{i=0}^n c_i^*\varphi_i(x)
φ∗(x)=∑i=0nci∗φi(x)
G
c
∗
=
[
(
φ
0
,
φ
0
)
(
φ
0
,
φ
1
)
⋯
(
φ
0
,
φ
n
)
(
φ
1
,
φ
0
)
(
φ
1
,
φ
1
)
⋯
(
φ
1
,
φ
n
)
⋯
⋯
⋯
⋯
(
φ
n
,
φ
0
)
(
φ
n
,
φ
1
)
⋯
(
φ
n
,
φ
n
)
]
[
c
0
∗
c
1
∗
⋮
c
n
∗
]
=
[
(
f
,
φ
0
)
(
f
,
φ
1
)
⋮
(
f
,
φ
n
)
]
G \pmb c^*= \begin{bmatrix} (\varphi_0,\varphi_0) & (\varphi_0,\varphi_1) &\cdots &(\varphi_0,\varphi_n) \\ (\varphi_1,\varphi_0) & (\varphi_1,\varphi_1) &\cdots &(\varphi_1,\varphi_n) \\ \cdots & \cdots &\cdots &\cdots \\ (\varphi_n,\varphi_0) & (\varphi_n,\varphi_1) &\cdots &(\varphi_n,\varphi_n) \\ \end{bmatrix} \begin{bmatrix} c_0^* \\c_1^*\\ \vdots \\ c_n^* \end{bmatrix}= \begin{bmatrix} (f,\varphi_0) \\ (f,\varphi_1) \\ \vdots \\ (f,\varphi_n) \\ \end{bmatrix}
Gc∗=
(φ0,φ0)(φ1,φ0)⋯(φn,φ0)(φ0,φ1)(φ1,φ1)⋯(φn,φ1)⋯⋯⋯⋯(φ0,φn)(φ1,φn)⋯(φn,φn)
c0∗c1∗⋮cn∗
=
(f,φ0)(f,φ1)⋮(f,φn)
-
最佳平方逼近误差
δ 2 = ∥ f ( x ) − φ ∗ ( x ) ∥ 2 2 = ∫ a b ρ ( x ) ( f − φ ∗ ) 2 d x = ∥ f ∥ 2 2 − ∑ k = 0 n c k ∗ ( f , φ k ) \delta^2 = \parallel f(x) - \varphi^*(x) \parallel_2^2 = \int_a^b \rho(x)(f-\varphi^*)^2dx = \parallel f\parallel_2^2 - \sum_{k=0}^n c^*_k (f,\varphi_k) δ2=∥f(x)−φ∗(x)∥22=∫abρ(x)(f−φ∗)2dx=∥f∥22−k=0∑nck∗(f,φk)
离散数据的最小二乘曲线拟合
-
插值函数曲线经过已知点,拟合曲线未必经过已知点
-
数据量很大时,插值函数中的参数个数很多,拟合函数中的参数个数很少
-
连续函数的 Gram \text{Gram} Gram矩阵对称正定,离散函数的 Gram \text{Gram} Gram矩阵对称半正定
和连续函数的最佳平方逼近类似,只需要将连续函数内积的定义改为离散内积的定义,即
- 离散内积: ( f , g ) = ∑ i = 0 m w i f ( x i ) g ( x i ) (f,g) =\sum_{i=0}^m w_i f(x_i)g(x_i) (f,g)=∑i=0mwif(xi)g(xi)
- 连续函数内积: ( f , g ) = ∫ a b w ( x ) f ( x ) g ( x ) d x (f,g) = \int_a^b w(x)f(x)g(x)dx (f,g)=∫abw(x)f(x)g(x)dx
数值积分与数值微分
数值积分的基本概念
数值求积公式的代数精度
I = ∫ a b f ( x ) d x ≈ ∑ i = 0 n A i f ( x i ) = I n I = \int_a^b f(x)dx \approx \sum_{i=0}^n A_if(x_i) = I_n I=∫abf(x)dx≈i=0∑nAif(xi)=In
E n ( f ) = ∫ a b f ( x ) d x − ∑ i = 0 n A i f ( x i ) (截断误差) E_n(f) = \int_a^b f(x) dx - \sum_{i=0}^n A_if(x_i) \tag{截断误差} En(f)=∫abf(x)dx−i=0∑nAif(xi)(截断误差)
数值求积公式 m m m次代数精度: ∀ f ( x ) ∈ P m → E n ( f ) = 0 , f ( x ) = x m + 1 → E n ( f ) ≠ 0 \forall f(x) \in P^m \rightarrow E_n(f)=0, f(x) = x^{m+1} \rightarrow E_n(f) \neq 0 ∀f(x)∈Pm→En(f)=0,f(x)=xm+1→En(f)=0。
等价于 ∀ f ( x ) = x i ( i < = m ) → E n ( f ) = 0 , f ( x ) = x m + 1 → E n ( f ) ≠ 0 \forall f(x) = x^i(i<=m) \rightarrow E_n(f) = 0,f(x) = x^{m+1}\rightarrow E_n(f) \neq 0 ∀f(x)=xi(i<=m)→En(f)=0,f(x)=xm+1→En(f)=0
在区间 [ a , b ] [a,b] [a,b]上给定 [ n + 1 ] [n+1] [n+1]个求积节点 { x i } i = 0 n \{x_i\}_{i=0}^n {xi}i=0n,则存在唯一的一组求积系数 { A i } i = 0 n \{A_i\}_{i=0}^n {Ai}i=0n使得求积公式至少具有 n n n次代数精度
求积公式的收敛性与稳定性
-
收敛性
∀ f ( x ) ∈ C [ a , b ] , lim n → ∞ h → 0 ∑ k = 0 n A k f ( x k ) = ∫ a b f ( x ) d x \forall f(x) \in C[a,b], \lim_{\begin{aligned} &n\rightarrow \infty \\ &h\rightarrow 0\end{aligned}} \sum_{k=0}^n A_k f(x_k) = \int_a^bf(x)dx ∀f(x)∈C[a,b],n→∞h→0limk=0∑nAkf(xk)=∫abf(x)dx -
稳定性
设计算 f ( x k ) f(x_k) f(xk)有舍入误差 e k = f ∗ ( x k ) − f ( x ) e_k = f^*(x_k) - f(x) ek=f∗(xk)−f(x)。有 e ( I n ) = ∑ k = 0 n A i e k ≤ ε ∑ k = 0 n A i e(I_n) = \sum_{k=0}^n A_ie_k \leq \varepsilon \sum_{k=0}^n A_i e(In)=∑k=0nAiek≤ε∑k=0nAi。如果 A i ≥ 0 A_i \geq 0 Ai≥0,具有代数精度为 0 0 0的求积公式系数之和为 b − a b-a b−a,则 ∣ e ( I n ) ∣ ≤ ε ( b − a ) \mid e(I_n) \mid \leq \varepsilon (b-a) ∣e(In)∣≤ε(b−a)。
因此,如果系数全非负,系数之和等于区间长度,求积公式稳定
插值型求积公式
插值型求积公式
使用插值函数近似原函数,进行求积运算
∫
a
b
f
(
x
)
d
x
≈
∫
a
b
p
n
(
x
)
d
x
\int_a^b f(x) dx \approx \int_a^b p_n(x)dx
∫abf(x)dx≈∫abpn(x)dx
E n ( f ) = ∫ a b f ( n + 1 ) ( ξ ) ( n + 1 ) ! ω n + 1 ( x ) d x E_n(f) = \int_a^b \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x) dx En(f)=∫ab(n+1)!f(n+1)(ξ)ωn+1(x)dx
- 具有 n + 1 n+1 n+1个求积节点的数值求积公式是插值型求积公式的充要条件为该求积公式至少具有 n n n次代数精度
Newton-Cotes求积公式
当求积节点等距分布时,有
∫
a
b
f
(
x
)
d
x
≈
∑
i
=
0
n
c
i
(
n
)
f
(
a
+
i
h
)
\int_a^b f(x) dx \approx \sum_{i=0}^n c_i^{(n)} f(a+ih)
∫abf(x)dx≈i=0∑nci(n)f(a+ih)
- 节点数 n n n为奇数, Newton-cotes \text{Newton-cotes} Newton-cotes代数精度至少 n n n次;为偶数时,代数精度至少 n − 1 n-1 n−1次
求积公式 | ∫ a b f ( x ) d x ≈ \int_a^b f(x)dx \approx ∫abf(x)dx≈ | 截断误差 | 代数精度 |
---|---|---|---|
一点中矩形求积公式 | $ (b-a) f(\frac{a+b}2)$ | ( b − a ) 3 24 f ′ ′ ( η ) \frac{(b-a)^3}{24}f^{\prime \prime}(\eta) 24(b−a)3f′′(η) | 1 |
两点梯形求积公式 | b − a 2 [ f ( a ) + f ( b ) ] \frac{b-a}2 [f(a) + f(b)] 2b−a[f(a)+f(b)] | − ( b − a ) 3 12 f ′ ′ ( η ) -\frac{(b-a)^3}{12}f^{\prime \prime}(\eta) −12(b−a)3f′′(η) | 1 |
三点抛物型/ S i m p s o n Simpson Simpson | b − a 6 [ f ( a ) + 4 f ( a + b 2 ) + f ( b ) ] \frac{b-a}6 [f(a) + 4f(\frac{a+b}2)+f(b)] 6b−a[f(a)+4f(2a+b)+f(b)] | − ( b − a ) 5 2880 f ( 4 ) ( η ) -\frac{(b-a)^5}{2880}f^{(4)}(\eta) −2880(b−a)5f(4)(η) | 3 |
其中一点中矩形求积公式的截断误差使用泰勒展开求解,其余使用插值求积公式的截断误差公式代入求解
∫
a
b
f
(
x
)
d
x
=
(
b
−
a
)
f
(
a
+
b
2
)
⟹
f
(
x
)
=
f
(
a
+
b
2
)
+
(
x
−
a
+
b
2
)
f
′
(
a
+
b
2
)
+
(
x
−
a
+
b
2
)
2
f
′
′
(
ξ
)
2
!
⟺
∫
a
b
f
(
x
)
d
x
=
(
b
−
a
)
f
(
a
+
b
2
)
+
∫
a
b
f
′
′
(
ξ
)
2
!
(
x
−
a
+
b
2
)
2
d
x
Since
f
′
′
(
x
)
is continuous on
[
a
,
b
]
,
(
x
−
a
+
b
2
)
2
unchanged sign
By the median integral theorem
,
∃
η
∈
[
a
,
b
]
⟹
∫
a
b
f
(
x
)
d
x
=
(
b
−
a
)
f
(
a
+
b
2
)
+
f
′
′
(
η
)
2
!
∫
a
b
(
x
−
a
+
b
2
)
2
d
x
⟹
E
(
f
)
=
f
′
′
(
η
)
2
!
∫
a
b
(
x
−
a
+
b
2
)
2
d
x
=
(
b
−
a
)
3
24
f
′
′
(
η
)
\begin{aligned} &\int_a^b f(x)dx = (b-a)f(\frac{a+b}2) \\ \Longrightarrow&f(x) = f(\frac{a+b}{2})+(x-\frac{a+b}{2})f^\prime(\frac{a+b}2) + (x-\frac{a+b}{2})^2\frac{f^{\prime \prime}(\xi)}{2!} \\ \iff& \int_a^b f(x)dx = (b-a)f(\frac{a+b}2) + \int_a^b \frac{f^{\prime \prime}(\xi)}{2!}(x-\frac{a+b}{2})^2dx \\ &\text{Since } f^{\prime \prime}(x) \text{ is continuous on } [a,b],(x-\frac{a+b}2)^2 \text{ unchanged sign }\\ &\text{By the median integral theorem},\exist\eta \in [a,b]\\ \Longrightarrow&\int_a^bf(x)dx = (b-a)f(\frac{a+b}2) + \frac{f^{\prime \prime}(\eta)}{2!}\int_a^b (x-\frac{a+b}{2})^2dx\\ \Longrightarrow& E(f) = \frac{f^{\prime \prime}(\eta)}{2!}\int_a^b (x-\frac{a+b}{2})^2dx = \frac{(b-a)^3}{24} f^{\prime \prime}(\eta) \end{aligned}
⟹⟺⟹⟹∫abf(x)dx=(b−a)f(2a+b)f(x)=f(2a+b)+(x−2a+b)f′(2a+b)+(x−2a+b)22!f′′(ξ)∫abf(x)dx=(b−a)f(2a+b)+∫ab2!f′′(ξ)(x−2a+b)2dxSince f′′(x) is continuous on [a,b],(x−2a+b)2 unchanged sign By the median integral theorem,∃η∈[a,b]∫abf(x)dx=(b−a)f(2a+b)+2!f′′(η)∫ab(x−2a+b)2dxE(f)=2!f′′(η)∫ab(x−2a+b)2dx=24(b−a)3f′′(η)
E 1 ( f ) = ∫ a b f ′ ′ ( ξ ) 2 ! ( x − a ) ( x − b ) d x = f ′ ′ ( η ) 2 ! ∫ a b ( x − a ) ( x − b ) d x = − ( b − a ) 3 12 f ′ ′ ( η ) E_1(f) = \int_a^b \frac{f^{\prime \prime }(\xi)}{2!} (x-a)(x-b)dx = \frac{f^{\prime \prime }(\eta)}{2!} \int_a^b(x-a)(x-b)dx = -\frac{(b-a)^3}{12}f^{\prime \prime}(\eta) E1(f)=∫ab2!f′′(ξ)(x−a)(x−b)dx=2!f′′(η)∫ab(x−a)(x−b)dx=−12(b−a)3f′′(η)
复化求积算法
将
[
a
,
b
]
n
[a,b]n
[a,b]n等份,子区间
[
x
i
,
x
i
+
1
]
(
i
=
0
,
1
,
⋯
,
n
−
1
)
[x_i,x_{i+1}](i=0,1,\cdots,n-1)
[xi,xi+1](i=0,1,⋯,n−1),步长
h
=
b
−
a
n
h=\frac{b-a}{n}
h=nb−a
I
(
f
)
=
∫
a
b
f
(
x
)
d
x
=
∑
i
=
0
n
−
1
∫
x
i
x
i
+
1
f
(
x
)
d
x
I(f) = \int_a^b f(x)dx = \sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}}f(x)dx
I(f)=∫abf(x)dx=i=0∑n−1∫xixi+1f(x)dx
复化梯形公式
T n ( f ) = ∫ a b f ( x ) d x = ∑ i = 0 n − 1 ∫ x i x i + 1 f ( x ) d x = ∑ i = 0 n − 1 h 2 [ f ( x i ) + f ( x i + 1 ) ] − h 3 12 ∑ i = 0 n − 1 f ′ ′ ( η i ) = h 2 [ f ( a ) + 2 ∑ i = 0 n − 1 + f ( b ) ] − ( b − a ) h 2 12 [ 1 n ∑ i = 0 n − 1 f ′ ′ ( η i ) ] \begin{aligned} T_n(f) &= \int_a^bf(x)dx = \sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}} f(x)dx\\ &=\sum_{i=0}^{n-1}\frac h2 [f(x_i) + f(x_{i+1})] - \frac{h^3}{12}\sum_{i=0}^{n-1}f^{\prime \prime}(\eta_i)\\ &=\frac{h}2[f(a)+2\sum_{i=0}^{n-1}+f(b)] - (b-a) \frac{h^2}{12} \Big[\frac1n\sum_{i=0}^{n-1}f^{\prime \prime }(\eta_i) \Big] \end{aligned} Tn(f)=∫abf(x)dx=i=0∑n−1∫xixi+1f(x)dx=i=0∑n−12h[f(xi)+f(xi+1)]−12h3i=0∑n−1f′′(ηi)=2h[f(a)+2i=0∑n−1+f(b)]−(b−a)12h2[n1i=0∑n−1f′′(ηi)]
当
f
′
′
(
x
)
f^{\prime \prime }(x)
f′′(x)在区间
[
a
,
b
]
[a,b]
[a,b]连续时,利用连续函数介值定理值,
∃
η
∈
(
a
,
b
)
\exists \eta \in (a,b)
∃η∈(a,b)
1
n
∑
i
=
0
n
−
1
f
′
′
(
η
i
)
=
f
′
′
(
η
)
\frac 1n \sum_{i=0}^{n-1} f^{\prime \prime }(\eta_i) = f^{\prime \prime }(\eta)
n1i=0∑n−1f′′(ηi)=f′′(η)
T n ( f ) = h 2 [ f ( a ) + 2 ∑ i = 1 n − 1 f ( x i ) + f ( b ) ] T_n(f) = \frac h2 \bigg[f(a) +2\sum_{i=1}^{n-1} f(x_i) + f(b) \bigg] Tn(f)=2h[f(a)+2i=1∑n−1f(xi)+f(b)]
E T ( f ) = − ( b − a ) h 2 12 f ′ ′ ( η ) E_T(f) = -\frac{(b-a)h^2}{12} f^{\prime \prime} (\eta) ET(f)=−12(b−a)h2f′′(η)
复化Simpson求积公式
在每个子区间
[
x
i
,
x
i
+
1
]
[x_i,x_{i+1}]
[xi,xi+1]插入中点
x
i
+
1
/
2
x_{i+1/2}
xi+1/2
S
n
(
f
)
=
h
6
[
f
(
a
)
+
2
∑
i
=
1
n
−
1
f
(
x
i
)
+
4
∑
i
=
0
n
−
1
f
(
x
i
+
1
/
2
)
+
f
(
b
)
]
S_n(f) = \frac h6 \bigg[ f(a) + 2 \sum_{i=1}^{n-1}f(x_i) + 4 \sum_{i=0}^{n-1} f(x_{i+1/2})+f(b)\bigg]
Sn(f)=6h[f(a)+2i=1∑n−1f(xi)+4i=0∑n−1f(xi+1/2)+f(b)]
E S ( f ) = − b − a 2880 h 4 f ( 4 ) ( η ) E_S(f) = -\frac{b-a}{2880}h^4 f^{(4)}(\eta) ES(f)=−2880b−ah4f(4)(η)
事后误差估计—区间逐次分半算法
I ≈ T 2 n + 1 3 ( T 2 n − T n ) = S n I \approx T_{2n} + \frac13 (T_{2n}-T_n) = S_n I≈T2n+31(T2n−Tn)=Sn
I ≈ S 2 n + 1 15 ( S 2 n − S n ) = C n I \approx S_{2n} + \frac 1{15}(S_{2n}-S_n) = C_n I≈S2n+151(S2n−Sn)=Cn
I ≈ C 2 n + 1 63 ( C 2 n − C n ) = R n I \approx C_{2n} +\frac{1}{63}(C_{2n}-C_n) = R_n I≈C2n+631(C2n−Cn)=Rn
Romberg求积算法
将事后误差估计 T n → S n → C n → R n T_n\rightarrow S_n \rightarrow C_n \rightarrow R_n Tn→Sn→Cn→Rn列成表格就是 Romberg \text{Romberg} Romberg求积公式
外推技巧
Q ∗ ( h ) = Q + ⋯ Q^*(h) = Q + \cdots Q∗(h)=Q+⋯
Q
∗
(
h
2
)
=
Q
+
⋯
Q^*(\frac h2) =Q+ \cdots
Q∗(2h)=Q+⋯
4
Q
∗
(
h
2
)
−
Q
∗
(
h
)
4
−
1
=
Q
+
⋯
\frac{4Q^*(\frac h2)-Q^*(h)}{4-1} = Q+\cdots
4−14Q∗(2h)−Q∗(h)=Q+⋯
Gauss型求积公式
对于插值求积公式 ∫ a b ρ ( x ) f ( x ) d x ≈ ∑ k = 0 n A k f ( x k ) \int_a^b \rho(x) f(x) dx \approx \sum_{k=0}^n A_k f(x_k) ∫abρ(x)f(x)dx≈∑k=0nAkf(xk)的代数精度至少为 n n n次。公式含有 2 n + 2 2n+2 2n+2个待定参数(节点和系数),可以设 2 n + 2 2n+2 2n+2个公式满足左右侧相等,其中 f = 1 , x , ⋯ , x 2 n + 1 f=1,x,\cdots,x^{2n+1} f=1,x,⋯,x2n+1,解出节点和系数。
- 插值型求积公式的代数精度为 2 n + 1 2n+1 2n+1次的求积公式为 Gauss \text{Gauss} Gauss型求积公式
- x k x_k xk为 Gauss \text{Gauss} Gauss点
- A k A_k Ak为 Gauss \text{Gauss} Gauss系数
- Gauss \text{Gauss} Gauss型求积公式是精度最高的求积公式
Gauss求积公式构造
-
插值型求积公式 ∫ a b ρ ( x ) f ( x ) d x ≈ ∑ k = 0 n A k f ( x k ) \int_a^b \rho(x) f(x) dx \approx \sum_{k=0}^n A_k f(x_k) ∫abρ(x)f(x)dx≈∑k=0nAkf(xk)中的节点为 Gauss \text{Gauss} Gauss点的充要条件是:以这些点为零点的 n + 1 n+1 n+1多项式 ω n + 1 ( x ) = ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n ) \omega_{n+1}(x) = (x-x_0)(x-x_1) \cdots (x-x_n) ωn+1(x)=(x−x0)(x−x1)⋯(x−xn)正交于空间 P n P^n Pn
∫ a b ρ ( x ) ω n + 1 ( x ) p n ( x ) d x = 0 , ∀ p n ∈ P n \int_a^b \rho(x) \omega_{n+1}(x) p_n(x)dx = 0,\forall p_n \in P^n ∫abρ(x)ωn+1(x)pn(x)dx=0,∀pn∈Pn -
Gauss \text{Gauss} Gauss型求积公式是插值公式,满足
A k = ∫ a b l k ( x ) d x A_k = \int_a^b l_k(x)dx Ak=∫ablk(x)dx
步骤
-
确定与不超过 n n n次多项式 p ( x ) p(x) p(x)都正交的 n + 1 n+1 n+1次首一多项式 ω n + 1 ( x ) \omega_{n+1}(x) ωn+1(x)
( p , ω n + 1 ) = ∫ a b p ( x ) ω n + 1 ( x ) ρ ( x ) d x = 0 (p,\omega_{n+1}) = \int_a^b p(x) \omega_{n+1}(x) \rho(x) dx=0 (p,ωn+1)=∫abp(x)ωn+1(x)ρ(x)dx=0 -
求 ω n + 1 ( x ) \omega_{n+1}(x) ωn+1(x)的所有零点
-
由代数精度的等价定义列方程组—待定系数法
截断误差
E n ( f ) = f ( 2 n + 2 ) ( η ) ( 2 n + 2 ) ! ∫ a b ρ ( x ) [ ω n + 1 ( x ) ] 2 d x E_n(f) = \frac{f^{(2n+2)}(\eta)}{(2n+2)!} \int_a^b \rho(x) \Big[\omega_{n+1}(x)\Big]^2 dx En(f)=(2n+2)!f(2n+2)(η)∫abρ(x)[ωn+1(x)]2dx
常见的Gauss求积公式
-
Gauss-Legendre \text{Gauss-Legendre} Gauss-Legendre求积公式
∫ − 1 1 f ( x ) d x ≈ ∑ i = 0 n A i f ( x i ) \int_{-1}^1f(x)dx \approx \sum_{i=0}^n A_if(x_i) ∫−11f(x)dx≈i=0∑nAif(xi)-
节点: n + 1 n+1 n+1次 Lengendre \text{Lengendre} Lengendre多项式 P n + 1 ( x ) P_{n+1}(x) Pn+1(x)的零点
-
系数:
- A i = ∫ − 1 1 l i ( x ) d x A_i = \int_{-1}^1 l_i(x)dx Ai=∫−11li(x)dx
- 待定系数法
- A i = 2 ( 1 − x i 2 ) [ p n + 1 ′ ( x i ) ] 2 A_i = \cfrac{2 }{(1-x_i^2)[p^\prime_{n+1}(x_i)]^2} Ai=(1−xi2)[pn+1′(xi)]22
-
截断误差:略
-
一点 Gauss-Legendre \text{Gauss-Legendre} Gauss-Legendre求积公式
∫ − 1 1 f ( x ) d x ≈ 2 f ( 0 ) \int_{-1}^1f(x)dx \approx 2f(0) ∫−11f(x)dx≈2f(0) -
两点 Gauss-Legendre \text{Gauss-Legendre} Gauss-Legendre求积公式
∫ − 1 1 f ( x ) d x ≈ f ( − 3 3 ) + f ( 3 3 ) \int_{-1}^1 f(x)dx \approx f(-\frac{\sqrt3}3) + f(\frac{\sqrt{3}}3) ∫−11f(x)dx≈f(−33)+f(33)
-
-
Gauss—Chebyshev \text{Gauss—Chebyshev} Gauss—Chebyshev求积公式
∫ − 1 1 1 1 − x 2 f ( x ) d x ≈ π n + 1 ∑ i = 0 n f ( cos 2 i + 1 2 n + 2 π ) \int_{-1}^{1} \frac1{\sqrt{1-x^2} }f(x)dx \approx \frac{\pi}{n+1} \sum_{i=0}^n f(\cos \frac{2i+1}{2n+2}\pi) ∫−111−x21f(x)dx≈n+1πi=0∑nf(cos2n+22i+1π) -
Gauss-Laguerre \text{Gauss-Laguerre} Gauss-Laguerre求积公式
∫ 0 + ∞ e − x f ( x ) d x ≈ ∑ k = 0 n A k f ( x k ) \int_0^{+\infty}e^{-x} f(x) dx \approx \sum_{k=0}^n A_k f(x_k) ∫0+∞e−xf(x)dx≈k=0∑nAkf(xk) -
Gauss-Hermite \text{Gauss-Hermite} Gauss-Hermite求积公式
∫ − ∞ + ∞ e − x 2 f ( x ) d x ≈ ∑ k = 0 n A k f ( x k ) \int_{-\infty}^{+\infty}e^{-x^2} f(x)dx \approx \sum_{k=0}^n A_k f(x_k) ∫−∞+∞e−x2f(x)dx≈k=0∑nAkf(xk)
数值微分
用离散数据 y i = f ( x i ) y_i=f(x_i) yi=f(xi)近似表示函数 y y y在一些节点的导数
如果插值区间
[
a
,
b
]
[a,b]
[a,b]上
∣
f
(
n
+
1
)
(
ξ
n
)
∣
\mid f^{(n+1)}(\xi_n) \mid
∣f(n+1)(ξn)∣有上界
M
M
M,则误差
R
n
′
(
x
i
)
=
f
(
n
+
1
)
(
ξ
n
)
(
n
+
1
)
!
ω
n
+
1
′
(
x
i
)
≤
M
(
n
+
1
)
!
(
b
−
a
)
n
R_n^\prime(x_i) = \frac{f^{(n+1)}(\xi_n)}{(n+1)!}\omega_{n+1}^\prime(x_i) \leq \frac{M}{(n+1)!}(b-a)^n
Rn′(xi)=(n+1)!f(n+1)(ξn)ωn+1′(xi)≤(n+1)!M(b−a)n
插值法
f ′ ( x i ) − p n ′ ( x i ) = R n ′ ( x i ) = f ( n + 1 ) ( ξ n ) ( n + 1 ) ! ω n + 1 ′ ( x i ) f^\prime(x_i ) - p_n^\prime(x_i) = R_n^\prime(x_i) = \frac{f^{(n+1)}(\xi_n)}{(n+1)!}\omega_{n+1}^\prime(x_i) f′(xi)−pn′(xi)=Rn′(xi)=(n+1)!f(n+1)(ξn)ωn+1′(xi)
-
一阶两点公式
f ′ ( x 0 ) = 1 h [ f ( x 1 ) − f ( x 0 ) ] − h 2 f ′ ′ ( ξ 1 ) f ′ ( x 1 ) = 1 h [ f ( x 1 ) − f ( x 0 ) ] + h 2 f ′ ′ ( ξ 2 ) \begin{aligned} &f^\prime(x_0) = \frac 1h[f(x_1)-f(x_0)] - \frac{h}{2} f^{\prime \prime}(\xi_1)\\ &f^\prime(x_1) = \frac 1h[f(x_1)-f(x_0)] + \frac{h}{2} f^{\prime \prime}(\xi_2) \end{aligned} f′(x0)=h1[f(x1)−f(x0)]−2hf′′(ξ1)f′(x1)=h1[f(x1)−f(x0)]+2hf′′(ξ2) -
一阶三点公式
f ′ ( x 0 ) = 1 2 h [ − 3 f ( x 0 ) + 4 f ( x 1 ) − f ( x 2 ) ] + h 2 3 f ( 3 ) ( ξ 1 ) f ′ ( x 1 ) = 1 2 h [ f ( x 2 ) − f ( x 0 ) ] − h 2 6 f ( 3 ) ( ξ 2 ) f ′ ( x 2 ) = 1 2 h [ f ( x 0 ) − 4 f ( x 1 ) + 3 f ( x 2 ) ] + h 2 3 f ( 3 ) ( ξ 3 ) \begin{aligned} &f^\prime(x_0) = \frac 1{2h}[-3f(x_0) + 4f(x_1) - f(x_2)] + \frac{h^2}3 f^{(3)}(\xi_1) \\ &f^\prime(x_1)=\frac 1{2h}[f(x_2)-f(x_0)] - \frac{h^2}{6} f^{(3)}(\xi_2) \\ &f^\prime(x_2) = \frac{1}{2h} [f(x_0)-4f(x_1)+3f(x_2)] + \frac {h^2}3f^{(3)}(\xi_3) \end{aligned} f′(x0)=2h1[−3f(x0)+4f(x1)−f(x2)]+3h2f(3)(ξ1)f′(x1)=2h1[f(x2)−f(x0)]−6h2f(3)(ξ2)f′(x2)=2h1[f(x0)−4f(x1)+3f(x2)]+3h2f(3)(ξ3) -
二阶三点公式
f ′ ′ ( x 0 ) = 1 h 2 [ f ( x 0 ) − 2 f ( x 1 ) + f ( x 2 ) ] − h 2 12 f ( 4 ) ( ξ ) \begin{aligned} f^{\prime \prime} (x_0) = \frac1{h^2} [f(x_0)-2f(x_1) + f(x_2)] - \frac{h^2}{12}f^{(4)}(\xi) \end{aligned} f′′(x0)=h21[f(x0)−2f(x1)+f(x2)]−12h2f(4)(ξ)
Taylor展开法
求哪一点的导数,就把右端函数值在哪一点 Taylor \text{Taylor} Taylor展开
数值微分的外推
从收敛性上讲,步长 h h h越小越好。 从稳定性上讲,步长 h h h因出现在分母,故不宜过小。
常微分方程初值问题的数值解法
{ d y d x = f ( x , y ) , a < x ≤ b y ( a ) = y 0 \begin{cases} \cfrac{dy}{dx} = f(x,y) ,a < x \leq b \\ y(a) = y_0 \end{cases} ⎩ ⎨ ⎧dxdy=f(x,y),a<x≤by(a)=y0
f
(
x
,
y
)
f(x,y)
f(x,y)连续,且关于
y
y
y满足
Lipschitz
\text{Lipschitz}
Lipschitz条件,则解存在唯一,且稳定(适定)
∣
f
(
x
,
y
1
)
−
f
(
x
,
y
2
)
∣
≤
L
∣
y
1
−
y
2
∣
(Lipschitz)
\mid f(x,y_1) - f(x, y_2) \mid \leq L \mid y_1 - y_2 \mid \tag{Lipschitz}
∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣(Lipschitz)
目的是求关于常微分方程初值问题的
y
(
x
1
)
,
y
(
x
2
)
,
⋯
,
y
(
x
n
)
y(x_1),y(x_2),\cdots,y(x_n)
y(x1),y(x2),⋯,y(xn)的近似值
y
1
,
y
2
,
⋯
,
y
n
y_1,y_2,\cdots,y_n
y1,y2,⋯,yn
几种简单的单步法
设节点 x n = x 0 + n h x_n=x_0+nh xn=x0+nh
显式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)
隐式Euler公式
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 ( 0 ) = y n + h f ( x n , y n ) 显格式提供初值 y n + 1 ( s + 1 ) = y n + h f ( x n + 1 , y n + 1 ( s ) ) 隐格式反复迭代 \begin{cases} y_{n+1}^{(0)} = y_n +hf(x_n,y_n) \text{显格式提供初值} \\ y_{n+1}^{(s+1)} = y_n +hf(x_{n+1},y_{n+1}^{(s)}) \text{隐格式反复迭代} \end{cases} {yn+1(0)=yn+hf(xn,yn)显格式提供初值yn+1(s+1)=yn+hf(xn+1,yn+1(s))隐格式反复迭代
- 当 h L < 1 hL<1 hL<1, Euler \text{Euler} Euler隐格式关于 s s s的迭代收敛
梯形方法
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)]
- 隐格式求解:与隐式 Euler \text{Euler} Euler类似, Euler \text{Euler} Euler显格式提供初值,梯形隐格式反复迭代
- 当 h L 2 < 1 \frac{hL}2 <1 2hL<1,梯形隐格式关于 s s s的迭代收敛
- Euler \text{Euler} Euler——梯形预校格式:显格式提供初值,隐格式只算一次
单步法的局部截断误差和阶
-
局部截断误差
假设 y n = y ( x n ) y_n=y(x_n) yn=y(xn), R n + 1 = y ( x n + 1 ) − y n + 1 = O ( h p + 1 ) R_{n+1} = y(x_{n+1})-y_{n+1}=O(h^{p+1}) Rn+1=y(xn+1)−yn+1=O(hp+1),则相应的方法为 p p p阶方法
-
整体截断误差
e n + 1 = y ( x n + 1 ) − y n + 1 e_{n+1} = y(x_{n+1})-y_{n+1} en+1=y(xn+1)−yn+1 -
局部—整体截断误差关系
R n + 1 = O ( h p + 1 ) → e n + 1 = O ( h p ) R_{n+1} = O(h^{p+1}) \rightarrow e_{n+1} = O(h^p) Rn+1=O(hp+1)→en+1=O(hp)
高阶单步法Runge-Kutta法
显式
Runge-Kutta
\text{Runge-Kutta}
Runge-Kutta法的一般形式
{
y
n
+
1
=
y
n
+
h
(
c
1
k
1
+
c
2
k
2
+
⋯
+
c
r
k
r
)
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
α
2
h
,
y
n
+
h
β
21
k
1
)
⋯
k
r
=
f
(
x
n
+
α
r
h
,
y
n
+
h
(
β
r
1
k
1
+
⋯
+
β
r
,
r
−
1
k
r
−
1
)
)
\begin{cases} y_{n+1} = y_n + h(c_1k_1+c_2k_2 + \cdots + c_rk_r) \\ k_1 = f(x_n,y_n) \\ k_2 = f(x_n+\alpha_2h,y_n+h \beta_{21}k_1) \\ \cdots \\ k_r = f(x_n+\alpha_rh,y_n+h(\beta_{r1}k_1+\cdots+\beta_{r,r-1}k_{r-1})) \end{cases}
⎩
⎨
⎧yn+1=yn+h(c1k1+c2k2+⋯+crkr)k1=f(xn,yn)k2=f(xn+α2h,yn+hβ21k1)⋯kr=f(xn+αrh,yn+h(βr1k1+⋯+βr,r−1kr−1))
-
二级二阶Runge-kutta的一种形式是Eluer预校格式
-
n级n阶Runge-kutta的稳定区域是
∣ 1 + λ h + 1 2 ! ( λ h ) 2 + 1 3 ! ( λ h ) 3 + ⋯ ∣ ≤ 1 ( n ≤ 4 ) \mid 1+ \lambda h + \frac{1}{2!}(\lambda h)^2 + \frac{1}{3!}(\lambda h)^3+ \cdots \mid \leq 1 \ \ \ \ \ \ \ \ \ \ \ \ (n \leq 4) ∣1+λh+2!1(λh)2+3!1(λh)3+⋯∣≤1 (n≤4)
单步法的收敛性、相容性与稳定性
- 收敛性: h → 0 , y n + 1 → y ( x n + 1 ) h\rightarrow 0,y_{n+1} \rightarrow y(x_{n+1}) h→0,yn+1→y(xn+1)
- 稳定性:考虑 d y d x = λ y \cfrac{dy}{dx}=\lambda y dxdy=λy,如果 y n y_n yn有舍入误差 δ n \delta_n δn,由此引起计算 y n + 1 y_{n+1} yn+1有误差 δ n + 1 \delta_{n+1} δn+1。如果 ∣ δ n + 1 ∣ ≤ ∣ δ n ∣ \mid \delta_{n+1} \mid \leq \mid \delta_n\mid ∣δn+1∣≤∣δn∣,则关于步长绝对稳定
λ h \lambda h λh | |
---|---|
Euler显格式 | [ − 2 , 0 ] [-2,0] [−2,0] |
Euler隐格式 | R R R |
梯形格式 | R R R |
二级二阶Runge-kutta | ∣ 1 + λ h + 1 2 ( λ h ) 2 ∣ < 1 \mid 1+\lambda h + \frac12 (\lambda h)^2 \mid<1 ∣1+λh+21(λh)2∣<1 |
线性多步法
y n + 1 = ∑ i = 0 k − 1 α i y n − i + h ∑ i = − 1 k − 1 β i f n − i y_{n+1} = \sum_{i=0}^{k-1} \alpha_i y_{n-i} + h \sum_{i=-1}^{k-1}\beta_i f_{n-i} yn+1=i=0∑k−1αiyn−i+hi=−1∑k−1βifn−i
- β − 1 = 0 \beta_{-1}=0 β−1=0时为显格式,否则为隐格式
y n + 1 = y n + h 24 ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) (Adams) y_{n+1} = y_n + \frac h{24}(55f_n-59f_{n-1}+37f_{n-2}-9f_{n-3}) \tag{Adams} yn+1=yn+24h(55fn−59fn−1+37fn−2−9fn−3)(Adams)
R n + 1 = 251 720 h 5 y ( 5 ) ( η ) R_{n+1} = \frac{251}{720}h^5y^{(5)}(\eta) Rn+1=720251h5y(5)(η)