精华内容
下载资源
问答
  • 常微分方程初值问题数值解法
    千次阅读
    2020-04-27 19:05:31

    一、目的与要求
    (一)目的通过设计、编制、调试1~2个求常微分方程初值问题的数值解解的程序,加深对其数值计算方法及有关的基础理论知识的理解。
    (二)要求 用编程语言实现用改进的欧拉(Euler)公式求解常微分方程初值问题、用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题的程序。
    二、示例
    1、问题用改进的欧拉(Euler)公式求解常微分方程初值问题。
    2、算法描述(略)
    3、程序中变量说明 (略)
    4、源程序清单及运行结果(略)
    5、按以上4点要求编写上机实验报告。
    三、实验题用编程语言编程实现以下算法:
    1、用改进的欧拉(Euler)公式求解常微分方程初值问题。2、用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题。
    数值分析实验真的都是直接套公式即可…
    做法:首先理解书上关于欧拉公式以及四阶龙贝格库塔方法,明白其中的道题,模拟程序流程。
    只是图画的有点丑…
    代码

    import numpy as np
    import matplotlib.pyplot as plt
    import math
    def func(x,y):
        return y - 2*x/y
    #模拟常微分方程为:
    #      f'(x)=y-2*x/y
    #      f(x0)=y0
    #另外:解析法可知f(x) = 根号下(2*x+1)
    x0 = 0
    y0 = 1
    h = 0.1
    N = 10
    # 改进的欧拉(Euler)公式
    def ImproEuler():
        x1 = np.zeros(N+1)
        y1 = np.zeros(N+1)
        for i in range(N+1):
            x1[i] = x0 + h*(i)
        # 计算y数组各值
        for i in range(N+1):
            if i==0:
                y1[i] = y0
            else:
                k1 = y1[i-1] + h*func(x1[i-1],y1[i-1])
                k2 = y1[i-1] + h*func(x1[i],k1)
                y1[i] = (k1+k2)*1/2
        return x1,y1
    # 四阶龙格-库塔(Runge-Kutta)方法
    def RungeKutta():
        x1 = np.zeros(N+1)
        y1 = np.zeros(N+1)
        for i in range(N+1):
            x1[i] = x0 + h*(i)
        for i in range(N+1):
            if i==0:
                y1[i] = y0
            else:
                k1 = func(x1[i-1],y1[i-1])
                k2 = func(x1[i-1]+h/2,y1[i-1]+k1*h/2)
                k3 = func(x1[i-1]+h/2,y1[i-1]+k2*h/2)
                k4 = func(x1[i-1]+h,y1[i-1]+h*k3)
                y1[i] = y1[i-1] + (k1+2*k2+2*k3+k4)*h/6
        print(x1,y1)
        return x1,y1
    def Draw():
    #三条线画在同一张图中 发现误差很小
        xEuler,yEuler = ImproEuler()
        xKutta,yKutta = RungeKutta()
        plt.plot(xEuler,yEuler,'*',xKutta,yKutta,'gx',xEuler,np.sqrt(2*xEuler+1),'r')
        plt.annotate(r'Euler Formula or Runge Kutta',xy=(xEuler[5],yEuler[5]),xytext=(xEuler[5],yEuler[5]-0.2),arrowprops=dict(facecolor='black',shrink=0.1,width=2))
        plt.annotate(r'Correct results',xy=(xEuler[10]-0.05,np.sqrt(2*xEuler[10]+1)),xytext=(xEuler[10]-0.05,np.sqrt(2*xEuler[10]+1)-0.2),arrowprops=dict(facecolor='black',shrink=0.1,width=2))
        plt.grid(True)
        plt.savefig('ODE test',dpi = 600)
        plt.show()
    if __name__ == '__main__':
        Draw()
    更多相关内容
  • 常微分方程初值问题数值解法.ppt
  • 常微分方程初值问题数值解法的比较.pdf
  • 常微分方程初值问题数值解法matlab程序及说明
  • 计算方法第九章常微分方程初值问题数值解法.ppt
  • (1)常微分方程初值问题数值解法 (2)解题步骤 (3)数值微分解法 (4)数值积分解法 2、所有知识点代码 3、结果---以三阶Runge-Kutta公式为例(其他的类示) 1、概述 (1)常微分方程初值问题数值解法 ...

    目录

    1、概述

    (1)常微分方程初值问题数值解法 

    (2)解题步骤

    (3)数值微分解法

     (4)数值积分解法

    2、所有知识点代码 

    3、结果---以三阶Runge-Kutta公式为例(其他的类似)


    1、概述

    (1)常微分方程初值问题数值解法 

    (2)解题步骤

                         

    (3)数值微分解法

     (4)数值积分解法

    2、所有知识点代码 

    import numpy as np
    import matplotlib.pyplot as plt
    
    def funEval(x,y):             #近似值
        fxy = (x*y-y**2)/x**2 #1
        #fxy=2*y/x+x**2*np.e**x #2
        #stablity
        #fxy=-30*y
        #fxy= np.e**(-x**2)
        #fxy=1-y
        #fxy=2*y/x+x**2*np.e**x
        return fxy
    def funtrue(x):            #真实值
        ft=x/(0.5+np.log(x)) #1
        #ft=x**2*(np.e**x-np.e) #2
        #stablityu
        #ft = np.e**(-30*x) 
        #ft = 1-np.e**(-x) 
        #ft=x**2*(np.e**x-np.e)
        return ft
    def Euler(a,b,f,y0,n):      #Euler公式
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] =y0
        x[0] = a
        for i in range(1,n,1):
            x[i]=a+i*h
            y[i]= y[i-1]+h*f(x[i-1],y[i-1])
        return x,y
    def ModEuler(a,b,f,y0,n):     #改进Euler公式
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
    
        y[0] =y0
        x[0] = a
        for i in range(1,n,1):
            x[i]=a+i*h
            y[i]= y[i-1]+h*f(x[i-1],y[i-1])
            y[i] = y[i-1]+h/2*(f(x[i-1],y[i-1])+f(x[i],y[i]))
        return x,y
    
    def Heun(a,b,f,y0,n):       #二阶Runge—Kutta方法:Heun公式
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] =y0
        x[0] = a
        K1,K2=0,0
        for i in range(1,n,1):
            x[i]=a+i*h
            K1 = f(x[i-1],y[i-1])
            K2 = f(x[i-1]+2/3*h,y[i-1]+2/3*h*K1)
            y[i] = y[i-1]+h/4*(K1+3*K2)
        return x,y
    
    def Ord3Kutta(a,b,f,y0,n):      #三阶Kutta公式
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] =y0
        x[0] = a
        K1,K2,K3=0,0,0
        for i in range(1,n,1):
            x[i]=a+i*h
            K1 = f(x[i-1],y[i-1])
            K2 = f(x[i-1]+1/2*h,y[i-1]+1/2*h*K1)
            K3 = f(x[i-1]+h,y[i-1]-h*K1+2*h*K2)
            y[i] = y[i-1]+h/6*(K1+4*K2+K3)
        return x,y
    
    def Ord3Heun(a,b,f,y0,n):     #三阶Heun公式
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] =y0
        x[0] = a
        K1,K2,K3=0,0,0
        for i in range(1,n,1):
            x[i]=a+i*h
            K1 = f(x[i-1],y[i-1])
            K2 = f(x[i-1]+1/3*h,y[i-1]+1/3*h*K1)
            K3 = f(x[i-1]+2/3*h,y[i-1]+2/3*h*K2)
            y[i] = y[i-1]+h/4*(K1+3*K2)
        return x,y
    
    def Ord4Kutta(a,b,f,y0,n):        #四阶古典Runge—Kutta公式
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] = y0
        x[0] = a
        K1,K2,K3,K4 = 0,0,0,0
        for i in range(1,n,1):
            x[i]=a+i*h
            K1 = f(x[i-1],y[i-1])
            K2 = f(x[i-1]+1/2*h,y[i-1]+1/2*h*K1)
            K3 = f(x[i-1]+1/2*h,y[i-1]+1/2*h*K2)
            K4 = f(x[i-1]+h,y[i-1]+h*K3)
            y[i] = y[i-1]+h/6*(K1+2*K2+2*K3+K4)
        return x,y
    
    def Ord4Kutta2(a,b,f,y0,n):        #四阶Kutta公式
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] = y0
        x[0] = a
        K1,K2,K3,K4 = 0,0,0,0
        for i in range(1,n,1):
            x[i]=a+i*h
            K1 = f(x[i-1],y[i-1])
            K2 = f(x[i-1]+1/3*h,y[i-1]+1/3*h*K1)
            K3 = f(x[i-1]+2/3*h,y[i-1]-1/3*h*K1+h*K2)
            K4 = f(x[i-1]+h,y[i-1]+h*K1-h*K2+h*K3)
            y[i] = y[i-1]+h/8*(K1+3*K2+3*K3+K4)
        return x,y
    
    
    def ExpAdamsK1(a,b,f,y0,n):    #二步显式Adams
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] = y0
        x[0] = a
        for i in range(1,n,1):
            x[i]=a+i*h
            if i ==1:
                y[i] =y[i-1]+h*f(x[i-1],y[i-1])
            else:
                y[i] = y[i-1]+h/2*(3*f(x[i-1],y[i-1])-f(x[i-2],y[i-2])) 
        return x,y
    
    def ExpAdamsK2(a,b,f,y0,n):     #三步显式Adams
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] = y0
        x[0] = a
        if n>=2:
            for i in range(1,n,1):
                x[i]=a+i*h
                if i <=2:
                    y[i] =y[i-1]+h*f(x[i-1],y[i-1])
                else:
                    y[i] = y[i-1]+h/12*(23*f(x[i-1],y[i-1])-16*f(x[i-2],y[i-2])+5*f(x[i-3],y[i-3])) 
        else:
            print("n must be larger than 2")
        return x,y
    
    def ModifiedAdamsK2(a,b,f,y0,n):   #预测-校正法:四阶Adams预测-校正法
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        
        y[0] = y0
        x[0] = a
        temp=0
    
        if n>4:
            x1,y1=Ord4Kutta2(a,h*4,f,y0,4)
            y[0:4]=y1
            for i in range(4,n,1):
                x[i]=a+i*h
                temp = y[i]+h/24*(55*f(x[i-1],y[i-1])-59*f(x[i-2],y[i-2])+37*f(x[i-31],y[i-3])-9*f(x[i-4],y[i-4]))
                y[i] = y[i-1]+h/24*(9*temp+19*f(x[i-1],y[i-1])-5*f(x[i-2],y[i-2])+f(x[i-3],y[i-3])) 
        else:
            print("n must be larger than 4")
    
        return x,y
    ##y'=-30y,y(0)=1, 0<=x<=0.6
    def ImplictEuler(a,b,n):    #显式Euler公式和隐式Euler法精确度的比较
        h=np.abs(b-a)/(n-1)
        y=np.zeros((n,1))
        x=np.zeros((n,1))
        y0 = 1
        y[0] = y0
        x[0] = a
        K1,K2,K3,K4 = 0,0,0,0
        for i in range(1,n,1):
            x[i]=a+i*h
            y[i] = y[i-1]/(1+30*h)
        yt = np.e**(-30*x)
        return x,y,yt
    def main():
    
        a,b = 1,1.5 #12
     
        n = [2,5,10,20,40,80,160,320,640]
        #n=[5]
        y0 = 2 #1
        #y0 = 0 #2
    
        #stablity
        #a,b = 1,2
        #y0 = 0
        # for i in n:
        #     x,y=Euler(a,b,funEval,y0,i)
        #     plt.plot(x,y,label='Euler-solution-'+str(i))
        #     plt.plot(x,funtrue(x),label='True-solution-'+str(i))
        #     plt.legend()
        #     plt.show() 
    
        for i in n:
            #x,y=ModEuler(a,b,funEval,y0,i)
            #x,y=ModifiedAdamsK2(a,b,funEval,y0,i)
            #x,y,yt=ImplictEuler(0,0.6,i)
            x,y=Ord3Kutta(a, b, funEval, y0, i)
            yt= funtrue(x)
            plt.plot(x,y,label='Euler-solution-'+str(i))
            plt.plot(x,yt,label='True-solution-'+str(i))
            plt.legend()
            plt.show()
            print(x[0:5],(y-yt)[0:5])
            print(y)
    
     
    
    
    
    
    if __name__ == '__main__':
        main()

    3、结果---以三阶Runge-Kutta公式为例(其他的类似)

    [[1. ]
     [1.5]] [[ 0.        ]
     [-0.03207385]]
    [[2.        ]
     [1.62453333]]

     

    展开全文
  • 数值分析 常微分方程初值问题数值解法PPT学习教案.pptx
  • 一、实验平台要求不限,程序语言采用基本高级语言(注:推荐使用C/C++,根据课上要 二、实验报告撰写格式:1)实验要求(实验题目和初始数据),2)算法描述(文字
  • 常微分方程初值问题数值解法PPT教案.pptx
  • 常微分方程初值问题数值解法PPT学习教案.pptx
  • 计算方法常微分方程初值问题数值解法PPT课件.pptx
  • 计算方法常微分方程初值问题数值解法PPT学习教案.pptx
  • 计算方法第六章常微分方程初值问题数值解法PPT课件.pptx
  • 计算方法第六章常微分方程初值问题数值解法PPT学习教案.pptx
  • 计算方法 常微分方程初值问题数值解法Euler公式龙格库塔法PPT课件.pptx
  • 常微分方程初值问题数值解法问题一、一阶常微分方程初值问题的有限差分方法与误差分析二、向前Euler法及误差分析1.向前Euler法2.误差分析3.后退Euler法三、改进欧拉公式四、龙格—库塔方法 问题 一阶常微分方程的...

    问题

    一阶常微分方程的初值问题:
    y ′ = f ( x , y ) y^{'}=f(x,y) y=f(x,y)
    y ( x 0 ) = y 0 y(x_0)=y_0 y(x0)=y0
    定理一
    设f在区域 D = ( x , y ) ∣ a ≤ x ≤ b , y ∈ R D={(x,y)|a\leq x\leq b,y\in R} D=(x,y)axb,yR上连续,关于y满足 ∣ f ( x , y 1 ) − f ( x , y 2 ) ∣ ≤ L ∣ y 1 − y 2 ∣ , ∀ x 0 ∈ [ a , b ] , y 0 ∈ R |f(x,y_1)-f(x,y_2)|\leq L|y_1-y_2|,\forall x_0\in[a,b],y_0\in R f(x,y1)f(x,y2)Ly1y2,x0[a,b],y0R(这条件叫做Lipschitz条件)常微分方程的初值问题存在唯一的连续可微解 y ( x ) y(x) y(x)

    一、一阶常微分方程初值问题的有限差分方法与误差分析

    只有一些特殊的常微分方程能求解析解,实际问题中主要靠数值解,所谓数值解法,就是寻求解y(x)在一系列离散节点 x 0 < x 1 < x 2 . . . . . . x_0<x_1<x_2...... x0<x1<x2......上的近似值 y 0 , y 1 , y 2 . . . . . . y_0,y_1,y_2...... y0,y1,y2......
    整体截断误差
    g i = ∣ y i − y ( x i ) ∣ g_i=|y_i-y(x_i)| gi=yiy(xi)
    其中 y ( x i ) y(x_i) y(xi)就是我们的精确解, y i y_i yi就是给定计算格式(迭代公式,初值是精确值)之后所计算出来的数值解。
    局部截断误差
    e i = ∣ y ( x i + 1 ) − w i + 1 ∣ e_i=|y(x_{i+1})-w_{i+1}| ei=y(xi+1)wi+1
    其中 w i + 1 w_{i+1} wi+1初值为 y ( x i ) y(x_i) y(xi) 时迭代出的数值解
    对于某种方法,局部截断误差满足 e i = O ( h p + 1 ) e_i=O(h^{p+1}) ei=O(hp+1),则称该方法具有p阶精度。

    二、向前Euler法及误差分析

    1.向前Euler法

    先设步长相等,大小为h,对常微分方程 y ′ = f ( x , y ) y^{'}=f(x,y) y=f(x,y)中的倒数用均差近似,即
    y ( x i + 1 ) − y ( x i ) h ≈ y ′ = f ( x i , y ( x i ) ) \frac{y(x_{i+1})-y(x_i)}{h}\approx y^{'}=f(x_i,y(x_i)) hy(xi+1)y(xi)y=f(xi,y(xi))
    若初值已知,利用迭代格式
    { y 0 = y ( x 0 ) y i + 1 = y i + h f ( x i , y i ) \begin{cases} y_0=y(x_0)\\ y_{i+1}=y_i+hf(x_i,y_i) \end{cases} {y0=y(x0)yi+1=yi+hf(xi,yi)
    可逐次算出 y 1 , y 2 . . . . . y_1,y_2..... y1,y2.....

    2.误差分析

    对于精确值做泰勒展开式
    y ( x i + 1 ) = y ( x i ) + h y ′ ( x i ) + h 2 2 y ′ ′ ( c ) y(x_{i+1})=y(x_i)+hy^{'}(x_i)+\frac{h^2}{2}y^{''}(c) y(xi+1)=y(xi)+hy(xi)+2h2y(c)其中 c ∈ ( x i , x i + 1 ) c \in(x_i,x_{i+1}) c(xi,xi+1),其中 y ′ ( x i ) = f ( x i , y ( x i ) ) y^{'}(x_i)=f(x_i,y(x_i)) y(xi)=f(xi,y(xi))
    对于欧拉法的迭代格式,由局部截断误差定义可得
    y i + 1 = y ( x i ) + h y ′ ( x i ) y_{i+1}=y(x_i)+hy^{'}(x_i) yi+1=y(xi)+hy(xi)
    局部截断误差 e i = h 2 2 y ′ ′ ( c ) e_i=\frac{h^2}{2}y^{''}(c) ei=2h2y(c)

    3.后退Euler法

    通过利用均差近似倒数 y ′ ( x i + 1 ) y^{'}(x_{i+1}) y(xi+1)

    y ( x i + 1 ) − y ( x i ) h ≈ y ′ ( x i + 1 ) = f ( x i + 1 , y ( x i + 1 ) ) \frac{y(x_{i+1})-y(x_i)}{h}\approx y^{'}(x_{i+1})=f(x_{i+1},y(x_{i+1})) hy(xi+1)y(xi)y(xi+1)=f(xi+1,y(xi+1))
    因此后退欧拉法的递推公式为
    y i + 1 = y i + h f ( x i + 1 , y i + 1 ) y_{i+1}=y_i+hf(x_{i+1},y_{i+1}) yi+1=yi+hf(xi+1,yi+1)(1式)

    相比欧拉法,后退欧拉法的递推公式并不能直接进行计算,而需要解方程,即此公式是隐式的,而欧拉法的递推公式是可以直接计算的,即为显式的。隐式方程通常用迭代法求解,而迭代过程的实质是逐步隐式化。
    { y i + 1 0 = y i + h f ( x i , y i ) y i + 1 1 = y i + h f ( x i , y i + 1 0 ) . . . . . . y i + 1 k + 1 = y i + h f ( x i , y i + 1 k ) ( 2 式 ) \begin{cases} y_{i+1}^{0}=y_i+hf(x_i,y_i)\\ y_{i+1}^{1}=y_i+hf(x_i,y_{i+1}^{0}) \\ ......\\ y_{i+1}^{k+1}=y_i+hf(x_i,y_{i+1}^{k}) (2式) \end{cases} yi+10=yi+hf(xi,yi)yi+11=yi+hf(xi,yi+10)......yi+1k+1=yi+hf(xi,yi+1k)2
    由于f(x,y)满足Lipschitz条件,利用(1)式减去(2)式得到
    ∣ y i + 1 k + 1 − y i + 1 ∣ = h ( f ( x i , y i + 1 k ) − h f ( x i + 1 , y i + 1 ) ) ≤ h L ∣ y i + 1 k − y i + 1 ∣ |y_{i+1}^{k+1}-y_{i+1}|=h(f(x_i,y_{i+1}^{k})-hf(x_{i+1},y_{i+1}))\leq hL|y_{i+1}^{k}-y_{i+1}| yi+1k+1yi+1=h(f(xi,yi+1k)hf(xi+1,yi+1))hLyi+1kyi+1
    h L < 1 hL<1 hL<1则迭代法收敛到解 y i + 1 y_{i+1} yi+1

    三、改进欧拉公式

    补充:梯形方法
    y i + 1 = y i + h 2 ( f ( x i + 1 , y i + 1 ) + f ( x i , y i ) ) y_{i+1}=y_i+\frac{h}{2}(f(x_{i+1},y_{i+1})+f(x_{i},y_{i})) yi+1=yi+2h(f(xi+1,yi+1)+f(xi,yi))

    先用欧拉法求得一个初步的近似值 y ‾ i + 1 \overline{y}_{i+1} yi+1,称为预报值,然后用它替代梯形法右端的 y i + 1 y_{i+1} yi+1,得到左端校正值 y i + 1 {y}_{i+1} yi+1,可以校正n次,但是改进欧拉法中,我们只校正一次,这样建立的预报-校正系统称为改进的欧拉格式:
    { y ‾ i + 1 = y i + h f ( x i , y i ) y i + 1 = y i + h 2 ( f ( x i + 1 , y ‾ i + 1 ) + f ( x i , y i ) ) \begin{cases} \overline{y}_{i+1}=y_i+hf(x_i,y_i)\\ y_{i+1}=y_i+\frac{h}{2}(f(x_{i+1},\overline{y}_{i+1})+f(x_{i},y_{i})) \end{cases} {yi+1=yi+hf(xi,yi)yi+1=yi+2h(f(xi+1,yi+1)+f(xi,yi))
    也可以写成下面一个形式:
    { y p = y n + h f ( x n , y n ) y c = y n + h f ( x n , y p ) y n + 1 = 1 2 ( y p + y c ) \begin{cases} {y}_{p}=y_n+hf(x_n,y_n)\\ {y}_{c}=y_n+hf(x_n,y_p)\\ y_{n+1}=\frac{1}{2}(y_p+y_c) \end{cases} yp=yn+hf(xn,yn)yc=yn+hf(xn,yp)yn+1=21(yp+yc)

    四、单步法局部截断误差与阶

    初值问题 y ′ = f ( x , y ) y^{'}=f(x,y) y=f(x,y) y ( x 0 ) = y 0 y(x_0)=y_0 y(x0)=y0的单步法可用一般形式表示为
    y n + 1 = y n + h φ ( x n , y n , y n + 1 , h ) y_{n+1}=y_{n}+h\varphi(x_n,y_n,y_{n+1},h) yn+1=yn+hφ(xnynyn+1h)
    其中 φ ( x n , y n , y n + 1 , h ) \varphi(x_n,y_n,y_{n+1},h) φ(xnynyn+1h)称为增量函数
    显式单步法的表示形式
    y n + 1 = y n + h φ ( x n , y n , h ) y_{n+1}=y_{n}+h\varphi(x_n,y_n,h) yn+1=yn+hφ(xnynh)
    例如:对欧拉法 φ ( x n , y n , h ) = f ( x , y ) \varphi(x_n,y_n,h)=f(x,y) φ(xnynh)=f(x,y)
    对于显式单步法的局部截断误差
    T n + 1 = y ( x n + 1 ) − y ( x n ) − h φ ( x n , y ( x n ) , h ) T_{n+1}=y(x_{n+1})-y(x_{n})-h\varphi(x_n,y(x_n),h) Tn+1=y(xn+1)y(xn)hφ(xny(xn)h)

    五、龙格—库塔方法

    1.定义

    对于方程 y ′ = f ( x , y ) y^{'}=f(x,y) y=f(x,y)的等价积分形式

    y ( x n + 1 ) − y ( x n ) = ∫ x n x n + 1 f ( x , y ( x ) )   d x y(x_{n+1})-y(x_{n})=\int_{x_{n}}^{x_{n+1}} {f(x,y(x))} \ dx y(xn+1)y(xn)=xnxn+1f(x,y(x)) dx

    方程右端可用求积公式表示为
    ∫ x n x n + 1 f ( x , y ( x ) )   d x ≈ h ∑ i = 1 r c i f ( x n + λ i h , y ( x n + λ i h ) ) \int_{x_{n}}^{x_{n+1}} {f(x,y(x))} \ dx \approx h\sum_{i=1}^rc_if(x_n+\lambda_ih,y(x_n+\lambda_ih)) xnxn+1f(x,y(x)) dxhi=1rcif(xn+λih,y(xn+λih))
    一般来说点数r越多精度越高,上式右端相当于增量函数 φ ( x n , y n , h ) \varphi(x_n,y_n,h) φ(xnynh),为了得到便于计算的显示方法,将公式表示为
    y n + 1 = y n + h φ ( x n , y n , h ) y_{n+1}=y_n+h\varphi(x_n,y_n,h) yn+1=yn+hφ(xn,yn,h)
    其中
    φ ( x n , y n , h ) = ∑ i = 1 n c i K i \varphi(x_n,y_n,h)=\sum_{i=1}^n c_iK_i φ(xn,yn,h)=i=1nciKi
    K 1 = f ( x n , y n ) K_1=f(x_n,y_n) K1=f(xn,yn)
    K 2 = f ( x n + λ 2 h , y n + h u 21 K 1 ) , λ 2 = h K_2=f(x_n+\lambda_2h,y_n+hu_{21}K_1),\lambda_2=h K2=f(xn+λ2h,yn+hu21K1)λ2=h

    K i = f ( x n + λ i h , y n + h × ∑ j = 1 i − 1 u i j K j ) , λ i = ∑ j = 1 i − 1 u i j K_i=f(x_n+\lambda_ih,y_n+h\times \sum_{j=1}^{i-1} u_{ij}K_j),\lambda_i=\sum_{j=1}^{i-1} u_{ij} Ki=f(xn+λih,yn+h×j=1i1uijKj)λi=j=1i1uij
    将此方法称为龙格库r级塔显示方法

    2.常用的龙格库塔方法

    常用二阶显式龙格库塔方法
    y n + 1 = y n + h K 2 y_{n+1}=y_n+hK_2 yn+1=yn+hK2
    K 1 = f ( x n , y n ) K_1=f(x_n,y_n) K1=f(xn,yn)
    K 2 = f ( x n + h 2 , y n + h 2 K 1 ) K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1) K2=f(xn+2h,yn+2hK1)
    常用三阶
    y n + 1 = y n + h 6 ( K 1 + c 2 K 2 ) y_{n+1}=y_n+\frac{h}{6}(K_1+c2K_2) yn+1=yn+6h(K1+c2K2)
    K 1 = f ( x n , y n ) K_1=f(x_n,y_n) K1=f(xn,yn)
    K 2 = f ( x n + h 2 , y n + h 2 K 1 ) K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1) K2=f(xn+2h,yn+2hK1)
    K 3 = f ( x n + h , y n − h K 1 + 2 h K 2 ) K_3=f(x_n+h,y_n-hK_1+2hK_2) K3=f(xn+h,ynhK1+2hK2)
    常用四阶(记下来,要考,哈哈哈,matlab计算内核所在)
    y n + 1 = y n + h 6 ( K 1 + 2 K 2 + 2 K 3 + K 4 ) y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4) yn+1=yn+6h(K1+2K2+2K3+K4)
    K 1 = f ( x n , y n ) K_1=f(x_n,y_n) K1=f(xn,yn)
    K 2 = f ( x n + h 2 , y n + h 2 K 1 ) K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1) K2=f(xn+2h,yn+2hK1)
    K 3 = f ( x n + h 2 , y n + h 2 K 2 ) K_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2) K3=f(xn+2h,yn+2hK2)
    K 4 = f ( x n + h , y n + h K 3 ) K_4=f(x_n+h,y_n+hK_3) K4=f(xn+h,yn+hK3)

    六、单步法的收敛性与稳定性

    1.收敛性与相容性

    定义:若一种数值方法(如单步法)对于固定的 x n = x 0 + n h x_n=x_0+nh xn=x0+nh,当 h → 0 h\rightarrow 0 h0 ,有 y n = y ( x n ) y_n=y(x_n) yn=y(xn),则称该方法收敛。
    定理:假设单步法具有p阶精度,且增量函数 φ ( x , y , h ) \varphi(x,y,h) φ(xyh)满足利普希茨条件
    ∣ φ ( x , y , h ) − φ ( x , y ‾ , h ) ∣ ≤ L φ ∣ y − y ‾ ∣ |\varphi(x,y,h)-\varphi(x,\overline{y},h)|\leq L_{\varphi}|y-\overline{y}| φ(xyh)φ(xyh)Lφyy
    又设初值准确,则其整体截断误差
    y ( x n ) − y n = O ( h p ) y(x_n)-y_n=O(h^p) y(xn)yn=O(hp)
    定义:若单步法 y n + 1 = y n + h φ ( x n , y n , h ) y_{n+1}=y_n+h{\varphi(x_n,y_n,h)} yn+1=yn+hφ(xn,yn,h)的增量函数满足 φ ( x , y , 0 ) = f ( x , y ) \varphi(x,y,0)=f(x,y) φ(x,y,0)=f(x,y),则称单步法与初值问题相容(前述的单步法是用了泰勒展开式略去高阶项近似,有用求积公式离散得到的数值方法,相容性指数值方法逼近微分方程,在 h → 0 h\rightarrow0 h0时,可得到 y ′ ( x ) = f ( x , y ) y^{'}(x)=f(x,y) y(x)=f(x,y)

    2.绝对稳定性与绝对稳定域

    定义:若一种数值方法在节点值 y n y_n yn上大小扰动为 δ \delta δ,于以后的节点值 y m ( m > n ) y_m(m>n) ym(m>n)上产生的偏差均不超过 δ \delta δ,则称该方法是稳定的
    定义:对于模型方程 y ′ = λ y y^{'}=\lambda y y=λy,使用单步法解模型方程 y n + 1 = E ( λ h ) y n y_{n+1}=E(\lambda h)y_n yn+1=E(λh)yn,若得到 ∣ E ( h λ ) ∣ < 1 |E(h\lambda )|<1 E(hλ)<1,则称该方法是绝对稳定的(以此可以来确定步长)。
    常用:对于四阶龙格库塔方法: − 2.78 < h λ < 0 -2.78<h\lambda<0 2.78<hλ<0,即 0 < h < − 2.78 / λ 0<h<-2.78/\lambda 0<h<2.78/λ

    展开全文
  • 计算方法 常微分方程初值问题数值解法Euler公式龙格库塔法PPT学习教案.pptx
  • 此书是计算数学领域中常微分方程初值问题的最基本的书,全英文的教材,对研究生的深入学习起到很好的作用。
  • 初值问题:即满足初值条件的常微分方程的解 首先,常微分方程得有解---有解条件----利普希茨条件--- 两类问题: ①单步法 ---计算下一个点的值只需要用到前面一个点的值 ②多步法---计算下一个点的值需要...
    • 初值问题:即满足初值条件的常微分方程的解

    首先,常微分 方程得有解-- -有解条件----利普希茨条件---

     

    两类问题:

    ①单步法 ---计算下一个点的值  只需要用到前面一个点的值 

    ②多步法---计算下一个点的值  需要用到前面  个点的值 

    1、欧拉法---下一个点的计算值等于前一个点的计算值加上步长乘以前一个点的函数值

    ①初值条件:给出因变量在某个点上的值
    • 具体过程
    第一步:将 初值条件带入微分方程,得到在该点的导数值
    第二步:在该点,用taylor进行展开, 舍去二次项,将一次函数近似函数y
    第三:计算 第二个点在直线上的值,用这个值近似函数  的在第二个点的值,依此类推,直到迭代完成

    一些批注:显式欧拉方程指下一步要计算的值,不在迭代方程中;隐式欧拉方程指下一步要计算的值,在迭代方程中。

    • 显示欧拉法与隐式欧拉法
    评价: ① 精度为1(整体) ②误差阶为2(局部)
    • 梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均,构造的计算公式.

    • 改进的欧拉方法---

    思想:因为梯形公式是隐式公式,将显式欧拉公式对下一步的计算值进行预估,用梯形公式对下一步的计算值进行校正.

    评价: ① 精度为2(整体) ②误差阶为3(局部)

    2、龙格-库塔方法

     

    • 2阶-龙格-库塔方法----类似改进的欧拉法

     

    • 4阶-龙格-库塔方法

    k1为区间左端点的斜率;

    k2为区间中点的斜率,通过欧拉法,用斜率k1确定函数f在中点上的斜率k2

    k3也是区间中点的斜率,通过欧拉法,用斜率k2确定函数f在中点上的斜率k3

    k4是区间右端点的斜率,通过欧拉法,用斜率k3确定函数f在端点上的斜率k4

    4阶-龙格-库塔方法---- 局部误差  ;  整体误差为 

     

     

    展开全文
  • 实验报告七常微分方程初值问题数值解法.pdf实验报告七常微分方程初值问题数值解法.pdf实验报告七常微分方程初值问题数值解法.pdf实验报告七常微分方程初值问题数值解法.pdf实验报告七常微分方程初值问题的...
  • 实验报告七 常微分方程初值问题数值解法.pdf实验报告七 常微分方程初值问题数值解法.pdf实验报告七 常微分方程初值问题数值解法.pdf实验报告七 常微分方程初值问题数值解法.pdf实验报告七 常微分方程初值...
  • 实验报告七常微分方程初值问题数值解法.docx实验报告七常微分方程初值问题数值解法.docx实验报告七常微分方程初值问题数值解法.docx实验报告七常微分方程初值问题数值解法.docx实验报告七常微分方程初值问题...
  • 实验报告七 常微分方程初值问题数值解法.docx实验报告七 常微分方程初值问题数值解法.docx实验报告七 常微分方程初值问题数值解法.docx实验报告七 常微分方程初值问题数值解法.docx实验报告七 常微分方程...
  • import numpy as np import matplotlib.pyplot as plt x = 0 y = 1 h = 0.1 # x(n+1) - x(n) X = np.zeros(100) Y = np.zeros(100) for i in range(0,100): y = 1.1 * y - 0.2 * (x / y) ...0.2 1
  • 注重算法与程序实现,强调理论知识与程序设计的紧密结合,既有理论性,也有实用性,对每个...配有大量图形,侧重从几何含义的角度直观地说明问题;最后一章是与所学内容紧密结合的上机实验与指导;附录有部分习题答案。
  • 本文主要介绍了解常微分方程初值问题的欧拉方法、龙格-库塔法、一阶方程组和高阶方程的数值解法、边值问题数值解法等内容
  • 第六章 常微分方程初值问题数值解法微分方程是模拟自然、生物、工程等现象的重要数学工具。包含自变量、未知函数的导数或微分的方程称为微分方程。在求解微分方程时,须

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 753
精华内容 301
热门标签
关键字:

常微分方程初值问题数值解法