梯形法、龙格库塔法求解微分方程 数值解

根据积分表达式,微分方程的数值解关键在于微分方程的初值及计算微分方程式在tm(上一时刻)与tm+1（下一时刻）与坐标轴围成面积，若这个面积计算得越准确则得到的数值解也就越精确。微分表达式中与坐标轴围成的面积可表示如下，再实施算法的时候可以结合这个图更加直观点：

 

图1 微分方程中微分表达式中tm 与tm+1与坐标中组成阴影面积图

 

1.梯形法求数值解

       显然阴影部分的面积的近似计算中，用梯形法近似比用矩形近似更为准确，梯形的斜边用tm时刻的f的值代替，即原函数在tm时刻的导数。即用梯形法的数值解的表达式变为y(tm+1) = y(tm) + ( f（ y(tm), tm) + f( y1(tm+1), tm+1) ) / 2 * (tm+1 - tm)

从这个表达式来看，如果要用梯形法求得y(tm+1)的值则需要先知道y1(tm+1)的值，这个值可以用欧拉法求得，再把用欧拉法求得的下一时刻的数值解带入上式校正一个更为精确的值。如此，梯形法的求解表达式就变为：

y1(tm+1) =y(tm)  + f((ym, tm)) * (tm+1 - tm)                                                 (1)     [欧拉法求解数值解]

y(tm+1)= y(tm) + ( f（ y(tm), tm) + f(y1(tm+1), tm+1) ) / 2 * (tm+1 - tm)    （2）  [梯形法校正数值解]

 

采用程序求解数值解时初值和步长依然是必不可少的。除此之外用另外一个一个变量来存下一刻值也是必须的。用下面的代码表示：

y_tm    = y_t0;for(t = a; t <= b; t += h){        //Eular value        y1_tm1  = y_tm + function( t - h ) * h;        //Ti xing fa        y_tm1   = y_tm + ( funtion( y_tm, t - h ) + function(y1_tm1, t) ) * h;        //Last value        y_tm    = y_tm1;}

 

       其中y_tm表示上一时刻的数值解，a，b表示上下限，h表示步长值为(tm+1 - tm),y1_tm1为用欧拉法求解的下一刻的数值解，function表示微分方程表达式，用来求tm和tm+1时刻的值作为梯形上底和下底值，y_tm1为用梯形法获得的数值解，一次计算完毕后，ytm1变为上一时刻的数值解，一次求解就能将[a, b]之间的数值解求完毕。

 

2.龙格库塔法求数值解

2.1二阶龙格库塔法

二阶龙格库塔法就是梯形法。它求解的表达式为：

y(tm+1) = y(tm)  + (k1 + k2) / 2 * h;

k1 = f(y(tm), tm)

k2 = f( y(tm+1), tm+1 )

所以求解二阶段的程序段也是一样的。

 

2.2四阶龙格库塔法

二阶龙格库塔法的截断误差正比于步长h的三次方，在计算精度要求比较高的情形下多使用四阶龙格库塔法。其计算公式为

y(tm+1) =y(tm)  + (k1 + 2k2 +2k3 + k4) / 6 *h      （1）

k1 =f( y(tm), tm)                                                   （2）

k2 =f( y(tm)  + k1 / 2 * h, tm + h / 2)                    （3）

k3 =f( y(tm)  + k2 / 2 * h, tm + h / 2)                    （4）

k4 =f( y(tm) + h * k3, tm + h)                               （5）

 

其中y(tm+1)表示下一时刻的数值解，y(tm)表示上一时刻数值解，f()表示微分方程表达式。

仔细分析一下这几个式子，看看是怎么来逼近 微分表达式积分段与坐标抽围成的面积的：

• k1 * h 是对阴影面积的一个逼近，欧拉似逼近。
• k2 * h 亦是对阴影面积的一个逼近，k1 / 2 * h表示用欧拉法求得的是从当前时刻算起跨步长为h / 2的面积，则y(tm)  + k1 / 2 * h表示在 tm + h/2处的数值近似解，则k2代表微分表达式tm + h / 2处的值，则2k2 * h表示计算了阴影面积的两倍。
• k3 *h 更是对阴影面积的逼近，而且比k2 * h更为准确，因为k2 * h基于欧拉法求得原函数在tm + h / 2处的值从而得到微分表达式在tm + h / 2处的值，k3也是微分表达式在tm + h / 2处的值，只是基于k2，更加精确。2k3 * h表示2倍的阴影面积。
• k4 * h是对阴影面积的一个逼近，采用y(tm) + h * k3得到原函数tm+1时刻的数值解从而求得k4 = f( y(tm) + h * k3, tm + h)  为微分表达式在tm + h时刻的解。则k4 * h是表示用tm+1时刻的微分表达式值求得的阴影部分面积，也是用矩阵法求得的，属欧拉法。

 

那么（k1 + 2k2 +2k3 + k4）就表示6倍与阴影面积，所以(k1 + 2k2+2k3 + k4) / 6 * h就表示求得平均值，数学直觉上讲，这样子的到的平均值往往更趋近于真实值。

 

若要用程序来表示龙格库塔法，首先初值和步长依然是导火线。然后就根据(1) ~ (5 )式子列就行了：

float k1, k2, k3, k4;float y_tm, y_tm1;y_tm    = y_t0;for(t = a; t < b; t += h){        //Compute kn        k1      = function(y_tm, t - h);        k2      = function(y_tm + k1 * h / 2, (t - h) + h / 2);        k3      = function(y_tm + k2 * h / 2, (t - h) + h / 2);        k4      = function(y_tm + h * k3, (t - h) + h);        //Next one value        y_tm1   = y_tm + (k1 + 2*k2 + 2*k3 + k4) * h / 6;        //Updata Last value        y_tm    = y_tm1;}

 

y+t0为初值常值，function为分为表达式，t - h表示用欧拉法计算矩形面积（前一时刻的值为矩形高度）。然后求解的过程跟分析的过程相符合。

 

3总结

      龙格库塔法是求解已知初值微分方程的主要数值算法。它还试用于非线性的微分方程的数值解。如果需要描述求解值的曲线，因为这些值都是很准确的，所有采用曲线拟合或图形界面就能将微分方程解的曲线绘制出来。

 

读《Matlab控制系统仿真》。

 

Note Over.