/**
***改进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;
}
