四阶Runge-Kutta法解常微分方程

/**

***四阶Runge-Kutta法***

经典格式：

y(n+1) = y(n) + h/6 ( K1 + 2*K2 + 2*K3 + K4 )

K1 = f( x(n) , y(n) )

K2 = f( x(n+1/2) , y(n) + h/2*K1 )

K3 = f( x(n+1/2) , y(n) + h/2*K2 )

K4 = f( x(n+1) , y(n) + h*K3 )

Runge-Kutta法是基于泰勒展开方法，因而需要所求解具有较好的光滑性。

属性：差分方法

《数值分析简明教程》-2 Editon -高等教育出版社- page 105 算法流程图

代码维护：2005.6.14 DragonLord

**/

#include<iostream.h>

#include<stdio.h>

#include<math.h>

/*

举例方程：

y'= y - 2*x / y ( 0<x<1 )

y(0) = 1

*/

double f(double x,double y)

{

double re;

if(x==0)re=1;

else re=y-2*x/y;

return re;

}

int main()

{

double x0,x1,y0,y1,h,k1,k2,k3,k4,y;

int N;

while(cin>>x0>>y0>>h>>N)

{

int n=0;

for(;n<N;n++)

{

x1=x0+h;

y=sqrt(1+2*x1);

k1=f(x0,y0);

k2=f(x0+h/2,y0+h*k1/2);

k3=f(x0+h/2,y0+h*k2/2);

k4=f(x1,y0+h*k3);

y1=y0+h*(k1+2*k2+2*k3+k4)/6;

printf("%.1f %.4f %.4f\n",x1,y1,y);

x0=x1;

y0=y1;

}

}

return 0;

}