改进Eular方法解常微分方程

/**

***改进Eular方法***

预报：*Y(n+1)=Y(n) + h * f( x(n) , y(n) );

改进：Y(n+1)=Y(n) + h/2 *［ f( x(n) , y(n) ) + f( x(n+1) , *Y(n+1) ) ］

h=x(n+1)-x(n)

平均化形式：

Yp= Yn + hf( X(n),Y(n) )

Yc= Yn + hf( X(n+1),Y(p) )

Y(n+1)= 1/2 * ( Yc+Yp )

属性：差分方法

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

代码维护：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,y0,x1,y1,y;

double yp,yc,h;

int N;

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

{

int n=0;

for(;n<N;n++)

{

x1=x0+h;

y=sqrt(1+2*x1); //精确值

yp=y0+h*f(x0,y0);

yc=y0+h*f(x1,yp);

y1=(yp+yc)/2.0;

printf("%.1f %.4f %.4f\n",x1,y1,y); //输出与精确值比较

x0=x1;

y0=y1;

}

}

return 0;

}