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

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

　　龙格-库塔(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=%f\r\n",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

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

查看本文来源