龙格_库塔方法

近日翻出大学的数值分析教材，对一些算法重新进行研读，在此记录一些自己的心得和体会吧。这里是第一篇，由于是随意阅读，这个是关于求解一阶常微分方程常用的高精度算法：龙格——库塔方法（Runge Kuttan)。
具体的数值分析的背景知识，完全可以从数值分析的教材获得，这里就不再叙述了。
龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分。
一阶常微分方程可以写作：y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）
Yn+1=Yn+h*f(Xn,Yn)
另外根据微分中值定理，存在0Yn+1=Yn+h*f(Xn+th,Y(Xn+th))
这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。
利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：
Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)；
K1＝f(Xn,Yn);
K2=f(Xn+h/2,Yn+(h/2)*K1);
K3=f(Xn+h/2,Yn+(h/2)*K2);
K4=f(Xn+h,Yn+h*K3);
相应的，我使用C/C++语言具体实现了这个算法，具体算法如下：
//[a,b]为x的取值范围，指针Y存放[a,b]范围内Y的结果数组，step为步进，即公式中的h，fptr为函数。
void Runge_Kutta(double a,double b,double* Y,double step,
double (*fptr)(double,double))
{
double K[4]; //存储4个平均斜率K
double X=a;
double halfstep=step*0.5; //避免多次作除2运算
double coef = 1/6*step; //h*(1/6)为常量
int i=0;
while(X{
K[0]=0; //方便对应，K[0]无实际意义
K[1]=fptr(X,*(Y+i));
K[2]=fptr(X+halfstep,*(Y+i)+halfstep*K[1]);
K[3]=fptr(X+halfstep,*(Y+i)+halfstep*K[2]);
K[4]=fptr(X+step,*(Y+i)+step*K[3]); //这里就具体实现了求K
X+=step; //步进
i++;
//这里就求出了Yn+1的值并将结果保存在数组中
*(Y+i)=*(Y+i-1)+(K[1]+(K[2]+K[3])*2+K[4])*coef;
}
return;
}


这个只是简单的将算法进行了实现，其中还有许多可以优化的地方，比如变步长等可以更好的获得高精度。测试无问题，符合算法提供的截断误差。