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

核心提示：一、 引言在工程计算中我们经常要去解一些常微分方程，虽然在高等数学和其他一些涉及微分方程的专业书籍中介绍了不少类型的常微分方程，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=%frn",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=%frn",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

三、 使用经典龙格-库塔算法进行高精度求解

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。同前几种算法一样，该算法也是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有：

　　yi＋１＝yi＋h*K1

　　K1=f(xi,yi)

　　当用点xi处的斜率近似值K1与右端点xi＋１处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式：

　　yi＋１＝yi＋h*( K1+ K2)/2

　　K1=f(xi,yi)

　　K2=f(xi+h,yi+h*K1)

　　依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：

　　yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6

　　K1=f(xi,yi)

　　K2=f(xi+h/2,yi+h*K1/2)

　　K3=f(xi+h/2,yi+h*K2/2)

　　K4=f(xi+h,yi+h*K3)

　　下面的具体程序实现同改进的欧拉算法类似，只需作些必要的改动，下面将该算法的关键部分代码清单列出：

……
for(float x=0;x<0.6;x+=0.1)
{
r=x+expf(-x);
K1=x-y[i]+1; file://求K1
K2=(x+(float)(0.1/2))-(y[i]+K1*(float)(0.1/2))+1; file://求K2
K3=(x+(float)(0.1/2))-(y[i]+K2*(float)(0.1/2))+1; file://求K3
K4=(x+0.1)-(y[i]+K3*0.1)+1; file://求K4
y[i+1]=y[i]+(float)(0.1*(K1+2*K2+2*K3+K4)/6); file://求yi+1
r=fabs(r-y[i]); file://计算误差
str.Format("y[%d]=%f r=%frn",i,y[i],r);
i++;
msg+=str;
}
AfxMessageBox(msg); file://报告计算结果及误差情况

　　从下表记录的程序运行结果来看，在步长为0.1的情况下所计算出来的常微分方程的数值解是非常精确的，用浮点数进行运算时由近似所引入的误差几乎不会对计算结果产生影响，这样高的精度是非常令人满意的。无论是从计算速度上还是从计算精度要求上，都能非常好的满足工程计算的需要。

xI（各分点）	yI （数值解）	y(xi) （真实值）	| y(xi)- yI | (误差)
0.0	1.000000	1.000000	0.000000
0.1	1.004838	1.004837	0.000001
0.2	1.018731	1.018731	0.000000
0.3	1.040818	1.040818	0.000000
0.4	1.070320	1.070320	0.000000
0.5	1.106531	1.106531	0.000000

　　小结：本文针对工程计算中遇到的常微分方程初值问题的求解，根据实际情况引如了计算机作为辅助计算工具，并根据高等数学的有关知识将欧拉公式、改进的欧拉公式、经典龙格－库塔方法等融入到计算机程序算法之中，充分利用了计算机的速度优势，大大减轻了工程技术人员的劳动强度，同时也使计算结果更加可靠、准确。但在实际应用时，针对不同的工程函数选择合适的求解方法需要有较高的要求，既要考虑到方法的简易，又要减少计算量，同时又不能让截断误差超出指定范围。一般来说经典龙格－库塔算法精确度高又利于计算机编程实现，稳定性也很好，可以考虑作为首选实现算法。