数值分析(三)

函数的最佳平方逼近与数据的最小二乘拟合

预备知识

  • 赋范线性空间: ( 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} f2=[abf2(x)dx]21

  • 两个函数的距离: d ( f , g ) = ∥ f − g ∥ d(f,g) = \parallel f-g \parallel d(f,g)=∥fg

  • 最佳逼近:

    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)∥=φϕminf(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)=0f=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,,φnV,称矩阵
    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,,φnV线性无关    ⟺    \iff 相应的 Gram \text{Gram} Gram矩阵 G G G非奇异

正交多项式系

V V V是内积空间, f , g ∈ V f,g \in V f,gV,如果内积 ( 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)+ck1fk1(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=xa,g2=x2bxc

常用的正交多项式

多项式区域权函数定义
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(x21)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)=1x2 1 { 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=2xTnTn1
Laguerre \text{Laguerre} Laguerre [ 0 , + ∞ ) [0,+\infty) [0,+) ρ = e − x \rho=e^{-x} ρ=ex 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(xnex)
Hermite \text{Hermite} Hermite ( − ∞ , + ∞ ) (-\infty,+\infty) (,+) ρ = e − x 2 \rho=e^{-x^2} ρ=ex2 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(ex2)

连续函数的最佳平方逼近

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=φϕminf(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) c0c1cn = (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=∥f22k=0nck(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)dxi=0nAif(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)dxi=0nAif(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)PmEn(f)=0,f(x)=xm+1En(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+1En(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],nh0limk=0nAkf(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 Ai0,具有代数精度为 0 0 0的求积公式系数之和为 b − a b-a ba,则 ∣ e ( I n ) ∣ ≤ ε ( b − a ) \mid e(I_n) \mid \leq \varepsilon (b-a) e(In)∣≤ε(ba)

    因此,如果系数全非负,系数之和等于区间长度,求积公式稳定

插值型求积公式

插值型求积公式

使用插值函数近似原函数,进行求积运算
∫ 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)dxabpn(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)dxi=0nci(n)f(a+ih)

  • 节点数 n n n为奇数, Newton-cotes \text{Newton-cotes} Newton-cotes代数精度至少 n n n次;为偶数时,代数精度至少 n − 1 n-1 n1
求积公式 ∫ 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(ba)3f′′(η)1
两点梯形求积公式 b − a 2 [ f ( a ) + f ( b ) ] \frac{b-a}2 [f(a) + f(b)] 2ba[f(a)+f(b)] − ( b − a ) 3 12 f ′ ′ ( η ) -\frac{(b-a)^3}{12}f^{\prime \prime}(\eta) 12(ba)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)] 6ba[f(a)+4f(2a+b)+f(b)] − ( b − a ) 5 2880 f ( 4 ) ( η ) -\frac{(b-a)^5}{2880}f^{(4)}(\eta) 2880(ba)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=(ba)f(2a+b)f(x)=f(2a+b)+(x2a+b)f(2a+b)+(x2a+b)22!f′′(ξ)abf(x)dx=(ba)f(2a+b)+ab2!f′′(ξ)(x2a+b)2dxSince f′′(x) is continuous on [a,b],(x2a+b)2 unchanged sign By the median integral theorem,η[a,b]abf(x)dx=(ba)f(2a+b)+2!f′′(η)ab(x2a+b)2dxE(f)=2!f′′(η)ab(x2a+b)2dx=24(ba)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′′(ξ)(xa)(xb)dx=2!f′′(η)ab(xa)(xb)dx=12(ba)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,,n1),步长 h = b − a n h=\frac{b-a}{n} h=nba
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=0n1xixi+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=0n1xixi+1f(x)dx=i=0n12h[f(xi)+f(xi+1)]12h3i=0n1f′′(ηi)=2h[f(a)+2i=0n1+f(b)](ba)12h2[n1i=0n1f′′(η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=0n1f′′(η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=1n1f(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(ba)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=1n1f(xi)+4i=0n1f(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)=2880bah4f(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 IT2n+31(T2nTn)=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 IS2n+151(S2nSn)=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 IC2n+631(C2nCn)=Rn

Romberg求积算法

将事后误差估计 T n → S n → C n → R n T_n\rightarrow S_n \rightarrow C_n \rightarrow R_n TnSnCnRn列成表格就是 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 414Q(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)dxk=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)dxk=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)=(xx0)(xx1)(xxn)正交于空间 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,pnPn

  • 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

步骤
  1. 确定与不超过 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

  2. ω n + 1 ( x ) \omega_{n+1}(x) ωn+1(x)的所有零点

  3. 由代数精度的等价定义列方程组—待定系数法

截断误差

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)dxi=0nAif(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=(1xi2)[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)dx2f(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)dxf(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) 111x2 1f(x)dxn+1πi=0nf(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+exf(x)dxk=0nAkf(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) +ex2f(x)dxk=0nAkf(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(ba)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<xby(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)∣≤Ly1y2(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,r1kr1))

  • 二级二阶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            (n4)

单步法的收敛性、相容性与稳定性

  • 收敛性: h → 0 , y n + 1 → y ( x n + 1 ) h\rightarrow 0,y_{n+1} \rightarrow y(x_{n+1}) h0,yn+1y(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=0k1αiyni+hi=1k1βifni

  • β − 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(55fn59fn1+37fn29fn3)(Adams)

R n + 1 = 251 720 h 5 y ( 5 ) ( η ) R_{n+1} = \frac{251}{720}h^5y^{(5)}(\eta) Rn+1=720251h5y(5)(η)

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

愤怒的卤蛋

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值