精华内容
下载资源
问答
  • 数值微分的python实现

    2020-09-10 17:25:08
    1. 数值微分常用公式 2. 插值型求导公式 实现代码 3. 基于数值积分的方法 实现代码 4. 数值微分的外推算法 实现代码

    1. 数值微分常用公式

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

    • 前向

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

    • 后向

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

    • 中心(中点公式)

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

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

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

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

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

    \qquad

    2. 插值型求导公式

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

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

    f(x)Pn(x)=f(n+1)(ξ)(n+1)!ωn+1(x)+ωn+1(x)(n+1)!ddxf(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)

    多项式插值的余项: Rn(x)=f(x)Pn(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]

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

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

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

    \qquad

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

    P1(x)=xx1x0x1f(x0)+xx0x1x0f(x1)\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)
    \qquad
    \qquad记步长 h=x1x0h=x_1-x_0,那么在插值节点 x0,x1x_0,x_1 的导数值为:

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

    \qquad\qquad \Longrightarrow{f(x0)P1(x0)=1h[f(x1)f(x0)]f(x1)P1(x1)=1h[f(x1)f(x0)]\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.

    \qquad

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

    P2(x)=(xx1)(xx2)(x0x1)(x0x2)f(x0)+(xx0)(xx2)(x1x0)(x1x2)f(x1)+(xx0)(xx1)(x2x0)(x2x1)f(x2)\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)

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

    P2(x)=P2(x0+th)=12(t1)(t2)f(x0)t(t2)f(x1)+12t(t1)f(x2)\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)

    \qquad两端对 tt 求导,那么在插值节点 x0,x1,x2x_0,x_1,x_2 的导数值为:

    P2(x0+th)=12h[(2t3)f(x0)(4t4)f(x1)+(2t1)f(x2)]\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)]

    \qquad\qquad \Longrightarrow{f(x0)P2(x0)=12h[3f(x0)+4f(x1)f(x2)],t=0f(x1)P2(x1)=12h[f(x0)+f(x2)],t=1f(x2)P2(x2)=12h[f(x0)4f(x1)+3f(x2)],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.

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

    \qquad

    • 更高阶的数值微分公式

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

    P2(x0+th)=1h2[f(x0)2f(x1)+f(x2)]\qquad\qquad P_2^{''}(x_0+th)=\dfrac{1}{h^2}[f(x_0)-2f(x_1)+f(x_2)]

    \qquad可得到:

    P2(x1)=1h2[f(x1h)2f(x1)+f(x1+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
    \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)=exf(x)=e^xx=0x=0 处在不同步长时的数值微分结果:

    f(x)=exf(x)=e^x h=0.001h=0.001 h=0.01h=0.01 h=0.1h=0.1
    前向差分 1.0005001667083846 1.005016708416795 1.0517091807564771
    后向差分 0.9995001666249781 0.9950166250831893 0.9516258196404048
    中心差分 1.0000001666666813 1.0000166667499921 1.001667500198441

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

    3. 基于数值积分的方法

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

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

    f(xk+1)f(xk1)=xk1xk+1f(x)dx=xk1xk+1φ(x)dx,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}

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

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

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

    xk1xk+1φ(x)dx=(xk+1xk1)φ(xk+1+xk12)=2hφ(xk)=f(xk+1)f(xk1)\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}

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

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

    \qquad
    \qquad

    • 采用辛普森公式求积分

    abf(x)dx=ba6[f(a)+4f(a+b2)+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]

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

    xk1xk+1φ(x)dx=xk+1xk16[φ(xk1)+4φ(xk)+φ(xk+1)]=h3[φ(xk1)+4φ(xk)+φ(xk+1)]=f(xk+1)f(xk1)\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}

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

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

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

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

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

    [4114114114][m1m2mn2mn1]=[3h[f(x2)f(x0)]f(x0)3h[f(x3)f(x1)]3h[f(xn1)f(xn3)]3h[f(xn)f(xn2)]f(xn)]\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]

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

    实现代码

    f(x)=xf(x)=\sqrt{x}x=101,102,103,104x=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(xh)2h\qquad\qquad f^{'}(x)\approx G(h)=\dfrac{f(x+h)-f(x-h)}{2h}

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

    f(x)=f(x0)+f(x0)(xx0)+f(x0)2!(xx0)2+f(3)(x0)3!(xx0)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(x0+h)=f(x0)+f(x0)h+f(x0)2!h2+f(3)(x0)3!h3+f(x0h)=f(x0)f(x0)h+f(x0)2!h2f(3)(x0)3!h3+\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.

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

    f(x0+h)f(x0h)=f(x0)2h2f(3)(x0)3!h32f(5)(x0)5!h5+\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)=f(x0+h)f(x0h)2h+2f(3)(x0)3!h2+2f(5)(x0)5!h4+\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(x)=G(h)+2f(3)(x)3!h2+2f(5)(x)5!h4+\qquad\qquad\Longrightarrow f^{'}(x)=G(h)+\dfrac{2f^{(3)}(x)}{3!}h^2+\dfrac{2f^{(5)}(x)}{5!}h^4+\cdots

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

    f(x)=G(h)+α1h2+α2h4+\qquad\qquad f^{'}(x)=G(h)+\alpha_1h^2+\alpha_2h^4+\cdots
    \qquad
    \qquad因此,可以采用《数值积分的Richardson\text{Richardson}递推》(4.3节) 的思路,利用Richardson\text{Richardson}外推公式,对初始步长 hh 逐次分半,记 G0(h)=G(h)G_0(h)=G(h),那么数值微分的外推算法为:
    \qquad
    Gm(h)=4mGm1(h2)Gm1(h)4m1 ,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
    \qquad

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

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

    G0G_0 G1G_1 G2G_2
    G0(h):0.4516049081G_0(h):0.4516049081
    G0(h2):0.4540761694G_0(\frac{h}{2}):0.4540761694 G1(h):0.4548999231G_1(h):0.4548999231
    G0(h22):0.4546926288G_0(\frac{h}{2^2}):0.4546926288 G1(h2):0.4548981152G_1(\frac{h}{2}):0.4548981152 G2(h):0.4548979947G_2(h):0.4548979947
    \vdots \vdots \vdots

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

    \qquad由表可知:采用步长 h=0.122h=\frac{0.1}{2^2} 仅有 33 位有效数字,外推 11 次之后达到了 55 位有效数字,外推 22 次就达到了 99 位有效数字。
    \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) 在初始步长为 h=0.1h=0.1 时的计算结果:

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

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

    展开全文
  • 高斯勒让德数值积分公式

    千次阅读 2019-12-15 23:06:00
    高斯勒让德数值积分公式1 引言2 高斯积分公式2.1 一维区间上高斯数值积分公式2.2 二维三角形上的高斯数值积分...在插值型积分公式中,高斯积分公式具有最高的代数精度,所以是常用数值积分公式。 2 高斯积分公式...

    1 引言

    在数值求解微分方程的时候,例如利用有限元方法求解微分方程时,经常会遇到计算给定区域上的积分的问题。 在实际问题中,由于被积函数以及积分区域的复杂性,通常采用数值方法进行求解。在插值型积分公式中,高斯积分公式具有最高的代数精度,所以是常用的数值积分公式。

    2 高斯积分公式

    2.1 一维区间上高斯数值积分公式

    区间[a,b][a,b]上带权函数的插值型数值积分公式的一般形式为
    abρ(x)f(x)dxi=0nAif(xi):=In \int^b_a\rho(x)f(x)dx\approx\sum^n_{i=0}A_if(x_i):=I_n
    其中xi(i=0,1,2,,n)x_i(i=0,1,2,\cdots,n)称为插值节点。

    这里我们需要先介绍下代数精度的概念。

    定义(代数精度)若数值积分公式 In=i=0nAif(xi)I_n=\sum^n_{i=0}A_if(x_i)mm次多项式精确成立,并且存在一个m+1m+1次多项式,使得数值积分公式不精确成立,则称该数值积分公式的代数精度为mm.

    由于连续函数可以由多项式函数一致逼近,我们可以认为代数精度越高的数值积分公式,近似效果越好。对于有n+1n+1个节点的数值积分公式InI_n, 通过恰当地选择节点xix_i和系数AiA_i, 能使InI_n的代数精度达到2n+12n+1. 如果InI_n的代数精度为2n+12n+1,那么称InI_n高斯数值积分公式。在这里集中讨论ρ(x)=1\rho(x)=1的情况,相应的数值积分公式为高斯勒让德积分公式

    现在的问题是如何得到高斯积分公式的积分节点xix_i和系数AiA_i。直接利用高斯公式的代数精度为2n+12n+1,求解关于xi,Aix_i, A_i的非线性方程组
    abxkdx=i=0nAixik,k=0,1,2,,2n+1(1) \int^b_ax^kdx=\sum^n_{i=0}A_ix^k_i, k=0,1,2,\cdots,2n+1 \tag{1}
    是非常困难的。 但是幸运的是,可以证明正交多项式的零点高斯积分公式的积分节点(称为高斯点),所以可以先通过正交多项式的零点,得到高斯点xix_i, 此时(1)式变为关于AiA_i的线性方程组,从而很容易求出。

    我们知道勒让德多项式是区间[1,1][-1,1]上的正交多项式,我们首先要将区间[a,b][a,b]上的积分化为[1,1][-1,1]上。这并不困难,通过变换
    x=ba2t+a+b2 x=\frac{b-a}{2}t+\frac{a+b}{2}
    即可。勒让德多项式Pn+1(x)P_{n+1}(x)的零点即是数值积分公式(2)
    11f(x)dxi=0nAif(xi)(2) \int^1_{-1}f(x)dx\approx\sum^n_{i=0}A_if(x_i)\tag{2}
    的高斯点。低阶的高斯勒让德积分公式的节点和系数如下表所示:(图片截自百度文库https://wenku.baidu.com/view/a767f6878762caaedd33d4f4.html)

    在这里插入图片描述

    2.2 二维三角形上的高斯数值积分公式

    对于任意的三角形KK,总是可以通过面积坐标,将其映射到一个标准的三角形K^\hat{K}, 如图所示在这里插入图片描述
    所以下面我们仅讨论在K^\hat{K}上的高斯积分公式。显然有
    I=K^f(x,y)dxdy=01dx01xf(x,y)dy I=\int_{\hat{K}}f(x,y)dxdy=\int^1_0dx\int^{1-x}_0f(x,y)dy
    做变量替换u=x,v=y1xu=x, v=\frac{y}{1-x}, 则有0u,v10\leq u,v\leq 1, (x,y)(u,v)=1u\frac{\partial(x,y)}{\partial(u,v)}=1-u, 那么
    I=0101f(u,(1u)v)(1u)dudv I=\int^1_0\int^1_0f(u,(1-u)v)(1-u)dudv
    进一步做变换
    u=1+ξ2,v=1+η2 u=\frac{1+\xi}{2}, v=\frac{1+\eta}{2}
    则有1ξ,η1-1\leq \xi,\eta\leq 1, (u,v)(ξ,η)=14\frac{\partial(u,v)}{\partial(\xi,\eta)}=\frac{1}{4}, 现在有
    I=1111f(1+ξ2,(1ξ)(1+η)4)1ξ8dξη I=\int^1_{-1}\int^1_{-1}f(\frac{1+\xi}{2},\frac{(1-\xi)(1+\eta)}{4})\frac{1-\xi}{8}d\xi\eta
    现在我们已经将积分区域映射到[1,1]×[1,1][-1,1]\times[-1,1], 这样三角形上的积分就化为[1,1][-1,1]上的二重积分,从而可以利用一维高斯积分公式做计算。


    参考资料
    1.《计算方法》,孙文瑜,杜其奎,陈金如. 科学出版社.
    2. H. T. Rathod, K. V. Nagaraja, B. Venkatesudu and N. L. Ramesh, “Gauss Legendre quadrature over a triangle”. J. Indian Inst. Sci., Sept.–Oct. 2004, 84, 183–188.

    展开全文
  • 直观地说,对于一个给定的正实值函数,在一个实数区间上的定积分可以理解为在坐标平面上,由曲线、直线以及轴围成的曲边梯形的面积值(一种确定的实数值)。积分的一个严格的数学定义由波恩哈德·黎曼给出(参见条目...

    http://blog.csdn.net/pipisorry/article/details/52200140

    微积分

    直观地说,对于一个给定的正实值函数,在一个实数区间上的定积分可以理解为在坐标平面上,由曲线、直线以及轴围成的曲边梯形的面积值(一种确定的实数值)。积分的一个严格的数学定义由波恩哈德·黎曼给出(参见条目“黎曼积分”)。

     一.基本初等函数求导公式

    函数的和、差、积、商的求导法则

    反函数求导法则

    复合函数求导法则

    皮皮blog

     

     

    二、基本积分表

      

    皮皮blog

     

     

    常用凑微分公式

    [常用的求导和定积分公式(完美)]

    分部积分

    不定积分的分部积分

     

    [分部积分法]

    定积分的分部积分

     

    皮皮blog

     

     

    微分方程

     

     

    级数收敛与发散

     

    发散级数

    收敛级数

    皮皮blog

     

    微分中值定理

    f(x)为连续且光滑,任取其上两点(a, f(a))(b, f(b))a < b,那么在这两端点之间必定存在一点(c, f(c)), a < c < b,使得过c的切线斜率等于该二端点割线的斜率,即f'(c) = \frac{f(b) - f(a)}{b - a}

    [wikipedia]

    from: http://blog.csdn.net/pipisorry/article/details/52200140

    ref: 

    展开全文
  • 常用微分运算法则

    万次阅读 2015-09-28 22:16:39
    机器学习涉及到较多的数学知识,在工程应用领域,这些数学知识不是必要的,其实...本文直接给出常用的微分运算法则,并运用这些法则来计算分类回归算法 (Logistic Regression) 预测模型 Sigmoid Function 的微分公式

    机器学习涉及到较多的数学知识,在工程应用领域,这些数学知识不是必要的,其实很多算法都是数值运算专家写好了的。然而知其然知其所以然,了解这些数学公式的来龙去脉是帮助理解算法的关键。本文直接给出常用的微分运算法则,并运用这些法则来计算分类回归算法 (Logistic Regression) 预测模型 Sigmoid Function 的微分公式。

    基础函数的微分运算法则

    • 幂函数法则
      ddxxn=nxn1
    • 指数函数法则
      ddxex=ex

      ddxax=ln(a)ax
    • 对数函数法则
      ddxln(x)=1x

      ddxloga(x)=1xln(a)
    • 三角函数法则
      ddxsin(x)=cos(x)

      ddxcos(x)=sin(x)

      ddxtan(x)=sin2(x)=1cos2(x)=1+tan2(x)
    • 反三角函数法则
      ddxarcsin(x)=11x2,1<x<1

      ddxarccos(x)=11x2,1<x<1

      ddxarctan(x)=11+x2

    组合函数的微分运算法则

    • 常数法则:如果 f(x)=n,n 是常数,则
      f=0
    • 加法法则
      (αf+βg)=αf+βg
    • 乘法法则
      (fg)=fg+fg
    • 除法法则
      (fg)=fgfgg2

      根据除法法则和指数法则,可以得出推论
      ddxex=ddx1ex=0exe2x=1ex=ex
    • 链接法则:如果 f(x)=h(g(x)),则
      f(x)=h(g(x))g(x)

    计算 Sigmoid Function 的微分

    g(x)=11+ex 是分类算法的预测函数,也称为 Sigmoid Function 或 Logistic Function。我们利用上文介绍的微分运算法则来证明 Sigmoid Function 的一个特性:

    ddxg(x)=g(x)(1g(x))

    方法一

    假设 f(x)=1x,则 f(g(x))=1g(x),根据除法法则得到

    f(g(x))=(1g(x))=1g(x)1g(x)g(x)2=g(x)g(x)2

    其中 (17) 是根据除法法则得出的结论,除数是常数函数 1,被除数是 g(x)。(18) 是根据常数法则得出的结论。

    另一方面,f(g(x))=1g(x)=1+ex,根据指数法则直接计算微分得到

    f(g(x))=ddx(1+ex)=ex=11g(x)=g(x)1g(x)

    (18) 和 (22) 两式是相等的,即

    g(x)g(x)2g(x)=g(x)1g(x)=g(x)(1g(x))

    这样就得到了我们的结果。

    方法二

    g(x)=11+ex 的定义可知

    (1+ex)g(x)=1ddx((1+ex)g(x))=0exg(x)+(1+ex)ddxg(x)=0ddxg(x)=g(x)ex1+exddxg(x)=g(x)(1+ex)11+exddxg(x)=g(x)[111+ex]ddxg(x)=g(x)(1g(x))

    (26) 两边取微分;(27) 根据微分的乘法法则。

    方法三

    根据除法法则直接计算微分:

    ddxg(x)=ddx(11+ex)=0(ex)(1+ex)2=ex(1+ex)2=1(1+ex)ex(1+ex)=1(1+ex)(1+ex)1(1+ex)=1(1+ex)[11(1+ex)]=g(x)(1g(x))

    (33) 是根据除法法则得出的,其中除数是常数 1,被除数是 1+ex

    参考资料

    展开全文
  • 欧拉法简介几何意义证明泰勒展开近似求导近似积分近似几种欧拉方式向前欧拉公式向后欧拉公式梯形公式改进欧拉法截断误差...该系列主要介绍一些常用的常微分方程的数值解法,主要包括: [常微分方程的数值解法系列一]
  • 实验五、常微分方程初值问题的数值解法 ...特别是四阶龙格-库塔公式,易于编程计算,精度较高,是常用的工程计算方法。线性多步方法是在用插值多项式代替被积函数的基础上构造的,它可利用前面若干步计算结果的
  • 导数与积分公式

    2020-10-28 08:51:33
    1 极限公式(系数不为0的情况) ...10 下列常用微分公式 11 补充下面几个积分公式 12 分部积分法公式 分部积分公式: 13 第二换元积分法中的三角换元公式 14 特殊角的三角函数值 15 三角函数公式
  • 龙格---库塔方法是求解微分方程比较常用的方法,在理解数学上是怎么一回事后,编制这个程是相当容易的,就是个迭代的过程.步长的选取也是很有讲究的,过小的步长反而会导致误差累积过大.相关的理论请参考相关的数值算法...
  • 微积分 直观地说,对于一个给定的正实值函数,在一个实数区间上的定积分可以理解为在坐标平面上,由曲线、直线以及轴围成的曲边梯形的面积值(一种确定的实数值)。积分的一个严格的数学定义...
  • 求和 中图分类号O1 文献标识码A 文章编号 1674-6708(2011)35-0091-02无穷级数求和是高等数学的一个重要组成部分,它在函数的表达�研究函数的性质�求函数值以及求解微分方程等都是非常有用的,而收敛级数的求和在...
  • CotesNewton−Cotes公式:三、复化求积法:四、RombergRombergRomberg公式 / RichardsonRichardsonRichardson外推加速法:五、Gauss−LegendreGauss-LegendreGauss−Legendre公式:六、数值微分常用方法: ...
  • 微积分公式

    2018-03-21 13:45:27
    微积分直观地说,对于一个给定的正实值...基本初等函数求导公式函数的和、差、积、商的求导法则反函数求导法则复合函数求导法则皮皮blog二、基本积分表 皮皮blog常用微分公式[常用的求导和定积分公式(完美)]分...
  • 数值微分是用离散的方法近似地计算函数y=f(x)在某点x=a处的导数值,前差公式: 后差公式: 中心差商: matlab中可以用diff函数实现前差公式数值微分: 函数形式:diff(X)求解X的一阶微分,X是一个向量也...
  • 目录1 常用数值计算方法2 现代智能计算方法2.1 分类2.2 遗传算法2.3 神经网络3 科技计算软件4 习题 1 常用数值计算方法 ...(4)数值微分求解 (5)常微分方程的数值求解 • Euler公式 • 后退的Euler法 • 改进的E
  • 三角函数的有理逼近习题六第七章 数值积分与数值微分§7.1 插值型求积公式与代数度§7.2 复合求积公式及算法§7.3 外推原理与龙贝格算法§7.4 高斯型求积公式及其复合公式§7.5 数值微分应用:通信卫星覆盖地球...
  • 1. 有限差分 Taylor Table这是在数值微分常用的技巧。可以用等距网格中给定几个点的线性组合表示某一点处的微分,乃至高阶微分,并可顺便求出该差分方法的精度。 是关于 的函数,在 轴上有 个等距网格,记 我们...
  • MATLAB常用算法

    热门讨论 2010-04-05 10:34:28
    CISimpson 辛普森数值微分法法求函数的导数 Richason 理查森外推算法求函数的导数 ThreePoint2 三点法求函数的二阶导数 FourPoint2 四点法求函数的二阶导数 FivePoint2 五点法求函数的二阶导数 Diff2BSample 三次...
  • CISimpson 辛普森数值微分法法求函数的导数 Richason 理查森外推算法求函数的导数 ThreePoint2 三点法求函数的二阶导数 FourPoint2 四点法求函数的二阶导数 FivePoint2 五点法求函数的二阶导数 Diff2BSample 三次...
  • ​ 插值法是函数逼近的一种重要方法,它是数值积分、微分方程数值解等数值计算的基础与工具,其中多项式插值是最常用和最基本的方法。拉格朗日插值多项式的优点是表达式简单明确,形式对称,便于记忆,它的缺点是...
  • 量子物理学中的常用算法与程序

    热门讨论 2011-07-01 10:26:14
    函数插值与微商,常微分方程,数值积分,本征问题,递推与迭代,蒙特卡罗方法,快速傅里叶变换等内容,同时还有与内容相对应的大小程序82个。  《量子物理学中的常用算法与程序》既可作为计算物理的入门书,又可...
  • 谱方法的数值分析 !!

    2009-09-02 11:31:55
    谱方法的数值分析 全文目录 前言 第一章 预备知识 1、1Hilbert空间和Banach空间初步 1、1、1基本概念 1、1、2投影定理 1、1、3Riesz表现定理 1、1、4线性算子 1、2Sobolev空间简介 1、2、1广义导数 1、2、2Sobolev...
  • 直线绘制算法

    2019-01-15 20:58:13
    数值微分法 中点画线法 Bresenham算法   三种常用直线绘制算法: 1)数值微分法(Digital Differential Analyzer,DDA) 2)中点画线法 3)Bresenham算法     数值微分法 1. 应用直线公式:F(x) = y =...
  • 1.灰色预测法 使用范围 数据样本点个数少,6-15个 数据呈现指数或曲线的形式 ...3.求解初始微分方程所对应的时间响应函数,得到数列预测的基础公式,对一次累加生成数列的预测值,可以求得原始数的还原值 4.计算原
  • Matlab实验 微分积分和常微分方程 数值积分常用的方法 1用不定积分计算定积分 牛顿莱布尼兹公式,但是有时得不到原函数 2定义法取近似和的极限 高等数学中不是重点内容 但数值积分的各种算法却是基于定义建立的 数值...
  • Matlab实验 微分积分和常微分方程 数值积分常用的方法 1用不定积分计算定积分 牛顿莱布尼兹公式,但是有时得不到原函数 2定义法取近似和的极限 高等数学中不是重点内容 但数值积分的各种算法却是基于定义建立的 数值...
  • 4.1 常微分方程的数值解法 ;一. 数值求解的基本概念;取前两项近似; 欧拉法的特点导出简单几何意义明显便于理解能说明构造数值解法一般计算公式的基本思想通常用它来说明有关的基本概念 ;例 设系统方程为 ;则;已知...
  • 18·4 数值微分 18·5 数值积分 18·6 常微分方程的数值解法 18·7 方程的近似解 18·8 解线性方程组的直接方法 18·9 解线性方程组的迭代法 18·10 矩阵的特征值与特征向量计算 第十九章 组合论 19·1 生成函数 19...
  • 三角函数是以角度(数学上最常用弧度制,下同)为自变量,角度对应任意角终边...在数学分析中,三角函数也被定义为无穷级数或特定微分方程的解,允许它们的取值扩展到任意实数值,甚至是复数值。 三角函数十组诱导公式

空空如也

空空如也

1 2
收藏数 36
精华内容 14
热门标签