精华内容
下载资源
问答
  • (二)要求 用编程语言实现用改进的欧拉(Euler)公式求解常微分方程初值问题、用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题的程序。 二、示例 1、问题用改进的欧拉(Euler)公式求解常微分方程初值问题...

    一、目的与要求
    (一)目的通过设计、编制、调试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()
    展开全文
  • 本文提出了一些基于Legendre多项式作为基函数的k步线性多步方法的... 利用导出的方法对常微分方程的一些数值例子进行了求解,以证明其有效性和准确性。 获得的数值结果表明,所提出的方法也可以有效地解决此类问题
  • 常微分方程初值问题数值解法 一、问题提出 科学计算中经常遇到微分方程(组)初值问题,需要利用Euler法,改进Euler法,Rung-Kutta方法求其数值解,诸如以下问题: (1) (3) (4)利用四阶标准R- K...

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

    一、问题提出

    科学计算中经常遇到微分方程(组)初值问题,需要利用Euler法,改进Euler法,Rung-Kutta方法求其数值解,诸如以下问题:

    (1)

    (3)

     

    (4)利用四阶标准R- K方法求二阶方程初值问题的数值解

    二、问题分析

        使用Euler法求解,运算程序简单,但是计算结果准确度不高。使用改进的Euler法求解过程相对复杂,但是准确度会更高。准确度最高的是四阶龙格库塔法,求解步骤也是最复杂的。问题(1)使用Euler求解,并与准确解对比。问题(3)使用改进的Euler法求解。问题(4)(I)(IV)使用四届标准龙格库塔法求解。

    三、实验程序及注释

    1.Euler法:

    View Code
    //differential.cpp
    #include <iostream>
    using namespace std;
    #include <math.h>

    float h, //步长
    x00, //x的取值范围
    xn, // x0< x < xn
    y00; //初值
    float x[21];//结果
    float y[21];

    void Euler()
    {
    init();
    x[0] = x00;
    y[0] = y00;
    for (int i=1; i<=(xn-x00)/h+1; i++)
    {
    x[i] = x[i-1] + h;
    y[i] = y[i-1] + h * f(x[i-1],y[i-1]);
    }
    }



    2.解问题(1):

    View Code
    //常微分方程y’= f(xn,yn)
    float f(float x, float y)
    {
    return 4*x/y-x*y;
    //return y-2*x/y;//课本上的例题,由于检测程序的正确性
    }
    //初始化初值和x的范围
    void init()
    {
    x00=0;
    xn=2;
    y00=3;
    }
    //求常微分方程的标准解
    float pf(float x)
    {
    return sqrt(4+5*exp(-x*x));
    }
    int main()
    {
    h=0.1;
    Euler();
    cout << "步长 h = "
    << h << endl;
    for (int i=1; i<=(xn-x00)/h+1; i++)
    {
    cout << "x" << i << "=" << x[i] << "\t"
    << "y" << i << "=" << y[i] << "\t"
    << "y("<< x[i] << ")=" << pf(x[i]) << endl;
    }
    h=0.2;
    Euler();
    cout << "步长 h = "
    << h << endl;
    for (i=1; i<=(xn-x00)/h+1; i++)
    {
    cout << "x" << i << "=" << x[i] << "\t"
    << "y" << i << "=" << y[i] << "\t"
    << "y("<< x[i] << ")=" << pf(x[i]) << endl;
    }
    h=0.4;
    Euler();
    cout << "步长 h = "
    << h << endl;
    for (i=1; i<=(xn-x00)/h+1; i++)
    {
    cout << "x" << i << "=" << x[i] << "\t"
    << "y" << i << "=" << y[i] << "\t"
    << "y("<< x[i] << ")=" << pf(x[i]) << endl;
    }
    return 0;
    }



    3.改进的Euler法:

    View Code
    #include <iostream>
    #include <iomanip>
    using namespace std;

    float h, //步长
    x00, //x的取值范围
    xn, // x0< x < xn
    y00; //初值
    float y10,y20,y30;//初值

    float x[21];//结果
    float y1[21],y2[21],y3[21];
    //改进的Euler法
    void newEuler()
    {
    init();
    x[0] = x00;
    y1[0] = y10;
    y2[0] = y20;
    y3[0] = y30;

    for (int i=1; i<=(xn-x00)/h+1; i++)
    {
    x[i] = x[i-1] + h;
    float yp=y1[i-1]+h*f1(x[i-1],y1[i-1],y2[i-1],y3[i-1]);
    float yc=y1[i-1]+h*f1(x[i],yp,y2[i-1],y3[i-1]);
    y1[i] = 0.5*(yp+yc);
    yp=y2[i-1]+h*f2(x[i-1],y1[i-1],y2[i-1],y3[i-1]);
    yc=y2[i-1]+h*f2(x[i],y1[i-1],yp,y3[i-1]);
    y2[i] = 0.5*(yp+yc);
    yp=y3[i-1]+h*f3(x[i-1],y1[i-1],y2[i-1],y3[i-1]);
    yc=y3[i-1]+h*f3(x[i],y1[i-1],y2[i-1],yp);
    y3[i] = 0.5*(yp+yc);
    }
    }



    4.解问题(3):

    View Code
    float f1(float x, float y1, float y2, float y3)
    {
    return y2;
    }
    float f2(float x, float y1, float y2, float y3)
    {
    return -y1;
    }
    float f3(float x, float y1, float y2, float y3)
    {
    return -y3;
    }
    void init()
    {
    x00=0;
    xn=0.15;
    y10=-1;
    y20=0;
    y30=1;
    }
    int main()
    {
    h=0.01;
    newEuler();
    cout << "步长 h = "
    << h << endl;
    cout.precision(6);
    for (int i=0; i<=(xn-x00)/h; i+=5)
    {
    cout << "x(" << setw(2) << i << ")=" << x[i] << ""
    << "y1(" << setw(2) << x[i] << ")=" << y1[i] << ""
    << "y2(" << setw(2) << x[i] << ")=" << y2[i] << ""
    << "y3(" << setw(2) << x[i] << ")=" << y3[i] << endl;
    }
    return 0;
    }



    5.龙格库塔法:

    View Code
    #include <iostream>
    #include <iomanip>
    using namespace std;

    float h, //步长
    x00, //x的取值范围
    xn; // x0< x < xn
    float y10;
    float y20;//初值
    float x[51];//结果
    float y1[51],y2[51];
    //龙格库塔法
    void L_k()
    {
    float L[5];
    init();
    x[0]=x00;
    y1[0]=y10;
    y2[0]=y20;
    for (int i=1; i<=(xn-x00)/h; i++)
    {
    x[i] = x[i-1] + h;
    L[1] = f2(x[i-1],y1[i-1],y2[i-1]);
    L[2] = f2(x[i-1]+h/2,y1[i-1]+h*y2[i-1]/2,y2[i-1]+h*L[1]/2);
    L[3] = f2(x[i-1]+h/2,y1[i-1]+h*y2[i-1]/2+h*h*L[1]/4,y2[i-1]+h*L[2]/2);
    L[4] = f2(x[i-1]+h,y1[i-1]+h*y2[i-1]+h*h*L[2]/2,y2[i-1]+h*L[3]);
    y1[i] = y1[i-1] + h*y2[i-1] + h*h*(L[1]+L[2]+L[3])/6;
    y2[i] = y2[i-1] + h*(L[1]+2*L[2]+2*L[3]+L[4])/6;
    }
    }



    6.解问题(4)(I):

    View Code
    //  y1'=y2
    // y2'=3*y2-*y1
    float f1(float x, float y1, float y2)
    {
    return y2;
    }
    float f2(float x, float y1, float y2)
    {
    return 3*y2-2*y1;
    }

    void init()
    {
    x00=0;
    xn=1;
    y10=0;
    y20=1;
    }

    int main()
    {
    h=0.02;
    L_k();
    cout << "步长 h = "
    << h << endl;
    cout.precision(6);
    for (int i=0; i<=(xn-x00)/h; i++)
    {
    cout << "x(" << setw(2) << i << ")=" << x[i] << ""
    << "y1(" << setw(2) << x[i] << ")=" << y1[i] << ""
    << "y2(" << setw(2) << x[i] << ")=" << y2[i] << endl;
    }
    return 0;
    }



    7.解问题(4)(IV):

    View Code
    //  y1'=y2
    // y2'=-sin(y1)
    float f1(float x, float y1, float y2)
    {
    return y2;
    }
    float f2(float x, float y1, float y2)
    {
    return -sin(y1);
    }

    void init()
    {
    x00=0;
    xn=4;
    y10=1;
    y20=0;
    }
    int main()
    {
    h=0.2;
    L_k();
    ……
    ……
    return 0;
    }



    四、实验数据结果及分析

    1.解问题(1):

         

    Euler法

    步长 h = 0.1

    x1=0.1 y1=3 y(0.1)=2.9917

    x2=0.2 y2=2.98333 y(0.2)=2.96714

    x3=0.3 y3=2.95048 y(0.3)=2.9274

    x4=0.4 y4=2.90264 y(0.4)=2.87415

    x5=0.5 y5=2.84166 y(0.5)=2.80963

    x6=0.6 y6=2.76995 y(0.6)=2.73649

    x7=0.7 y7=2.6904 y(0.7)=2.65766

    x8=0.8 y8=2.60615 y(0.8)=2.57613

    x9=0.9 y9=2.52044 y(0.9)=2.49485

    x10=1 y10=2.43643 y(1)=2.41648

    x11=1.1 y11=2.35697 y(1.1)=2.34329

    x12=1.2 y12=2.28438 y(1.2)=2.27698

    x13=1.3 y13=2.22038 y(1.3)=2.21869

    x14=1.4 y14=2.16592 y(1.4)=2.16894

    x15=1.5 y15=2.12124 y(1.5)=2.12767

    x16=1.6 y16=2.08591 y(1.6)=2.0944

    x17=1.7 y17=2.05898 y(1.7)=2.0683

    x18=1.8 y18=2.03922 y(1.8)=2.04837

    x19=1.9 y19=2.02523 y(1.9)=2.03353

    x20=2 y20=2.01571 y(2)=2.02276

    步长 h = 0.2

    x1=0.2 y1=3 y(0.2)=2.96714

    x2=0.4 y2=2.93333 y(0.4)=2.87415

    x3=0.6 y3=2.80776 y(0.6)=2.73649

    x4=0.8 y4=2.64178 y(0.8)=2.57613

    x5=1 y5=2.46136 y(1)=2.41648

    x6=1.2 y6=2.29411 y(1.2)=2.27698

    x7=1.4 y7=2.16199 y(1.4)=2.16894

    x8=1.6 y8=2.07467 y(1.6)=2.0944

    x9=1.8 y9=2.02774 y(1.8)=2.04837

    x10=2 y10=2.0079 y(2)=2.02276

    步长 h = 0.4

    x1=0.4 y1=3 y(0.4)=2.87415

    x2=0.8 y2=2.73333 y(0.8)=2.57613

    x3=1.2 y3=2.32696 y(1.2)=2.27698

    x4=1.6 y4=2.03513 y(1.6)=2.0944

    x5=2 y5=1.99055 y(2)=2.02276

     

        2.解问题(3):

    步长 h = 0.01

    x( 0)=0 y1( 0)=-1 y2( 0)=0 y3( 0)=1

    x( 5)=0.05 y1(0.05)=-0.999 y2(0.05)=0.04999 y3(0.05)=0.95123

    x(10)=0.1 y1(0.1)=-0.995502 y2(0.1)=0.09988 y3(0.1)=0.904839

    x(15)=0.15 y1(0.15)=-0.989514 y2(0.15)=0.149545 y3(0.15)=0.86071

    参考结果:

    对比看,结果相差不大,参考结果应该是使用四阶龙格库塔方法求解的。

    3.解问题(4)(I):

    步长 h = 0.02

    x( 0)=0 y1( 0)=0 y2( 0)=1

    x( 1)=0.02 y1(0.02)=0.0206094 y2(0.02)=1.06142

    x( 2)=0.04 y1(0.04)=0.0424763 y2(0.04)=1.12576

    x( 3)=0.06 y1(0.06)=0.0656603 y2(0.06)=1.19316

    x( 4)=0.08 y1(0.08)=0.0902238 y2(0.08)=1.26373

    x( 5)=0.1 y1(0.1)=0.116232 y2(0.1)=1.33763

    x( 6)=0.12 y1(0.12)=0.143752 y2(0.12)=1.415

    x( 7)=0.14 y1(0.14)=0.172856 y2(0.14)=1.49599

    x( 8)=0.16 y1(0.16)=0.203617 y2(0.16)=1.58074

    x( 9)=0.18 y1(0.18)=0.236112 y2(0.18)=1.66944

    x(10)=0.2 y1(0.2)=0.270422 y2(0.2)=1.76225

    x(11)=0.22 y1(0.22)=0.30663 y2(0.22)=1.85934

    x(12)=0.24 y1(0.24)=0.344825 y2(0.24)=1.9609

    x(13)=0.26 y1(0.26)=0.385097 y2(0.26)=2.06712

    x(14)=0.28 y1(0.28)=0.427543 y2(0.28)=2.17821

    x(15)=0.3 y1(0.3)=0.47226 y2(0.3)=2.29438

    x(16)=0.32 y1(0.32)=0.519353 y2(0.32)=2.41583

    x(17)=0.34 y1(0.34)=0.56893 y2(0.34)=2.54281

    x(18)=0.36 y1(0.36)=0.621104 y2(0.36)=2.67554

    x(19)=0.38 y1(0.38)=0.675991 y2(0.38)=2.81427

    x(20)=0.4 y1(0.4)=0.733716 y2(0.4)=2.95926

    x(21)=0.42 y1(0.42)=0.794405 y2(0.42)=3.11077

    x(22)=0.44 y1(0.44)=0.858192 y2(0.44)=3.26909

    x(23)=0.46 y1(0.46)=0.925216 y2(0.46)=3.43451

    x(24)=0.48 y1(0.48)=0.995622 y2(0.48)=3.60732

    x(25)=0.5 y1(0.5)=1.06956 y2(0.5)=3.78784

    x(26)=0.52 y1(0.52)=1.14719 y2(0.52)=3.97641

    x(27)=0.54 y1(0.54)=1.22867 y2(0.54)=4.17335

    x(28)=0.56 y1(0.56)=1.31418 y2(0.56)=4.37903

    x(29)=0.58 y1(0.58)=1.40389 y2(0.58)=4.59383

    x(30)=0.6 y1(0.6)=1.498 y2(0.6)=4.81811

    x(31)=0.62 y1(0.62)=1.59668 y2(0.62)=5.0523

    x(32)=0.64 y1(0.64)=1.70016 y2(0.64)=5.2968

    x(33)=0.66 y1(0.66)=1.80863 y2(0.66)=5.55205

    x(34)=0.68 y1(0.68)=1.92232 y2(0.68)=5.81851

    x(35)=0.7 y1(0.7)=2.04145 y2(0.7)=6.09665

    x(36)=0.72 y1(0.72)=2.16626 y2(0.72)=6.38696

    x(37)=0.74 y1(0.74)=2.29701 y2(0.74)=6.68996

    x(38)=0.76 y1(0.76)=2.43395 y2(0.76)=7.00617

    x(39)=0.78 y1(0.78)=2.57735 y2(0.78)=7.33617

    x(40)=0.8 y1(0.8)=2.72749 y2(0.8)=7.68052

    x(41)=0.82 y1(0.82)=2.88467 y2(0.82)=8.03984

    x(42)=0.84 y1(0.84)=3.04919 y2(0.84)=8.41474

    x(43)=0.86 y1(0.86)=3.22137 y2(0.86)=8.8059

    x(44)=0.88 y1(0.88)=3.40154 y2(0.88)=9.21397

    x(45)=0.9 y1(0.9)=3.59004 y2(0.9)=9.63969

    x(46)=0.92 y1(0.92)=3.78725 y2(0.92)=10.0838

    x(47)=0.94 y1(0.94)=3.99352 y2(0.94)=10.547

    x(48)=0.96 y1(0.96)=4.20926 y2(0.96)=11.0302

    x(49)=0.98 y1(0.98)=4.43487 y2(0.98)=11.5342

    x(50)=1 y1( 1)=4.67077 y2( 1)=12.0598



    3.解问题(4)(IV):

    步长 h = 0.2

    x( 0)=0 y1( 0)=1 y2( 0)=0

    x( 1)=0.2 y1(0.2)=0.983201 y2(0.2)=-0.167682

    x( 2)=0.4 y1(0.4)=0.933176 y2(0.4)=-0.331607

    x( 3)=0.6 y1(0.6)=0.851087 y2(0.6)=-0.48757

    x( 4)=0.8 y1(0.8)=0.739008 y2(0.8)=-0.630607

    x( 5)=1 y1( 1)=0.60009 y2( 1)=-0.754957

    x( 6)=1.2 y1(1.2)=0.438688 y2(1.2)=-0.854407

    x( 7)=1.4 y1(1.4)=0.260389 y2(1.4)=-0.923023

    x( 8)=1.6 y1(1.6)=0.0718545 y2(1.6)=-0.956155

    x( 9)=1.8 y1(1.8)=-0.119534 y2(1.8)=-0.951378

    x(10)=2 y1( 2)=-0.306182 y2( 2)=-0.90905

    x(11)=2.2 y1(2.2)=-0.480844 y2(2.2)=-0.832224

    x(12)=2.4 y1(2.4)=-0.637103 y2(2.4)=-0.725968

    x(13)=2.6 y1(2.6)=-0.769673 y2(2.6)=-0.596373

    x(14)=2.8 y1(2.8)=-0.874507 y2(2.8)=-0.449594

    x(15)=3 y1( 3)=-0.948738 y2( 3)=-0.291209

    x(16)=3.2 y1(3.2)=-0.990535 y2(3.2)=-0.125976

    x(17)=3.4 y1(3.4)=-0.998943 y2(3.4)=0.0420493

    x(18)=3.6 y1(3.6)=-0.973777 y2(3.6)=0.209154

    x(19)=3.8 y1(3.8)=-0.915597 y2(3.8)=0.371509

    x(20)=4 y1( 4)=-0.825779 y2( 4)=0.524739



    五、实验结论

       使用c语言写代码虽然比Matlab的代码长,但是还是c语言适合用来做练习,因为用c写更有结构。Matlab隐藏了许多细节部分,不适合对方法的深入理解,但是写代码速度快,代码更简练。

    转载于:https://www.cnblogs.com/hanxi/archive/2011/12/02/2272597.html

    展开全文
  • 初值问题:即满足初值条件的常微分方程的解 首先,常微分方程得有解---有解条件----利普希茨条件---定理1 (利普希茨条件)若存在正数L,使得对任意 ,有 定理2 (解存在性) ①若函数 在方区域 连续,②函数 关于y...
    • 初值问题:即满足初值条件的常微分方程的解

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

    定理1 (利普希茨条件)若存在正数L,使得对任意

    ,有

    定理2 (解存在性) ①若函数

    在方区域
    连续,②函数
    关于y满足利普希茨条件,

    则对任意

    常微分方程存在唯一连续可微数值解.
    由于原函数是无法精确求解出来的,我们只需要求解在某点的值就好

    两类问题:

    ①单步法 ---计算下一个点的值

    只需要用到前面一个点的值

    ②多步法---计算下一个点的值

    需要用到前面
    个点的值

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

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

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

    怎么计算隐式欧拉方程----要借助显示欧拉迭代计算---一般用迭代法

    -----迭代---将微分方程在区间

    进行积分,然后函数f进行近似,即可得到迭代方程-----迭代方程收敛性?由函数关于y满足利普希茨条件,可以推出迭代公式收敛。
    • 局部截断误差:假设前n步误差为0,我们计算第n+1步的误差,将次误差称为局部截断误差,且局部误差为
    • p阶精度:由理论证明若局部误差阶的时间复杂度为
      ,则整体误差阶为
      我们称公式精度为p。
    • 显示欧拉法与隐式欧拉法
    评价: ① 精度为1(整体) ②误差阶为2(局部)
    • 梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均,构造的计算公式.
    • 改进的欧拉方法---

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

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

    2、龙格-库塔方法

    思想:根据Lagrange中值定理,下一次的计算值可以用前一次的计算值加上h乘以前一个点的斜率;而这个斜率用该区间上的多个点的斜率的算数平均来逼近

    注意:怎么计算任意斜率

    ?第
    个点的斜率
    有微分方程可以算出

    所以要算的

    值,由欧拉法即可算出,
    • 2阶-龙格-库塔方法----类似改进的欧拉法

    根据Lagrange中值定理,下一次的计算值可以用前一次的计算值加上h乘以斜率;而这个斜率用区间上的端点中点的斜率的算数平均来逼近

    其中,

    2阶龙格-库塔方法--局部误差为
    (误差阶为3);整体误差
    (精度为2)
    • 4阶-龙格-库塔方法

    42c8d73178c734fe159702218f07cfb5.png

    k1为区间左端点的斜率;

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

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

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

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

    2.1 单步法收敛性与稳定性

    • 收敛性: 对任意固定的点,函数y在该点的真实值,随着迭代步数增加,所用公式的计算值都会趋于所产生的近似值都趋于真实值,称该方法收敛
    公式收敛性判别:函数 f 对y满足利普希茨条件(充分条件),则数值微分求解公式是收敛的、
    相容:欧拉法的增量函数与函数 f 相等,则称 欧拉法与初值问题是相容的
    • 稳定性:微分方程数值法在某点的计算值发生细微的扰动
      ,之后的计算值所
      产生的误差不超过
      ,
      则称该方法稳定
    绝对稳定:稳定的微分数值解法中,当某一步计算值有误差时,以后计算值的误差不会扩大,则称这种稳定是绝对稳定。绝对稳定域:就是小区间的步长,使得微分数值公式是绝对稳定,而这些步长构成的区间称为绝对稳定域.

    3、线性多步法

    思想:下一次的计算值等于
    个的计算值
    的算法平均加上步长乘以多个斜率的算数平均
    • 阿当姆斯外插法
    思想:下一次的计算值等于前一次的计算值加上过前
    个的计算值的
    插值函数在区间
    上的积分

    • 阿当姆斯内插法
    思想:下一次的计算值等于前一次的计算值加上过前
    的计算值的插值函数在区间上的积分

    两者区别:下一次计算值在不在迭代公式中....

    展开全文
  • 1、用改进的欧拉(Euler)公式求解常微分方程初值问题。 2、用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题。 欧拉公式求解y’=y-2x/y和y(0)=1在[0,1]上的结果 sy3_1.m clear; clc; x0=0; y0=1; h=0.1; N=...

    实验目的

    编程实现以下算法:
    1、用改进的欧拉(Euler)公式求解常微分方程初值问题。
    2、用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题。

    欧拉公式求解y’=y-2x/y和y(0)=1在[0,1]上的结果
    sy3_1.m

    clear;
    clc;
    x0=0;
    y0=1;
    h=0.1;
    N=10;
    n=1;
    X=[];
    Y=[];
    while n<=N
        x1=x0+h;
        yp=y0+h*(y0-2*x0/y0);
        yc=y0+h*(yp-2*x1/yp);
        y1=(yp+yc)/2;
        X=[X x1];
        Y=[Y y1];
        x0=x1;
        y0=y1;
        n=n+1;
    end
    X
    Y
    

    标准结果:
    标准结果欧拉公式求解结果:
    欧拉公式

    四阶龙格-库塔方法求解y’=y-2x/y和y(0)=1在[0,1]上的结果
    sy3_2.m

    clear;
    clc;
    x0=0;
    y0=1;
    h=0.1;
    N=10;
    n=1;
    X=[];
    Y=[];
    while n<=N;
        x1=x0+h;
        K1=y0-2*x0/y0;
        K2=y0+h*K1/2-(2*x0+h)/(y0+h*K1/2);
        K3=y0+h*K2/2-(2*x0+h)/(y0+h*K2/2);
        K4=y0+h*K3-2*(x0+h)/(y0+h*K3);
        y1=y0+(K1+2*K2+2*K3+K4)*h/6;
        X=[X x1];
        Y=[Y y1];
        x0=x1;
        y0=y1;
        n=n+1;
    end
    X
    Y
    

    四阶龙格-库塔方法求解结果:
    四阶龙格-库塔方法

    展开全文
  • 常微分方程组数值求解;主要内容 数值求解常微分方程组函数概述 非刚性/刚性常微分方程问题求解 隐式微分方程组求解 微分代数方程(DAE)与延迟微分方程(DDE)求解 边值问题求解 ;...第二节 非刚性/刚性常微分方程初值问题
  • 问题:用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题。 算法描述 算法的程序框图: 具体算法: (1)读取a,b,n,f (2)计算步长h = (b-a)/n,x0=a,y0=f (3)从i = 1开始每次i++,重复以下操作,直到i=n...
  • 我们的重点是开发和实施新的两步混合方法,用于直接求解一般的二阶常微分方程。 在该方法的开发中,采用幂级数作为基础函数,并且将出现的方程组微分系统并置在所有电网和离网点。 所得方程在选定点处插值。 然后,...
  • //用RKG法求解微分方程 #include #include #include using namespace std; class rkg { private:  int i, n;  double a, b, c, d, f, h, k1, k2, k3, k4, x, xf, y; public:  double func(double z, double t)
  • \qquad常微分方程是描述连续变化的...考虑一阶常微分方程初值问题 y′=f(x,y),x∈[x0,b],(1)y(x0)=y0(2) y^{&#x27;}=f(x,y),\quad x\in [x_{0},b],\quad(1)\\ y(x_{0})=y_{0}\quad(2) y′=f(x,y),x∈[x0​,...
  • 实验五、常微分方程初值问题的数值解法 ​ 常微分方程的求解问题在实践中经常遇到,因此研究常微分方程的数值解法就很有必要。欧拉方法是最简单、最基本的方法,利用差商代替微商,就可得到一系列欧拉公式。这些...
  • //用Runge_Kutta_Fehlberg法求解微分方程 #include #include #include #include using namespace std; class rkf { private:  int flag;  double eps, error, f, h, hnew, x, xf, y, yold;  double k1, k2, k
  • y0为初值。 %返回yn、xn向量,并绘制曲线 yn=zeros(1,((b-a)/h)+1); yn(1)=y0; xn=a:h:b; for i=1:((b-a)/h) k1=h*f(xn(i),yn(i)); k2=h*f(xn(i)+h/2,yn(i)+k1/2); k3=h*f(xn(i)+h/2,yn(i)+k2/2); k4=h*f(xn
  • (1)利用Euler方法和改进的Euler方法求解初值问题 (2)利用Runge-Kutta方法求解初值问题 步骤: import math def f1(x, y): if x == 0: return 0 else: return ((4 * x) / y) - (x * y) def f2(x, y): ...
  • //用RKG法求解常微分方程组 #include #include #include #include using namespace std; class s_ode { private:  int i, j, l, n;  double a, b, c, d, h, x, x1, xi, xf;  double *b1, *b2, *b3, *g, *y, *...
  • //实现Milne预报-校正法 #include #include #include using namespace std; class milne { private:  int i, n;  double f, h, x, x_last, y, y1, y2, y3, yc, yp;  double f0, f1, f2;... double fun
  • 我们采用了Duan-Rach-Wazwaz改进的Adomian分解方法来数值求解非线性常微分方程组的初值问题。 为了确认该方法的实用性,鲁棒性和可靠性,在某些情况下,我们将改进的Adomian分解方法与MATHEMATICA解决方案以及四阶...
  • //实现多步Euler法 #include #include #include using namespace std; class multi_euler { private:  int i, n;... double f, h, x, x_last, yc, yc_old, yp;... double func(double z, double t) ...
  • 求解常微分方程初值问题就是寻找函数y(t)使之满足如下方程: y’=f(t,y),to≤t≤b y(to)=yo所谓其数值解法,就是求y(t)在离散结点t处的函数近似值y。的方法,ya~y(xn)。 这些近似值称为常微分方程初值...
  • //实现Euler预报_校正法 #include #include #include using namespace std; class euler { private:  int i, n;... double f, h, x, x_last, yc, yp;... double func(double z, double t) ... f = (1
  • 采用欧拉法、改进欧拉法、龙格库塔法(经典RK法)求解常微分方程初值问题的自编MATLAB代码。所有函数均独立成文件便于移植。代码的使用结合一个具体题目说明,题目来源为浙江大学数值计算方法作业。
  • 结果如图,当[x1,y1]=euler(f,0,1/3,1,20)时,函数值很大,求解错误。 ![图片说明](https://img-ask.csdn.net/upload/202006/13/1592032393_549155.jpg) 当[x1,y1]=euler(f,0,1/3,0.3,20)时,近似正确 问题原题...
  • 常微分方程数值求解常微分方程数值求解的一般概念常微分方程数值求解的一般概念常微分方程...求解常微分方程初值问题就是寻找函数y(t)使之满足如下方程: y’=f(t,y),t0≤t≤b ,y(t0)=y0 所谓其数值解法,就是求...
  • 常微分方程的数值求解

    千次阅读 2018-04-13 11:59:32
    常微分方程 首先理解一下什么是常微分方程,简单的说就是只有一个未知数的微分方程,具体定义如下: 凡含有参数,未知函数和未知...一阶常微分方程初值问题是: {y′=f(x,y),y(x0)=y0,{y′=f(x,y),y(x0)=y0, \beg...
  • 本节将介绍常微分方程初值问题的数值求解,主要内容分为三个部分:常微分方程数值求解的概念、求解函数及刚性问题。 一、常微分方程数值求解的一般概念 首先,凡含有参数,未知函数和未知函数导数 (或微分) 的方程,...
  • import numpy as np from matplotlib import pyplot as plt plt.rcParams['font.sans-serif']=['SimHei'] plt.rcParams['axes.unicode_minus'] = False def myfun1(x, y1, y2): return y2 ...

空空如也

空空如也

1 2 3 4 5 ... 8
收藏数 147
精华内容 58
关键字:

常微分方程初值问题求解