欧拉法求解已知初值微分方程解

1.原理

数值积分算法是求解知初值的微分方程的重要方法。

如已知微分方程

d（y)/d(t) = f(y, t)

y(t0) = y0

方程两边对t积分就会有

 

此式表示原函数t1时刻的解y(t1)为原函数初值y(t0)与微分方程表达式（即f(y, t)）的初值时刻t0到t1时刻的积分值之和。而积分的含义表示积分段微分表达式曲线与坐标轴围成的面积。即原函数t1时刻的解y(t1)为原函数初值y(t0)与微分方程表达式的初值时刻t0到t1时刻与坐标轴围成的面积之和。原函数值与微分方程表达式可表示为如下形式

 

图 原函数值与微分方程表达式关系

 

2.欧拉法

  如果现在采用计算机来求解微分方程的解，由于计算机只能计算离散值，所以只能采用计算机计算的离散值来无限的逼近连续值。如果让t0与t1之间的间隔趋近于无限小，那么途中红色部分面积（t0到t1积分）就可以近似的计算：s =f0 * (t1 - t0)（途中红色部分为梯形，可用梯形面积计算公式计算）然后得到，根据（1）式得到y(t1) = y(t0)+ s，这就是传说中的欧拉法。

那么在计算机中用程序求取已知初值的微分方程解就成了求积分段（积分段与坐标轴）的近似值了。欧拉采用的是用矩形计算法计算红色部分的面积，矩形的高为f0。

 

那么我们就可以编写程序来求解微分方程了。

 

3. 编写程序求解微分方程（C语言）

3.1源代码

#include <stdio.h>//Function decalrefloat function(float t);int Eular(float init, float a, float b, float h, float *result);//Mainint main(){        float a, b, h;        float yt0, result;        //Initiazation        a       = 0;        b       = 3;        h       = 0.01;        yt0     = 0;        Eular(yt0, a, b, h, &result);        return 0;}//Function expressionfloat function(float t){        return (2 * t);}//Eular method//@para a: low limit, b: up limit, h: step//@result compute resutlint Eular(float init, float a, float b, float h, float *result){        int k;        float t;        if(a > b){                return -1;        }        *result = init;        k       = 1;        //Eular method        for(t = a; t <= b; t += h){                *result +=  function(t - h) * h;                printf("%f----%d vaule is %f\n", t, k, *result);                k++;        }        return 0;}

(1) 程序中function函数是微分方程方程式。

(2) Eular函数是用欧拉方法求解微分方程的数值解。当然题目中y(t) = t的平方曲线特殊，可以直接用梯形得到准确解，为了体现欧拉法，给函数function传值的时候需要将横坐标传为t - h，这样f 的值才会对应f0.然后采用y(t1) = y(t0) + s的方式连续的求取下一刻原函数的值。程序中用迭代语句*result  += function(t - h) * h;语句来实现。并打印出响应的数值解。

(3) 关于函数高的选择。题目中式选取的欧拉值即选取t0前一刻的值作为高，也可以选t1后一刻的值作为高，也可以选t0 和 t1的终点作为高。当然还可以用梯形的面积代替积分面积（梯形法）。

 

3.2 程序运行结果

 

h = 0.01时求得数值解

程序中的h值表示步长，当然h的值越小得到的值将越精确。可计算出原函数当t = 3时，真实值为9.而程序中得到的值为8.969793，这是步长h = 0.01时的结果。其实已经很逼近了，如果对程序的精确度还不满意还可以缩小步长。如当h = 0.001时程序的运行效果为：

 

h = 0.001时求的数值解

当，h = 0.0001时，程序求得的数值解更为精确，误差已经到了10的-3次了。当然步长越小，虽然数值解越精确，但是程序耗时越多的。

 

读《Matlab控制系统仿真》课本。

Note Over.