精华内容
下载资源
问答
  • 常微分方程初值问题数值解法.ppt
  • 常微分方程初值问题数值解法 一、问题提出 科学计算中经常遇到微分方程(组)初值问题,需要利用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隐藏了许多细节部分,不适合对方法的深入理解,但是写代码速度快,代码更简练。

    展开全文
  • 常微分方程初值问题数值解法问题一、一阶常微分方程初值问题的有限差分方法与误差分析二、向前Euler法及误差分析1.向前Euler法2.误差分析3.后退Euler法三、改进欧拉公式四、龙格—库塔方法 问题 一阶常微分方程的...

    问题

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

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

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

    二、向前Euler法及误差分析

    1.向前Euler法

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

    2.误差分析

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

    3.后退Euler法

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

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

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

    三、改进欧拉公式

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

    先用欧拉法求得一个初步的近似值yi+1\overline{y}_{i+1},称为预报值,然后用它替代梯形法右端的yi+1y_{i+1},得到左端校正值yi+1{y}_{i+1},可以校正n次,但是改进欧拉法中,我们只校正一次,这样建立的预报-校正系统称为改进的欧拉格式:
    {yi+1=yi+hf(xi,yi)yi+1=yi+h2(f(xi+1,yi+1)+f(xi,yi)) \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}
    也可以写成下面一个形式:
    {yp=yn+hf(xn,yn)yc=yn+hf(xn,yp)yn+1=12(yp+yc) \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}

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

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

    五、龙格—库塔方法

    1.定义

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

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

    方程右端可用求积公式表示为
    xnxn+1f(x,y(x)) dxhi=1rcif(xn+λih,y(xn+λih))\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))
    一般来说点数r越多精度越高,上式右端相当于增量函数φ(xnynh)\varphi(x_n,y_n,h),为了得到便于计算的显示方法,将公式表示为
    yn+1=yn+hφ(xn,yn,h)y_{n+1}=y_n+h\varphi(x_n,y_n,h)
    其中
    φ(xn,yn,h)=i=1nciKi\varphi(x_n,y_n,h)=\sum_{i=1}^n c_iK_i
    K1=f(xn,yn)K_1=f(x_n,y_n)
    K2=f(xn+λ2h,yn+hu21K1)λ2=hK_2=f(x_n+\lambda_2h,y_n+hu_{21}K_1),\lambda_2=h

    Ki=f(xn+λih,yn+h×j=1i1uijKj)λi=j=1i1uijK_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}
    将此方法称为龙格库r级塔显示方法

    2.常用的龙格库塔方法

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

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

    1.收敛性与相容性

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

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

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

    展开全文
  • 第9章-常微分方程初值问题数值解法9.1 引言9.2 简单的数值方法9.3 龙格-库塔方法9.4 单步法的收敛性与稳定性9.5 线性多步法9.6 线性多步法的收敛性与稳定性9.7 一阶方程组与刚性方程组 9.1 引言 9.2 简单的数值方法 ...
    第8章 回到目录


    9.1 引言

    一阶常微分方程
    y=f(x,y),x[x0,b](1.1)y'=f(x,y),\quad x\in[x_0,b] \tag{1.1}

    y(x0)=y0(1.2)y(x_0)=y_0 \tag{1.2}


    利普希茨 (Lipschitz) 条件 / 利普希茨常数

    定理1 解的存在唯一性定理

    ff 在区域 D={(x,y)axb,yR}D=\{(x,y)|a\le x \le b, y\in\mathbb{R}\} 上连续,关于 yy 满足利普希茨条件,则对任意 x0[a,b],y0Rx_0\in[a,b], y_0\in\mathbb{R},常微分方程初值问题(1.1)式和(1.2)式当 x[a,b]x \in[a,b] 时存在唯一的连续可微解 y(x)y(x)


    定理2 解对初值依赖的敏感性

    9.2 简单的数值方法

    积分曲线上一点 (x,y)(x,y) 的切线斜率等于函数 f(x,y)f(x,y) 的值。


    欧拉 (Euler) 方法

    y(xn+1)y(xn)hy(xn)=f(xn,y(xn))\frac{y(x_{n+1})-y(x_n)}{h} \approx y'(x_n) = f(x_n, y(x_n))


    9.3 龙格-库塔方法


    9.4 单步法的收敛性与稳定性


    9.5 线性多步法


    9.6 线性多步法的收敛性与稳定性


    9.7 一阶方程组与刚性方程组

    展开全文
  • \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​,...

    \qquad常微分方程是描述连续变化的数学语言,微分方程的求解就是确定给定方程的可微方程y(t)y(t)。考虑一阶常微分方程的初值问题
    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) \qquad如果存在实数L&gt;0L&gt;0,使得
    f(x,y1)f(x,y2)Ly1y2,y1,y2R, |f(x,y_{1})-f(x,y_{2})|\leq L|y_{1}-y_{2}|,\forall y_{1},y_{2}\in\mathbb{R}, 则称ff关于yy满足利普西茨(LipschiteLipschite)条件LL称为ff的利普西茨常数(简称Lisp.Lisp.常数)
    \qquad定理1ff在区域D={(x,y)axb,yR}D=\{(x,y)|a\leq x\leq b,y\in\mathbb{R}\}上连续,关于yy满足利普西茨条件,则对任意x0[a,b]y0Rx_{0}\in [a,b],y_{0}\in\mathbb{R},常微分方程初值问题满足(1)(2)式当x[a,b]x\in [a,b]时存在唯一的连续可微解y(x)y(x)
    \qquad所谓数值解法,就是寻求解y(x)y(x)在一系列离散节点
    x1&lt;x2&lt;&lt;xn&lt; x_{1}&lt;x_{2}&lt;\cdots&lt;x_{n}&lt;\cdots \quad上的近似值y1,y2,&ThinSpace;,yn,y_{1},y_{2},\cdots,y_{n},\cdots相邻两节点的间距hn=xn+1xnh_{n}=x_{n+1}-x_{n}称为步长。求解过程首先对常微分方程离散化,建立求数值解的递推公式。当计算yn+1y_{n+1}时只用到前一点yny_{n},称为单步法。用到前面kk个点的值yn,yn1,y_{n},y_{n-1},\cdots时称为kk步法
    \qquad欧拉法 设做出折线(由切线相连(x0,x1,xnx_{0},x_{1},\cdots x_{n}处作切线))的顶点PnP_{n},过Pn(xn,yn)P_{n}(x_{n},y_{n})依方向场的方向推进到Pn+1(xn+1,yn+1)P_{n+1}(x_{n+1},y_{n+1}),显然两个顶点Pn,Pn+1P_{n},P_{n+1}的坐标有关系
    yn+1ynxn+1xn=f(xn,yn) \frac{y_{n+1}-y_{n}}{x_{n+1}-x_{n}}=f(x_{n},y_{n}) \qquad对微分方程式(1)中的导数用均差近似,即
    yn+1ynhy(xn)=f(xn,y(xn)) \frac{y_{n+1}-y_{n}}{h}\approx y^{&#x27;}(x_{n})=f(x_{n},y(x_{n})) \qquad均差近似逐次计算。为分析计算公式的精度,使用泰勒展开将yn+1y_{n+1}xnx_{n}处展开,则有
    y(xn+1)=y(xn+h)=y(xn)+y(xn)h+h22y(ξn),ξn(xn,xn+1) y(x_{n+1})=y(x_{n}+h)=y(x_{n})+y^{&#x27;}(x_{n})h+\frac{h^{2}}{2}y^{&#x27;&#x27;}(\xi_{n}),\quad \xi_{n}\in(x_{n},x_{n+1}) \qquadyn=yxny_{n}=y_{x_{n}}的前提下,f(xn,yn)=f(xn,y(xn))=y(xn)f(x_{n},y_{n})=f(x_{n},y(x_{n}))=y^{&#x27;}(x_{n}),于是欧拉法的误差
    y(xn+1)yn+1=h22y(ξn)h22y(xn) y(x_{n+1})-y_{n+1}=\frac{h^{2}}{2}y^{&#x27;&#x27;}(\xi_{n})\approx\frac{h^{2}}{2}y^{&#x27;&#x27;}(x_{n}) 称为此方法的局部截断误差。显式单步法的局部截断误差一般性的表示方法为:
    Tn+1=y(xn+1)y(xn)hφ(xn,y(xn),h)=O(hp+1) T_{n+1}=y(x_{n+1})-y(x_{n})-h\varphi(x_{n},y(x_{n}),h)=O(h^{p+1}) 则称该方法具有pp阶精度。
    \qquad如果对微分方程从xnx_{n}xn+1x_{n+1}积分,得
    y(xn+1)=y(xn)+xnxn+1f(t,y(t))dt y(x_{n+1})=y(x_{n})+\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt 右端积分使用右矩形公式hf(xn+1,y(xn+1))hf(x_{n+1},y(x_{n+1}))近似得到另一个公式
    yn+1=yn+hf(xn+1,yn+1) y_{n+1}=y_{n}+hf(x_{n+1},y_{n+1}) 通常使用迭代法求解(实现xnxn+1x_{n}\to x_{n+1}的运算),设用欧拉公式
    yn+1(0)=yn+hf(xn,yn) y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n}) 给出迭代初值yn+1(0)y_{n+1}^{(0)}带入上式得
    yn+1(1)=yn+hf(xn+1,yn+1(0)) y_{n+1}^{(1)}=y_{n}+hf(x_{n+1},y^{(0)}_{n+1}) 反复迭代得到
    yn+1(k+1)=yn+hf(xn+1,yn+1(k)),k=0,1, y_{n+1}^{(k+1)}=y_{n}+hf(x_{n+1},y^{(k)}_{n+1}),k=0,1,\cdots 由利普西茨条件得(\color{#F00}{证明迭代法收敛})
    yn+1(k+1)yn+1=hf(xn+1,yn+1(k))f(xn+1,yn+1)hLyn+1(k)yn+1 |y^{(k+1)}_{n+1}-y_{n+1}|=h|f(x_{n+1},y^{(k)}_{n+1})-f(x_{n+1},y^{}_{n+1})|\leq hL|y^{(k)}_{n+1}-y_{n+1}| 由此,只要hL&lt;1hL&lt;1使用迭代法就能收敛到解yn+1y_{n+1}
    \qquad梯形方法(增加精度) 公式为yn+1=yn+h2[f(xn,yn)+f(xn+1,yn+1)]y_{n+1} = y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},y_{n+1})],迭代公式为:
    yn+1(0)=yn+hf(xn,yn)yn+1(k+1)=yn+h2(f(xn,yn)+f(xn+1,yn+1(k)),k=0,1,2, y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n})\\ y_{n+1}^{(k+1)}=y_{n}+\frac{h}{2}(f(x_{n},y_{n})+f(x_{n+1},y_{n+1}^{(k)}),\quad k=0,1,2,\cdots \qquad改进欧拉公式(增加了欧拉公式的精度,减小梯形算法迭代带来的计算次数),建立预测-校正系统。

    预测 yˉn+1=yn+hf(xn,yn)\bar{y}_{n+1}=y_{n}+hf(x_{n},y_{n})

    校正 yn+1=yn+h2[f(xn,yn)+f(xn+1,yˉn+1)]y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},\bar{y}_{n+1})]
    \qquad龙格-库塔法,由公式
    y(xn+1)=y(xn)+xnxn+1f(t,y(t))dt y(x_{n+1})=y(x_{n})+\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt 要使近似精度提高就是提高右式的近似精度,近似公式为:
    xnxn+1f(t,y(t))dthi=1rcif(xn+λih,y(xn+λih)) \int_{x_{n}}^{x_{n+1}}f(t,y(t))dt\approx h\sum_{i=1}^{r}c_{i}f(x_{n}+\lambda _{i}h,y(x_{n}+\lambda _{i}h)) 分割的段越多,rr越大,上式右端(相当于增量函数φ(xn,y(xn)h)\varphi(x_{n},y(x_{n})h))的近似精度越高。由公式
    yn+1=yn+hφ(xn,yn,h) y_{n+1} = y_{n} + h\varphi(x_{n},y_{n},h) 其中,
    φ(xn,yn,h)=i=1rciKiK1=f(xn,yn)Ki=f(xn+λih,yn+j=1i1μijKj)i=2,3,&ThinSpace;,r \varphi(x_{n},y_{n},h)=\sum_{i=1}^{r}c_{i}K_{i}\\ K_{1}=f(x_{n},y_{n})\\ K_{i}=f(x_{n}+\lambda_{i}h,y_{n}+\sum_{j=1}^{i-1}\mu_{ij}K_{j})\quad i=2,3,\cdots,r 通过给定局部截断误差精度求解未知系数。对于步长的选择,由给定的精度ε\varepsilon和偏差Δ=yn+1(h2)yn+1(h)\Delta=|y_{n+1}^{(\frac{h}{2})}-y_{n+1}^{(h)}|(\color{#F00}{折半前后两次的偏差})决定。当精度大于偏差时,减小步长(对折),反之则反。
    \qquad单步法的收敛性 假设单步法具有pp阶精度,且增量函数φ(x,y,h)\varphi(x,y,h)关于yy满足利普西茨条件
    φ(x,y,h)φ(x,yˉ,h)Lφyyˉ |\varphi(x,y,h)-\varphi(x,\bar{y},h)|\leq L_{\varphi}|y-\bar{y}| 又设初值y0y_{0}准确的,即y0=y(x0)y_{0}=y(x_{0}),其整体的截断误差
    y(xn)yn=O(hp) y(x_{n})-y_{n}=O(h^{p}) 稳定性 由模型方程y=λyy^{&#x27;}=\lambda y的欧拉公式为:
    yn+1=(1+hλ)yn y_{n+1}=(1+h\lambda)y_{n} 设在节点值yny_{n}上有一个扰动值εn\varepsilon_{n},它的传播使节点yn+1y_{n+1}产生εn+1\varepsilon_{n+1}的扰动值,扰动值满足
    εn+1=(1+hλ)εn \varepsilon_{n+1}=(1+h\lambda)\varepsilon_{n} 为使差分方程的解非增(解稳定),选择hh充分小,使得1+hλ1|1+h\lambda|\leq1,即2&lt;λ&lt;0-2&lt;\lambda&lt;0

    展开全文
  • 此书是计算数学领域中常微分方程初值问题的最基本的书,全英文的教材,对研究生的深入学习起到很好的作用。
  • (一)目的通过设计、编制、调试1~2个求常微分方程初值问题数值解解的程序,加深对其数值计算方法及有关的基础理论知识的理解。 (二)要求 用编程语言实现用改进的欧拉(Euler)公式求解常微分方程初值问题、用四...
  • 常微分方程初值问题数值解法 本文参考书为马东升著《数值计算方法》 引言 未知函数为一元函数的微分方程叫常微分方程,讨论一阶常微分方程的初值问题 {y′=f(x,y)y(x0)=y0 \begin{cases}y'=f(x,y)\\y(x_0)=y_0\...
  • 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=...
  • 注重算法与程序实现,强调理论知识与程序设计的紧密结合,既有理论性,也有实用性,对每个...配有大量图形,侧重从几何含义的角度直观地说明问题;最后一章是与所学内容紧密结合的上机实验与指导;附录有部分习题答案。
  • 实验五、常微分方程初值问题数值解法 ​ 常微分方程的求解问题在实践中经常遇到,因此研究常微分方程的数值解法就很有必要。欧拉方法是最简单、最基本的方法,利用差商代替微商,就可得到一系列欧拉公式。这些...
  • 数值分析(数值计算) 数学建模 实验报告 matlab程序
  • 本文主要介绍了解常微分方程初值问题的欧拉方法、龙格-库塔法、一阶方程组和高阶方程的数值解法、边值问题数值解法等内容
  • 科技信息SCIENCE & TECHNOLOGY INFORMATION2012 年 第 7 期 0 引言 在自然科学和经济生活的许多领域中,常常会遇到一阶常微分方程的初值问题如下:y'(t)=f(t,y)y(t0 )=y0 0 , t0≤t (1... 1 常微分方程初值问题...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 522
精华内容 208
关键字:

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