精华内容
下载资源
问答
  • 数值微分
    千次阅读
    2020-09-10 17:25:08

    1. 数值微分常用公式

    \qquad 按照导数的定义,数值微分可以简单地用“差商”来近似,常用的数值微分公式有:

    • 前向

      f ′ ( x ) = f ( x + h ) − f ( x ) h \qquad\qquad\qquad f^{'}(x)=\dfrac{f(x+h)-f(x)}{h} f(x)=hf(x+h)f(x)

    • 后向

      f ′ ( x ) = f ( x ) − f ( x − h ) h \qquad\qquad\qquad f^{'}(x)=\dfrac{f(x)-f(x-h)}{h} f(x)=hf(x)f(xh)

    • 中心(中点公式)

      f ′ ( x ) = f ( x + h ) − f ( x − h ) 2 h \qquad\qquad\qquad f^{'}(x)=\dfrac{f(x+h)-f(x-h)}{2h} f(x)=2hf(x+h)f(xh)

    \qquad 其中, h h h 为步长,直接决定了数值微分的精度:

    ( 1 ) \qquad(1) (1) 从截断误差的角度——步长越小,计算结果越精确

    ( 2 ) \qquad(2) (2) 从舍入误差的角度——步长不宜过小

    \qquad   当步长 h h h 非常小, f ( x + h ) f(x+h) f(x+h) f ( x − h ) f(x-h) f(xh) 的值非常接近,直接相减可能会造成有效数字的损失

    \qquad

    2. 插值型求导公式

    \qquad 由于多项式的求导比较容易,可以利用函数插值来建立插值多项式 P n ( x ) ≈ f ( x ) P_n(x)\approx f(x) Pn(x)f(x),从而建立数值微分的公式:

    f ′ ( x ) ≈ P n ′ ( x ) \qquad\qquad\qquad f^{'}(x)\approx P_n^{'}(x) f(x)Pn(x)
    \qquad
    \qquad 依据插值余项原理, f ′ ( x ) ≈ P n ′ ( x ) f^{'}(x)\approx P_n^{'}(x) f(x)Pn(x) 的余项为:

    f ′ ( x ) − P n ′ ( x ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ω n + 1 ′ ( x ) + ω n + 1 ( x ) ( n + 1 ) ! d d x f ( n + 1 ) ( ξ ) \qquad\qquad\qquad f^{'}(x)-P_n^{'}(x)=\dfrac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}^{'}(x)+\dfrac{\omega_{n+1}(x)}{(n+1)!}\dfrac{d}{dx}f^{(n+1)}(\xi) f(x)Pn(x)=(n+1)!f(n+1)(ξ)ωn+1(x)+(n+1)!ωn+1(x)dxdf(n+1)(ξ)

    多项式插值的余项: R n ( x ) = f ( x ) − P n ( x ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ω n + 1 ( x ) , ξ ∈ ( a , b ) ,   ∀ x ∈ [ a , b ] R_n(x)=f(x)-P_n(x)=\dfrac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x),\quad\xi\in(a,b),\ \forall x\in[a,b] Rn(x)=f(x)Pn(x)=(n+1)!f(n+1)(ξ)ωn+1(x),ξ(a,b), x[a,b]

    \qquad 显然,对于任意给定的点 x x x,误差 f ′ ( x ) − P n ′ ( x ) f^{'}(x)-P_n^{'}(x) f(x)Pn(x) 无法预估。

    \qquad 如果待求解的是插值节点 x k   ( k = 0 , 1 , ⋯   , n ) x_k\ (k=0,1,\cdots,n) xk (k=0,1,,n) 上的导数值,由于 ω n + 1 ( x ) = ∏ i = 0 n ( x − x i ) \omega_{n+1}(x)=\prod\limits_{i=0}^n(x-x_i) ωn+1(x)=i=0n(xxi) 在插值节点的值 ω n + 1 ( x k ) = 0 \omega_{n+1}(x_k)=0 ωn+1(xk)=0,插值余项可以简化为

    f ′ ( x k ) − P n ′ ( x k ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ω n + 1 ′ ( x k ) \qquad\qquad\qquad f^{'}(x_k)-P_n^{'}(x_k)=\dfrac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}^{'}(x_k) f(xk)Pn(xk)=(n+1)!f(n+1)(ξ)ωn+1(xk)

    \qquad

    • 两点公式:已知两个节点 x 0 , x 1 = x 0 + h x_0,x_1=x_0+h x0,x1=x0+h 的函数值为 f ( x 0 ) f(x_0) f(x0) f ( x 1 ) f(x_1) f(x1)

    P 1 ( x ) = x − x 1 x 0 − x 1 f ( x 0 ) + x − x 0 x 1 − x 0 f ( x 1 ) \qquad\qquad P_1(x)=\dfrac{x-x_1}{x_0-x_1}f(x_0)+\dfrac{x-x_0}{x_1-x_0}f(x_1) P1(x)=x0x1xx1f(x0)+x1x0xx0f(x1)
    \qquad
    \qquad 记步长 h = x 1 − x 0 h=x_1-x_0 h=x1x0,那么在插值节点 x 0 , x 1 x_0,x_1 x0,x1 的导数值为:

    P 1 ′ ( x ) = 1 h [ − f ( x 0 ) + f ( x 1 ) ] \qquad\qquad P_1^{'}(x)=\dfrac{1}{h}[-f(x_0)+f(x_1)] P1(x)=h1[f(x0)+f(x1)]

    ⟹ \qquad\qquad \Longrightarrow { f ′ ( x 0 ) ≈ P 1 ′ ( x 0 ) = 1 h [ f ( x 1 ) − f ( x 0 ) ] f ′ ( x 1 ) ≈ P 1 ′ ( x 1 ) = 1 h [ f ( x 1 ) − f ( x 0 ) ] \quad\left\{\begin{aligned}f^{'}(x_0)\approx P_1^{'}(x_0)=\dfrac{1}{h}[f(x_1)-f(x_0)]\\ \\f^{'}(x_1)\approx P_1^{'}(x_1)=\dfrac{1}{h}[f(x_1)-f(x_0)]\end{aligned}\right. f(x0)P1(x0)=h1[f(x1)f(x0)]f(x1)P1(x1)=h1[f(x1)f(x0)]

    \qquad

    • 三点公式:已知三个等距节点 x 0 , x 1 = x 0 + h , x 2 = x 0 + 2 h x_0,x_1=x_0+h,x_2=x_0+2h x0,x1=x0+h,x2=x0+2h 的函数值为 f ( x 0 ) , f ( x 1 ) , f ( x 2 ) f(x_0),f(x_1),f(x_2) f(x0),f(x1),f(x2)

    P 2 ( x ) = ( x − x 1 ) ( x − x 2 ) ( x 0 − x 1 ) ( x 0 − x 2 ) f ( x 0 ) + ( x − x 0 ) ( x − x 2 ) ( x 1 − x 0 ) ( x 1 − x 2 ) f ( x 1 ) + ( x − x 0 ) ( x − x 1 ) ( x 2 − x 0 ) ( x 2 − x 1 ) f ( x 2 ) \qquad\qquad P_2(x)=\dfrac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+\dfrac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1)+\dfrac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2) P2(x)=(x0x1)(x0x2)(xx1)(xx2)f(x0)+(x1x0)(x1x2)(xx0)(xx2)f(x1)+(x2x0)(x2x1)(xx0)(xx1)f(x2)

    \qquad
    \qquad 令等距节点 x = x 0 + t h   ( t = 0 , 1 , 2 ) x=x_0+th\ (t=0,1,2) x=x0+th (t=0,1,2)

    P 2 ( x ) = P 2 ( x 0 + t h ) = 1 2 ( t − 1 ) ( t − 2 ) f ( x 0 ) − t ( t − 2 ) f ( x 1 ) + 1 2 t ( t − 1 ) f ( x 2 ) \qquad\qquad P_2(x)=P_2(x_0+th)=\dfrac{1}{2}(t-1)(t-2)f(x_0)-t(t-2)f(x_1)+\dfrac{1}{2}t(t-1)f(x_2) P2(x)=P2(x0+th)=21(t1)(t2)f(x0)t(t2)f(x1)+21t(t1)f(x2)

    \qquad 两端对 t t t 求导,那么在插值节点 x 0 , x 1 , x 2 x_0,x_1,x_2 x0,x1,x2 的导数值为:

    P 2 ′ ( x 0 + t h ) = 1 2 h [ ( 2 t − 3 ) f ( x 0 ) − ( 4 t − 4 ) f ( x 1 ) + ( 2 t − 1 ) f ( x 2 ) ] \qquad\qquad P_2^{'}(x_0+th)=\dfrac{1}{2h}[(2t-3)f(x_0)-(4t-4)f(x_1)+(2t-1)f(x_2)] P2(x0+th)=2h1[(2t3)f(x0)(4t4)f(x1)+(2t1)f(x2)]

    ⟹ \qquad\qquad \Longrightarrow { f ′ ( x 0 ) ≈ P 2 ′ ( x 0 ) = 1 2 h [ − 3 f ( x 0 ) + 4 f ( x 1 ) − f ( x 2 ) ] , t = 0 f ′ ( x 1 ) ≈ P 2 ′ ( x 1 ) = 1 2 h [ − f ( x 0 ) + f ( x 2 ) ] , t = 1 f ′ ( x 2 ) ≈ P 2 ′ ( x 2 ) = 1 2 h [ f ( x 0 ) − 4 f ( x 1 ) + 3 f ( x 2 ) ] , t = 2 \quad\left\{\begin{aligned}f^{'}(x_0)&\approx P_2^{'}(x_0)=\dfrac{1}{2h}[-3f(x_0)+4f(x_1)-f(x_2)],&t=0\\ f^{'}(x_1)&\approx P_2^{'}(x_1)=\dfrac{1}{2h}[-f(x_0)+f(x_2)] ,&t=1\\f^{'}(x_2)&\approx P_2^{'}(x_2)=\dfrac{1}{2h}[f(x_0)-4f(x_1)+3f(x_2)],&t=2\end{aligned}\right. f(x0)f(x1)f(x2)P2(x0)=2h1[3f(x0)+4f(x1)f(x2)],P2(x1)=2h1[f(x0)+f(x2)],P2(x2)=2h1[f(x0)4f(x1)+3f(x2)],t=0t=1t=2

    \qquad 其中 f ′ ( x 1 ) = 1 2 h [ − f ( x 0 ) + f ( x 2 ) ] f^{'}(x_1)=\dfrac{1}{2h}[-f(x_0)+f(x_2)] f(x1)=2h1[f(x0)+f(x2)] 由于少用一个函数值而应用广泛,也就是“中点公式”。

    \qquad

    • 更高阶的数值微分公式

    \qquad 令等距的插值节点 x = x 0 + t h   ( t = 0 , 1 , 2 ) x=x_0+th\ (t=0,1,2) x=x0+th (t=0,1,2),在插值节点 x 0 , x 1 , x 2 x_0,x_1,x_2 x0,x1,x2 的二阶导数值为

    P 2 ′ ′ ( x 0 + t h ) = 1 h 2 [ f ( x 0 ) − 2 f ( x 1 ) + f ( x 2 ) ] \qquad\qquad P_2^{''}(x_0+th)=\dfrac{1}{h^2}[f(x_0)-2f(x_1)+f(x_2)] P2(x0+th)=h21[f(x0)2f(x1)+f(x2)]

    \qquad 可得到:

    P 2 ′ ′ ( x 1 ) = 1 h 2 [ f ( x 1 − h ) − 2 f ( x 1 ) + f ( x 1 + h ) ] , t = 1 \qquad\qquad P_2^{''}(x_1)=\dfrac{1}{h^2}[f(x_1-h)-2f(x_1)+f(x_1+h)],\quad t=1 P2(x1)=h21[f(x1h)2f(x1)+f(x1+h)],t=1
    \qquad

    实现代码

    import numpy as np
    
    def diff1_f(f,x,h):#前向
        return (f(x+h)-f(x))/h
    
    def diff1_b(f,x,h):#后向
        return (f(x)-f(x-h))/h
    
    def diff1_c(f,x,h):#中心差分
        return (f(x+h)-f(x-h))/(2*h)
    
    def diff2_c(f,x,h):#等距节点的二阶差分
        return (f(x-h)+f(x+h)-2*f(x))/h**2
    
    if __name__ == '__main__':    
        f = lambda x: np.exp(x)
        h = 0.01
        val1 = diff1_f(f,0,h)
        val2 = diff1_b(f,0,h)
        val3 = diff1_c(f,0,h)
    

    f ( x ) = e x f(x)=e^x f(x)=ex x = 0 x=0 x=0 处在不同步长时的数值微分结果:

    f ( x ) = e x f(x)=e^x f(x)=ex h = 0.001 h=0.001 h=0.001 h = 0.01 h=0.01 h=0.01 h = 0.1 h=0.1 h=0.1
    前向差分1.00050016670838461.0050167084167951.0517091807564771
    后向差分0.99950016662497810.99501662508318930.9516258196404048
    中心差分1.00000016666668131.00001666674999211.001667500198441

    利用这些公式,可以画出一些函数的数值微分结果:
    \qquad 在这里插入图片描述

    3. 基于数值积分的方法

    \qquad 微分是积分的逆运算,可以利用积分公式来求微分值。

    \qquad φ ( x ) = f ′ ( x ) \varphi(x)=f^{'}(x) φ(x)=f(x),考虑步长 h = b − a n h=\dfrac{b-a}{n} h=nba 时的等距节点: x k = a + k h ,   k = 0 , 1 , ⋯   , n x_k=a+kh,\ k=0,1,\cdots,n xk=a+kh, k=0,1,,n,那么

    f ( x k + 1 ) − f ( x k − 1 ) = ∫ x k − 1 x k + 1 f ′ ( x ) d x = ∫ x k − 1 x k + 1 φ ( x ) d x , k = 0 , 1 , ⋯   , n \qquad\qquad\begin{aligned} f(x_{k+1})-f(x_{k-1})&=\displaystyle\int_{x_{k-1}}^{x_{k+1}}f^{'}(x)dx\\ &=\displaystyle\int_{x_{k-1}}^{x_{k+1}}\varphi(x)dx\quad,\quad k=0,1,\cdots,n\end{aligned} f(xk+1)f(xk1)=xk1xk+1f(x)dx=xk1xk+1φ(x)dx,k=0,1,,n

    • 采用中矩形公式求解积分

    ∫ a b f ( x ) d x = ( b − a ) f ( a + b 2 ) \qquad\qquad\displaystyle\int_a^bf(x)dx=(b-a)f\left(\dfrac{a+b}{2}\right) abf(x)dx=(ba)f(2a+b)

    \qquad 由于是等距节点,那么

    ∫ x k − 1 x k + 1 φ ( x ) d x = ( x k + 1 − x k − 1 ) φ ( x k + 1 + x k − 1 2 ) = 2 h φ ( x k ) = f ( x k + 1 ) − f ( x k − 1 ) \qquad\qquad\begin{aligned}\displaystyle\int_{x_{k-1}}^{x_{k+1}}\varphi(x)dx&=(x_{k+1}-x_{k-1})\varphi\left(\dfrac{x_{k+1}+x_{k-1}}{2}\right)\\ &=2h\varphi(x_k)\\ &=f(x_{k+1})-f(x_{k-1})\end{aligned} xk1xk+1φ(x)dx=(xk+1xk1)φ(2xk+1+xk1)=2hφ(xk)=f(xk+1)f(xk1)

    \qquad 实际上就是“中点公式”:

    f ′ ( x k ) = φ ( x k ) = f ( x k + 1 ) − f ( x k − 1 ) 2 h \qquad\qquad f^{'}(x_k)=\varphi(x_k)=\dfrac{f(x_{k+1})-f(x_{k-1})}{2h} f(xk)=φ(xk)=2hf(xk+1)f(xk1)

    \qquad
    \qquad

    • 采用辛普森公式求积分

    ∫ a b f ( x ) d x = b − a 6 [ f ( a ) + 4 f ( a + b 2 ) + f ( b ) ] \qquad\qquad\displaystyle\int_a^bf(x)dx=\dfrac{b-a}{6}\left[f(a)+4f\left(\dfrac{a+b}{2}\right)+f(b)\right] abf(x)dx=6ba[f(a)+4f(2a+b)+f(b)]

    \qquad 由于是等距节点,那么

    ∫ x k − 1 x k + 1 φ ( x ) d x = x k + 1 − x k − 1 6 [ φ ( x k − 1 ) + 4 φ ( x k ) + φ ( x k + 1 ) ] = h 3 [ φ ( x k − 1 ) + 4 φ ( x k ) + φ ( x k + 1 ) ] = f ( x k + 1 ) − f ( x k − 1 ) \qquad\qquad\begin{aligned}\displaystyle\int_{x_{k-1}}^{x_{k+1}}\varphi(x)dx&=\dfrac{x_{k+1}-x_{k-1}}{6}\left[\varphi(x_{k-1})+4\varphi(x_k)+\varphi(x_{k+1})\right]\\ &=\dfrac{h}{3}\left[\varphi(x_{k-1})+4\varphi(x_k)+\varphi(x_{k+1})\right]\\ &=f(x_{k+1})-f(x_{k-1})\end{aligned} xk1xk+1φ(x)dx=6xk+1xk1[φ(xk1)+4φ(xk)+φ(xk+1)]=3h[φ(xk1)+4φ(xk)+φ(xk+1)]=f(xk+1)f(xk1)

    \qquad 得到了一组递推关系式:

    f ′ ( x k − 1 ) + 4 f ′ ( x k ) + f ′ ( x k + 1 ) = 3 h [ f ( x k + 1 ) − f ( x k − 1 ) ] \qquad\qquad f^{'}(x_{k-1})+4f^{'}(x_k)+f^{'}(x_{k+1})=\dfrac{3}{h}[f(x_{k+1})-f(x_{k-1})] f(xk1)+4f(xk)+f(xk+1)=h3[f(xk+1)f(xk1)]

    \qquad
    \qquad
    \qquad 若记 φ ( x k ) = f ′ ( x k ) = m k \varphi(x_k)=f^{'}(x_k)=m_k φ(xk)=f(xk)=mk,那么

    m k − 1 + m k + m k − 1 = 3 h [ f ( x k + 1 ) − f ( x k − 1 ) ] \qquad\qquad m_{k-1}+m_k+m_{k-1}=\dfrac{3}{h}[f(x_{k+1})-f(x_{k-1})] mk1+mk+mk1=h3[f(xk+1)f(xk1)]

    \qquad 可列出“系数矩阵严格对角占优”的线性方程组(可用追赶法求解):

    [ 4 1 1 4 1 ⋱ ⋱ ⋱ 1 4 1 1 4 ] [ m 1 m 2 ⋮ m n − 2 m n − 1 ] = [ 3 h [ f ( x 2 ) − f ( x 0 ) ] − f ′ ( x 0 ) 3 h [ f ( x 3 ) − f ( x 1 ) ] ⋮ 3 h [ f ( x n − 1 ) − f ( x n − 3 ) ] 3 h [ f ( x n ) − f ( x n − 2 ) ] − f ′ ( x n ) ] \qquad\qquad\left[\begin{matrix}4&1\\1&4&1\\&\ddots&\ddots&\ddots&\\&&1&4&1\\&&&1&4\end{matrix}\right] \left[\begin{matrix}m_1\\m_2\\\vdots\\m_{n-2}\\m_{n-1}\end{matrix}\right] =\left[\begin{matrix}\frac{3}{h}[f(x_2)-f(x_0)]-f^{'}(x_0)\\\\\frac{3}{h}[f(x_3)-f(x_1)]\\\vdots\\\frac{3}{h}[f(x_{n-1})-f(x_{n-3})]\\\\\frac{3}{h}[f(x_n)-f(x_{n-2})]-f^{'}(x_n) \end{matrix}\right] 4114114114m1m2mn2mn1=h3[f(x2)f(x0)]f(x0)h3[f(x3)f(x1)]h3[f(xn1)f(xn3)]h3[f(xn)f(xn2)]f(xn)

    第一个和最后一个线性方程中包含的 f ′ ( x 0 ) f^{'}(x_0) f(x0) f ′ ( x n ) f^{'}(x_n) f(xn) 可采用“中点公式”来求解

    实现代码

    f ( x ) = x f(x)=\sqrt{x} f(x)=x x = 101 , 102 , 103 , 104 x=101,102,103,104 x=101,102,103,104 上的一阶导数

    import numpy as np
    
    def diff1_c(f,x,h):
        return (f(x+h)-f(x-h))/(2*h)
    
    def chasing(f,a,b,h):   # solve linear equations by chasing method
        xp = np.arange(a,b+h,h)
        tf = f(xp)
        bm = (tf[2:] - tf[:len(xp)-2])*3/h
        bm[0] = bm[0] - diff1_c(f,xp[0],h)
        bm[-1] = bm[-1] - diff1_c(f,xp[-1],h)
    
        n = len(bm)
        Am = np.eye(n)*4
        for i in range(n-1):
            Am[i+1,i] = 1
            Am[i,i+1] = 1
        
        L = np.zeros(Am.shape)
        L[0,0] = Am[0,0]
        U = np.eye(n)
        U[0,1] = Am[0,1]/L[0,0]
        for i in range(1,n):
            L[i,i-1] = Am[i,i-1]                     # r_i
            L[i,i] = Am[i,i] - L[i,i-1]*U[i-1,i]     # alpha_i = b_i - r_i*beta_{i-1}
            if i==n-1:
                break
            U[i,i+1] = Am[i,i+1]/L[i,i]              # beta_i = c_i/alpha_i      
           
        # Ly = bm
        rows = len(bm)
        y = np.zeros(rows)
        y[0] = bm[0]/L[0,0]
        for k in range(1,rows):
            y[k] = (bm[k] - np.sum(L[k,:k]*y[:k]))/L[k,k]   
        # Ux = y (back substitution)     
        x = np.zeros(rows)
        k = rows-1
        x[k] = y[k]/U[k,k] 
        for k in range(rows-2,-1,-1):
            x[k] = (y[k] - np.sum(x[k+1:]*U[k,k+1:]))/U[k,k]           
        return x        
    
    if __name__ == '__main__':     
        f = lambda x: np.sqrt(x)
        a = 100
        b = 105
        h = 1
        x =chasing(f,a,b,h) # 必须满足b=a+kh
        print(x)
    

    计算结果:
    [0.0497516947 0.0495074114 0.0492664916 0.0490288885]

    \qquad

    4. 数值微分的外推算法

    \qquad 利用中心差分公式计算导数值:

    f ′ ( x ) ≈ G ( h ) = f ( x + h ) − f ( x − h ) 2 h \qquad\qquad f^{'}(x)\approx G(h)=\dfrac{f(x+h)-f(x-h)}{2h} f(x)G(h)=2hf(x+h)f(xh)

    \qquad f ( x ) f(x) f(x) 做泰勒级数展开:

    f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) 2 ! ( x − x 0 ) 2 + f ( 3 ) ( x 0 ) 3 ! ( x − x 0 ) 3 + ⋯ \qquad\qquad f(x)=f(x_0)+f^{'}(x_0)(x-x_0)+\dfrac{f^{''}(x_0)}{2!}(x-x_0)^2+\dfrac{f^{(3)}(x_0)}{3!}(x-x_0)^3+\cdots f(x)=f(x0)+f(x0)(xx0)+2!f(x0)(xx0)2+3!f(3)(x0)(xx0)3+

    ⟹ { f ( x 0 + h ) = f ( x 0 ) + f ′ ( x 0 ) h + f ′ ′ ( x 0 ) 2 ! h 2 + f ( 3 ) ( x 0 ) 3 ! h 3 + ⋯ f ( x 0 − h ) = f ( x 0 ) − f ′ ( x 0 ) h + f ′ ′ ( x 0 ) 2 ! h 2 − f ( 3 ) ( x 0 ) 3 ! h 3 + ⋯ \qquad\Longrightarrow\left\{\begin{aligned} f(x_0+h)=f(x_0)+f^{'}(x_0)h+\dfrac{f^{''}(x_0)}{2!}h^2+\dfrac{f^{(3)}(x_0)}{3!}h^3+\cdots \\ f(x_0-h)=f(x_0)-f^{'}(x_0)h+\dfrac{f^{''}(x_0)}{2!}h^2-\dfrac{f^{(3)}(x_0)}{3!}h^3+\cdots \end{aligned}\right. f(x0+h)=f(x0)+f(x0)h+2!f(x0)h2+3!f(3)(x0)h3+f(x0h)=f(x0)f(x0)h+2!f(x0)h23!f(3)(x0)h3+

    \qquad 两个式子相减,可得到:

    f ( x 0 + h ) − f ( x 0 − h ) = f ′ ( x 0 ) 2 h − 2 f ( 3 ) ( x 0 ) 3 ! h 3 − 2 f ( 5 ) ( x 0 ) 5 ! h 5 + ⋯ \qquad\qquad f(x_0+h)-f(x_0-h)=f^{'}(x_0)2h-\dfrac{2f^{(3)}(x_0)}{3!}h^3-\dfrac{2f^{(5)}(x_0)}{5!}h^5+\cdots f(x0+h)f(x0h)=f(x0)2h3!2f(3)(x0)h35!2f(5)(x0)h5+

    ⟹ f ′ ( x 0 ) = f ( x 0 + h ) − f ( x 0 − h ) 2 h + 2 f ( 3 ) ( x 0 ) 3 ! h 2 + 2 f ( 5 ) ( x 0 ) 5 ! h 4 + ⋯ \qquad\qquad\Longrightarrow f^{'}(x_0)=\dfrac{f(x_0+h)-f(x_0-h)}{2h}+\dfrac{2f^{(3)}(x_0)}{3!}h^2+\dfrac{2f^{(5)}(x_0)}{5!}h^4+\cdots f(x0)=2hf(x0+h)f(x0h)+3!2f(3)(x0)h2+5!2f(5)(x0)h4+

    ⟹ f ′ ( x ) = G ( h ) + 2 f ( 3 ) ( x ) 3 ! h 2 + 2 f ( 5 ) ( x ) 5 ! h 4 + ⋯ \qquad\qquad\Longrightarrow f^{'}(x)=G(h)+\dfrac{2f^{(3)}(x)}{3!}h^2+\dfrac{2f^{(5)}(x)}{5!}h^4+\cdots f(x)=G(h)+3!2f(3)(x)h2+5!2f(5)(x)h4+

    \qquad 也就是说, f ′ ( x ) f^{'}(x) f(x) 展开成泰勒级数后可以表示为:

    f ′ ( x ) = G ( h ) + α 1 h 2 + α 2 h 4 + ⋯ \qquad\qquad f^{'}(x)=G(h)+\alpha_1h^2+\alpha_2h^4+\cdots f(x)=G(h)+α1h2+α2h4+
    \qquad
    \qquad 因此,可以采用《数值积分的 Richardson \text{Richardson} Richardson递推》(4.3节) 的思路,利用 Richardson \text{Richardson} Richardson外推公式,对初始步长 h h h 逐次分半,记 G 0 ( h ) = G ( h ) G_0(h)=G(h) G0(h)=G(h),那么数值微分的外推算法为:
    \qquad
    G m ( h ) = 4 m G m − 1 ( h 2 ) − G m − 1 ( h ) 4 m − 1   , m = 1 , 2 , ⋯ \qquad\qquad G_m(h)=\dfrac{4^mG_{m-1}(\frac{h}{2})-G_{m-1}(h)}{4^m-1}\ ,\quad m=1,2,\cdots Gm(h)=4m14mGm1(2h)Gm1(h) ,m=1,2,
    \qquad

    Richardson外推方法的误差为: f ′ ( x ) − G m ( h ) = O ( h 2 ( m + 1 ) ) f^{'}(x)-G_m(h)=O(h^{2(m+1)}) f(x)Gm(h)=O(h2(m+1))

    \qquad
    \qquad f ( x ) = x 2 e − x f(x)=x^2e^{-x} f(x)=x2ex 为例,求 f ′ ( 0.5 ) f^{'}(0.5) f(0.5) 的值,初始步长 h = 0.1 h=0.1 h=0.1

    G 0 G_0 G0 G 1 G_1 G1 G 2 G_2 G2
    G 0 ( h ) : 0.4516049081 G_0(h):0.4516049081 G0(h):0.4516049081
    G 0 ( h 2 ) : 0.4540761694 G_0(\frac{h}{2}):0.4540761694 G0(2h):0.4540761694 G 1 ( h ) : 0.4548999231 G_1(h):0.4548999231 G1(h):0.4548999231
    G 0 ( h 2 2 ) : 0.4546926288 G_0(\frac{h}{2^2}):0.4546926288 G0(22h):0.4546926288 G 1 ( h 2 ) : 0.4548981152 G_1(\frac{h}{2}):0.4548981152 G1(2h):0.4548981152 G 2 ( h ) : 0.4548979947 G_2(h):0.4548979947 G2(h):0.4548979947
    ⋮ \vdots ⋮ \vdots ⋮ \vdots

    只需要根据中心差分公式,求出第1列中步长为 h , h 2 , h 2 2 , h 2 3 , ⋯ h,\frac{h}{2},\frac{h}{2^2},\frac{h}{2^3},\cdots h,2h,22h,23h, 的微分值
    由外推算法,表中第2列开始的所有表项的值都是直接通过(与左、左上表项的)递推关系计算得到。

    \qquad 由表可知:采用步长 h = 0.1 2 2 h=\frac{0.1}{2^2} h=220.1 仅有 3 3 3 位有效数字,外推 1 1 1 次之后达到了 5 5 5 位有效数字,外推 2 2 2 次就达到了 9 9 9 位有效数字。
    \qquad

    实现代码

    import numpy as np
    import matplotlib.pyplot as plt
    
    def diff1_c(f,x,h):
        return (f(x+h)-f(x-h))/(2*h)
    
    def diff1_richardson(f,x,h,m):
        val_t = np.zeros(m+1)
        for i in range(m+1):
            val_t[i] = diff1_c(f,x,h/(2**i)) 
        np.set_printoptions(precision=10)        
        # print(val_t[-1::-1].flatten())    
        
        for i in range(1,m+1):
            t0 = 4**i
            t1 = np.concatenate((np.ones(1),val_t[0:-1]))
            t2 = (t0*val_t - t1)/(t0-1)
            val_t = t2[1:]
            # print(val_t[-1::-1].flatten())
        # print('[%.9f]' % (val_t[0]))    
        return val_t[0]     
       
    if __name__ == '__main__':    
        f = lambda x: np.exp(-x)*x**2
        h = 0.1
        x = 0.5
        val = diff1_richardson(f,x,h,2)
        print(val)
    

    f ′ ( 0.5 ) f^{'}(0.5) f(0.5) 在初始步长为 h = 0.1 h=0.1 h=0.1 时的计算结果:

    [0.4546926288 0.4540761694 0.4516049081]
    [0.4548981152 0.4548999231]
    [0.4548979947]

    可以画出一些函数的数值微分结果:
    \qquad 在这里插入图片描述
    参考
    数值积分的python实现——NewtonCotes、复化求积、Romberg、richardson递推

    更多相关内容
  • Matlab数值微分法汇总:MidPoint 中点公式求取导数ThreePoint 三点法求函数的导数FivePoint 五点法求函数的导数DiffBSample 三次样条法求函数的导数SmartDF 自适应法求函数的导数CISimpson 辛普森数值微分法法求函数...
  • 绝对实用的数值计算方法,差商型求导公式、插值型求导公式的数值微分,插值型求积公式、复化求积公式、复化求积公式的数值积分
  • 第7章 MATLAB数值微分与积分_习题答案.doc.pdf
  • 信息与计算科学专业——一篇适用于数值计算方法期末考试的复习笔记
  • 1 一元函数的数值微分可以根据函数在某点处导数的定义来实现代码,首先回顾一下函数在一点处导数的相关定义定义 设函数在点的某个邻域内有定义,当自变量处取得增量仍在
  • 计算方法第四章数值积分和数值微分.ppt
  • M a t l a b 数 值 积 分 与 数 值 微 分 Matlab 数值积分 1. 一重数值积分的实现方法 变步长辛普森法高斯 - 克朗罗德法梯形积分法 1.1 变步长辛普森法 Matlab 提供了 quad 函数和 quadl 函数用于实现变步长 辛普森...
  • 数值积分和数值微分

    2018-06-29 16:01:36
    数值求积,代数精度,插值型求积公式,求积公式的收敛性与稳定性等。
  • 计算机图形学之直线段的扫描转换C++实现 包含DDA数值微分算法,中点画线法,Bresenham算法
  • 1. 完成坐标变换,将坐标原点移动到(400,300)处,并使X轴正方向水平向右,使Y轴正方向垂直向上(原来的坐标原点位于绘图...2. 分别实现数值微分算法和Bresenham算法,不能使用系统API生成直线而只能使用SetPixel函数。
  • 向量 Y 相对于通道 T 的数值微分 (即 Delta Y / Delta T) 对任意偶数阶 N 用中心差分法。
  • MATLAB数值微分

    2013-04-16 22:34:33
    MATLAB进行数值微分程序,主要包括微分函数,适合简单的编程
  • 工程数值计算matlab实验报告——辛普森数值微分数值积分.pdf工程数值计算matlab实验报告——辛普森数值微分数值积分.pdf工程数值计算matlab实验报告——辛普森数值微分数值积分.pdf工程数值计算matlab实验报告——...
  • 工程数值计算matlab实验报告——辛普森数值微分数值积分.docx工程数值计算matlab实验报告——辛普森数值微分数值积分.docx工程数值计算matlab实验报告——辛普森数值微分数值积分.docx工程数值计算matlab实验报告...
  • 差分方程和数值微分,离散时段上描述变化过程的数学模型
  • 6.1 matlab数值微分与数值积分

    千次阅读 2021-12-10 14:58:53
    1、数值微分 (1)数值差分与差商 任意函数f(x)在x0点的导数是通过极限定义的: 如果去掉极限定义中h趋向于0的极限过程,得到函数在x0点处以h为步长的向前差分。 当步长h充分小时得到函数在x0点处,以h为步长的向前...

    数值微积分适合求解没有或很准求出微分或积分表达式的问题的计算。
    1、数值微分
    (1)数值差分与差商
    任意函数f(x)在x0点的导数是通过极限定义的:
    在这里插入图片描述
    如果去掉极限定义中h趋向于0的极限过程,得到函数在x0点处以h为步长的向前差分。
    在这里插入图片描述当步长h充分小时得到函数在x0点处,以h为步长的向前差商
    在这里插入图片描述函数f(x)在x0点的微分接近于函数在该点的差分。而f(x)在x0点的导数接近于函数在该点的差商。

    (2)数值微分的实现
    MATLAB提供了求 向前差分的函数diff,其调用格式有三种:
    (1)dx=diff(x): 计算向量x的一阶向前差分,dx(i)=x(i+1)-x(i),i=1,
    2,…,n-1。n是向量x的元素的个数

    (2)dx=diff(x,n): 计算向量x的n阶向前差分。例如,diff(x,2)=diff(diff(x)):x的二阶差分等于x的一阶差分再求一阶差分

    (3)dx=diff(A,n,dim): 计算矩阵A的n阶差分,dim=1时(默认状态)按列计算差分;dim=2,按行计算差分。

    注意:dff函数计算的是向量元素间的差分,故差分向量元素的个数比原向量少了一个;同样对于矩阵来说差分后的矩阵比原矩阵少了一行或一列。另外计算差分之后,可以用f(x)在某点处的差商作为其导数的近以值。

    例1设f(x)=sin x,在[0,2π]范围内随机采样,计算f’(x)的近似值,并与
    理论值f’(x)=cos x进行比较。

    x = [0,sort(2*pi*rand(1,5000)),2*pi];   %确定x向量,其首尾元素分别是02π,中闻是02π开区间的随机数并按照从小到大排序
    y = sin(x); 
    f1 = diff(y)./diff(x);   %利用diff函数求一阶向前差分,y的差分除以x的差分,得到差商向量f1。差商作为导数的近似值
    f2 = cos(x(1:end-1));   %求各点导数的理论值向量f2
    plot(x(1:end-1),f1,'--',x(1:end-1),f2,'r:');   %绘制近似值和理论值曲线
    legend('近似值曲线','理论值曲线','location','southeast')
    
    d = norm(f1-f2)   %求近似值和理论值之差的范数,其值较小,说明近似值和理论值比较接近
    
    %输出结果
    d =
    
        0.0439
    

    在这里插入图片描述
    2、数值积分

    (1)数值积分基本原理
    牛顿—莱布尼兹(Newton-Leibniz)公式:
    在这里插入图片描述当被积函数的原函数无法用初等函数表示或被积函数是用表格形式给出的,这时,就需要用数值解法来求定积分的近似值。求定积分的数值方法多种多样:如梯形法、辛普森法、高斯求积公式等。它们的基本思想都是将积分区间[a b]分成n个子区间[xi xi+1],这里i从1变化到n,其中x1等于a ,xn+1等于b,这样求定积分问题就分解为下面的求和问题:
    在这里插入图片描述(2)数值积分的实现
    ①基于自适应辛普森方法
    [l,n]=quad(filename,a,b,tol,trace)

    ②基于自适应Gauss-Lobatto方法
    [l,n]=quadl(filename,a,b,tol,trace)

    其中,filename是被积函数名;a和b分别是定积分的下限和上限,积分限[a,b]必须是有限的,不能为无穷大(Inf ) ; tol用来控制积分精度,默认时取tol=10-6;trace控制是否展现积分过程,若取非0则展现积分过程,取0则不展现,默认时取trace=0;返回参数l即定积分的值,n为被积函数的调用次数。

    ③基于全局自适应积分方法
    l=integral(filename,a,b)
    其中,l是计算得到的积分; filename是被积函数;a和b分别是定积分的下限和上限,积分限可以为无穷大。

    ④基于自适应高斯-克朗罗德方法 (求震荡函数的积分)
    [l,err]=quadgk(filename,a,b)

    其中,err返回近似误差范围,其他参数的含义和用法与quad函数相同。积分上下限可以是无穷大(-lnf或Inf),也可以是复数。如果积分上下限是复数,则quadgk函数在复平面上求积分。
    ⑤基于梯形积分法
    l=trapz(x,y)

    其中,向量x、y定义函数关系y=f(x)。x,y是两个等长的向量,从x1变化到xn,y对应从y1变化到yn,并且x1小于x2一直到小于xn,积分区间是从x1到xn

    例2:分别用quad函数和quadl函数求定积分的近似值,并在相同的积分精度下,比较被积函数的调用次数。
    在这里插入图片描述

    >>format long;
    f = @(x)4./(1+x.^2);   %定义被积函数
    [I1,n1] = quad(f,0,1,1e-8)  %调用quad函数求定积分,精度取10-8次方
    [I2,n2] = quadl(f,0,1,1e-8)  %调用quadl函数求定积分
    y = (atan(1)-atan(0))*4   %输出积分的理论值
    format short;
    
    
    
    %输出结果
    I1 =
       3.141592653733437
    
    n1 =
        61
    
    I2 =
       3.141592653589806
    
    n2 =
        48
    
    y =
       3.141592653589793
    

    例3:求定积分。

    在这里插入图片描述

    >>f = @(x)1./(x.*sqrt(1-log(x).^2));   %定义被积函数
    I = integral(f,1,exp(1))     %调用integral函数求定积分
    
    %输出结果
    I =
        1.5708
    

    例4:求定积分。

    在这里插入图片描述

    >>f = @(x)sin(1./x)./(x.^2);   %定义被积函数
    I = quadgk(f,2/pi,+Inf)   %调用函数quadgk求定积分
    
    
    %输出结果
    I =
        1.0000
    

    (3)多重定积分的数值求解
    ①求二重积分的数值解:
    在这里插入图片描述

    l=integral2(filename,a,b,c,d)
    l=quad2d(filename,a,b,c,d)
    l=dblquad(filename,a,b,c,d,tol)

    ②求三重积分的数值解:
    在这里插入图片描述
    l=integral3(filename,a,b,c,d,e,f)
    l=triplequad(filename,a,b,c,d,e,f,tol)

    例6:分别求二重积分和三重积分。

    >> f1 = @(x,y)sin(x.^2+y).*exp(-x.^2/2);   %定义被积函数f1
    f2 = @(x,y,z)4*x.*z.*exp(-z.^2.*y-x.^2);  %定义被积函数f2
    I1 = quad2d(f1,-2,2,-1,1)   %调用函数quad2d求二重定积分
    I2 = integral3(f2,0,pi,0,pi,0,1)   %调用函数integral3求三重定积分
    
    
    %输出结果
    I1 =
        1.5745
    
    I2 =
        1.7328
    
    展开全文
  • 数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx数值分析第四章外推法计算数值微分...
  • 数值分析第四章外推法计算数值微分MATLAB计算实验报告.pdf数值分析第四章外推法计算数值微分MATLAB计算实验报告.pdf数值分析第四章外推法计算数值微分MATLAB计算实验报告.pdf数值分析第四章外推法计算数值微分MATLAB...
  • 数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx.pdf数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx.pdf数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx.pdf数值分析第四章外推法...
  • 数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx.docx数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx.docx数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx.docx数值分析第四章外推法...
  • 这是一个 GUI,它在多个等距点上执行函数的数值微分。 还有一个代码,用于授予用于数值微分的系数。 图片和示例应该足以了解如何使用该文件。 例子: npoints=3; 订单=1; d=datnum(npoints,Order) d= -1.5 2 -...
  • [数值微分]数值微分的误差分析

    千次阅读 2019-09-26 18:22:43
    上篇文章中提到不知道时间间隔...具体内容资料很多不再赘述,然而数值微分对于误差采用的是截断或者舍入的方法,因此数值微分一定存在误差。 然而误差的绝对存在不代表近似值一定小于真实值,随着h值的减小,截断...

    上篇文章中提到不知道时间间隔deltaT的微分形式是如何计算出来的,多方查找找到了答案,如下记录。

    对于离散数据求微分,不能直接使用求导公式来计算,需要利用数值方法中的知识。即数值微分(可参考《数值方法》[美]安妮·戈林鲍姆等)。具体内容资料很多不再赘述,然而数值微分对于误差采用的是截断或者舍入的方法,因此数值微分一定存在误差。
    三种数值微分方法
    然而误差的绝对存在不代表近似值一定小于真实值,随着h值的减小,截断误差f(h)与舍入误差g(h)会逐渐接近,当两个近似值相近时,数值微分达到最大精度,可以证明,此时误差约等于机器精度的平方根,当h继续减小,由于舍入误差的存在,此时误差会逐渐扩大。用matlab可以计算如下。
    在这里插入图片描述
    由于我电脑中未安装matlab,用python计算了在sin(0)处随h逐渐减小的数值微分结果。

    import pandas as pd
    import numpy as np
    
    
    def get_diff(h):
    	x = np.arange(0,1,h)
    	y = np.sin(x)
    	df = np.diff(y)
    	res = df/h
    	return res[0]-1
    def print_diff(h):
    	for i in range(10):
    		h=h/10
    		print(get_diff(h))
    		
    print_diff(1)
    

    或者

    import pandas as pd
    import numpy as np
    
    
    def get_diff(h):
    	t = np.sin(0+h)-np.sin(0)
    	res = t/h
    	return res-1
    def print_diff(h):
    	for i in range(10):
    		h=h/10
    		print(get_diff(h))
    		
    print_diff(1)
    

    两种结果一致,思路略有差别,另外不能使用yield。原理略过。
    附上结果
    在这里插入图片描述

    展开全文
  • 数值微分和数值积分数值微分(1)数值差分与差商(2)数值微分的实现数值积分原理函数调用 数值微分 (1)数值差分与差商 可以用差商近似计算导数值: (2)数值微分的实现 diff() 向前差分函数 dx=diff(x):...

    数值微分

    (1)数值差分与差商

    在这里插入图片描述
    在这里插入图片描述

    可以用差商近似计算导数值:在这里插入图片描述

    (2)数值微分的实现

    diff() 向前差分函数

    • dx=diff(x):计算向量x的一阶向前差分,dx[i]=x[i+1]-x[i],i=1,2,3…n-1.
    • dx=diff(x,n):计算向量x的n阶向前差分。
    • dx=diff(A,n,dim):计算矩阵A的n阶差分。dim=1时按列计算差分,=2时按行计算差分。

    举个例子
    在这里插入图片描述

    x=[0,sort(2*pi*rand(1,5000)),2*pi];
    y=sin(x);
    f1=diff(y)./diff(x)%利用差商来近似导数
    f2=cos(x(1:end-1))%实际导数
    plot(x(1:end-1),f1,x(1:end-1),f2,'r:')
    

    两者较接近
    在这里插入图片描述

    数值积分

    原理

    在这里插入图片描述

    函数调用

    • [l,n]=quad(filename,a,b,tol,trace)基于自适应辛普森方法
    • [l,n]=quadl(filename,a,b,tol,trace)基于自适应GL方法(效率比较高)
      参数含义:在这里插入图片描述
    • I=integral(filename,a,b) 基于全局自适应积分方法。I是计算得到的积分。
    • [I,err]=quadgk(filename,a,b) 基于自适应高斯—克朗罗得方法。err返回近似误差。
    • I=trapz(x,y)x自变量,y因变量。在这里插入图片描述
      举个例子
      在这里插入图片描述
    format long
    f=@(x) 4./(1+x.^2);
    [I,n]=quad(f,0,1,1e-8)
    [I,n]=quadl(f,0,1,1e-8)
    (atan(1)-atan(0))*4
    结果
    I =3.141592653733437
    n =61
    I =3.141592653589806
    n =48
    ans =3.141592653589793
    

    在这里插入图片描述

    f=@(x) 1./(x.*sqrt(1-log(x).^2));
    I=integral(f,1,exp(1))
    结果:
    I =1.570796326795570
    

    在这里插入图片描述

    f=@(x) sin(1./x)./x.^2;
    I=quadgk(f,2/pi,+Inf)
    结果:
    I = 1.000000000000000
    

    在这里插入图片描述

    x=1:6;
    y=[6 8 11 7 5 2];
    plot(x,y,'-ko');
    grid on
    axis([1 6 0 11]);
    I1=trapz(x,y)
    I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)
    结果:
    I1 =35
    I2 =35
    

    下图是x—y图像
    在这里插入图片描述

    展开全文
  • 数值分析梯形法的matlab代码数值微分.jl 关于这个包 该模块使用 Tikhonov 或 Total Variation 方法的正则化来实现潜在噪声输入的数值微分。 它需要通过多种方式进行扩展: 我无法通过预处理在大型系统上为 TV 方法做...
  • 数值分析第七章知识点总结——数值积分与数值微分,PDF版知识点总结 具体内容详见:https://blog.csdn.net/qq_36770651/article/details/109837575
  • matlab数值微分

    千次阅读 2018-10-03 11:35:46
    文章目录标签(空格分隔): matlab 数值 微分 算法@[toc]matlab数值微分1 测试数值微分函数midD2Diff()2 中心二阶差商微分 midD2Diff()3 中心差商微分 midDiff()4 前向差商微分 forwardDiff()5 后向差商微分 ...
  • 数值计算求微分,采用辛普森数值微分法法求函数的导数。
  • Matlab实现数值微分

    2021-05-03 10:33:52
    基本证明就不说啦,基本上就是运用泰勒公式实现对函数的近似运算。对于f(x)=-u‘’(x),u=sinx这个例子,代码如下a(1)=2;a(2)=-1;for i=3:99a(i)=0;endfor i=1:97for j=1:99C(i,j)=0;endendfor i=1:97C(i,i)=-1;...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 32,569
精华内容 13,027
关键字:

数值微分