Visual C++常微分方程初值问题求解

　　摘要：本文讲述了以计算机为辅助计算工具，分别使用欧拉算法、改进欧拉算法以及经典龙格-库塔算法对常微分方程的初值问题进行数值求解的实现算法。

　　一、 引言

　　在工程计算中我们经常要去解一些常微分方程，虽然在高等数学和其他一些涉及微分方程的专业书籍中介绍了不少类型的常微分方程，及各自的解法。但工程技术人员在工程和科学研究中所关心的往往只是常微分方程的近似数值解，而非从事数学研究的技术人员所注重的"过程"。采用常规的人工推导、求解无疑是效率非常低下的，而且工程上的常微分方程往往结构非常复杂，要给出一般方程解的表达式也是非常困难的。实际上到目前为止，我们只能对有限的几种特殊类型的方程求精确解，这远不能满足工程需要，对那些不能用初等函数来表达的方程就只能去求其近似的数值解，而且这样还可以借助于运算速度快的计算机来进行辅助求解，大大提高求解的速度和精度，修改也比较灵活。

　　二、 使用欧拉算法及其改进算法进行一般求解

　　所谓的数值求解，就是求问题的解y(x)在一系列点上的值y(xi)的近似值yi。欧拉(Euler)算法是其中最基本、最简单的算法，但其求解精度较低，一般不在工程中单独进行计算。其实现的依据是用向前差商来近似代替导数。对于常微分方程：

　　dy/dx=f(x,y)，x∈[a,b]
　　y(a)=y0

　　可以将区间[a,b]分成n段，那么方程在第xI点有y'(xI)=f(xI,y(xI))，再用向前差商近似代替导数则为：(y(xI+1)-y(xI))/h= f(xI,y(xI))，因此可以根据xI点和yI点的数值计算出yI+1来：

　　yI+1= yI+h*f(xI ,yI)

　　下面就在Visual C++ 6.0编程环境下对一个简单的常微分方程

　　y'=x-y+1,x∈[0,0.5]
　　y(0)=1

　　求近似数值解，由于该简单方程可以用数学方法求得其精确描述式y(x)=x+e-x，所以可以据此检验近似数值解同真实解的误差情况。对于其他一些结构复杂的常微分方程的数值解实现方法也是一样的。下面就通过代码来实现上述算法，并对计算结果作了比较：

float y[6]; file://用于存放计算出的常微分方程数值解 float r; file://同真实解的误差情况 memset(y,0,sizeof(float)*6);//清零 y[0]=1; file://y(0)=1 …… for(float x=0;x<0.6;x+=0.1) file://区间分5段，步长为0.1 { r=x+expf(-x); file://真实解y(x)=x+e-x y[i+1]=y[i]+0.1*(x-y[i]+1); file://数值解（近似） r=fabs(r-y[i]); file://误差 str.Format("y[%d]=%f r=%f\r\n",i,y[i],r); i++; msg+=str; } AfxMessageBox(msg); ……

　　经过程序计算，得出y(xi)在各点的近似数值解及各自同真实解的误差，现列表如下，以兹对照：

xI（各分点）	yI （数值解）	y(xi) （真实值）	| y(xi)- yI | (误差)
0.0	1.000000	1.000000	0.000000
0.1	1.000000	1.004837	0.004837
0.2	1.010000	1.018731	0.008731
0.3	1.029000	1.040818	0.011818
0.4	1.056100	1.070320	0.014220
0.5	1.090490	1.106531	0.016041

　　虽然从实验结果看误差不算太大，但这仅仅是一个用于实验的非常简单的常微分方程，对于实际工程中应用的结构复杂的方程其求解结果的误差要远比此大的多，由于还存在着局部截断误差和整体截断误差，有必要采取措施来抑制、减少误差，尽量使结果精确。在构造欧拉公式时采取的一个重要步骤--用向前差商来代替导数，如将其改为向后差商也是行的通的。此时的欧拉公式就变成了：yI+1= yI+h*f(xI+1,yI+1)，由于该式是一个隐式公式，所以可用迭代法进行计算，直至获取到满足精度要求的yI+1。从数学上可以证明，该式的局部截断误差和前面的欧拉公式的截断误差在主部上之相差正负号，所以只要将显示和隐式的两个欧拉公式相加后再行求解会大大减少误差。可以解得改进后的欧拉公式的表达式为：

　　yI+1= yI+h*(f(xI, yI)+f(xI+1, yI+hf(xI,yI)))/2

　　对此式进行编程，就要比前面的代码要麻烦些，需要分步对其进行计算，以达到最高的运算效率，减少运算量：

…… for(float x=0;x<0.6;x+=0.1) { r=x+expf(-x); T1=y[i]+0.1*(x-y[i]+1); file://分步进行计算 T2=y[i]+0.1*((x+0.1)-T1+1); y[i+1]=(T1+T2)/2; r=fabs(r-y[i]); str.Format("y[%d]=%f r=%f\r\n",i,y[i],r); i++; msg+=str; } AfxMessageBox(msg);

　　从下表得出的实验数据可以看出，这种经过改进的欧拉算法所存在的误差已大为减少，可以直接单独应用于实际的工程计算。误差的减少主要是由于先利用了欧拉公式对yI+1的值进行了预估，然后又利用梯形公式对预估值作了校正，从而在预估--校正的过程中减少了误差。

xI（各分点）	yI （数值解）	y(xi) （真实值）	| y(xi)- yI | (误差)
0.0	1.000000	1.000000	0.000000
0.1	1.005000	1.004837	0.000163
0.2	1.019025	1.018731	0.000294
0.3	1.041218	1.040818	0.000400
0.4	1.070802	1.070320	0.000482
0.5	1.107076	1.106531	0.000545