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

一、问题提出

科学计算中经常遇到微分方程（组）初值问题，需要利用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隐藏了许多细节部分，不适合对方法的深入理解，但是写代码速度快，代码更简练。

posted @ 2011-12-02 18:58 涵曦阅读(7995) 评论(3) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

博客园

涵曦~博客园

专注技术 - 技术乐园

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

《 Skynet 游戏服务器开发实战》

公告