#include<stdio.h>
#include<math.h>
#include<iostream.h>
#include<iomanip.h>
double f(double x,double y,int flag)    
{	double t;
	if (flag==1) t=-x*y*y;
	if (flag==2) t=-x*x+y*y;
	return(t);
}
void DevEuler(double a,double b,double x0,double y0,double h,int flag)   //改進的歐拉方法
{
    cout<<"用改進的歐拉方法解得:"<<endl;
	int k,N;
	double x,y,Yp,Yc;
	N=(int )((b-a)/h+0.5);
	cout<<setw(5)<<"x"<<setw(15)<<"y"<<endl;
	for(k=0;k<N;k++)
	{
		x=x0+h;
		Yp=y0+h*f(x0,y0,flag);
		Yc=y0+h*f(x,Yp,flag);
		y=(Yp+Yc)/2;
	    cout<<setw(5)<< x <<setw(15)<< y <<endl;
		x0=x;
		y0=y;
	}
}

void Runge_Kutta(double a,double b,double x0,double y0,double h,int flag)     //四階龍格-庫塔法
{
	cout<<endl<<"四階龍格-庫塔法解得:"<<endl;
	cout<<setw(5)<<"x"<<setw(15)<<"y"<<endl;
	int k,N;
	double x,y,K1,K2,K3,K4;
	N=(int)((b-a)/h+0.5);
	for(k=0;k<N;k++)
	{
		x=x0+h;
		K1=f(x0,y0,flag);
		K2=f(x0+h/2,y0+K1*h/2,flag);
		K3=f(x0+h/2,y0+K2*h/2,flag);
		K4=f(x,y0+K3*h,flag);
		y=y0+(K1+2*K2+2*K3+K4)*h/6;
	    cout<<setw(5)<< x <<setw(15)<< y <<endl;
		x0=x;
		y0=y;
	}
}


void main()
{
	double h,a,b,x0,y0;
	cout<<"第六單元 常微分方程數值解法"<<endl;
	cout<<"用改進的歐拉方法和四階龍格-庫塔法解:"<<endl;
	cout<<"       (1) y`=-x*y*y, y(0)=2;x屬于(0,5)  "<<endl;
	cout<<"       (2) y`=-x*x+y*y, y(0)=0"<<endl;
	cout<<"第一題求解："<<endl;
	a=0;b=5;
	x0=0;y0=2;
	h=0.5;
    DevEuler(a,b,x0,y0,h,1);
	Runge_Kutta(a,b,x0,y0,h*2,1);
	cout<<endl<<"第二題求解："<<endl;
	a=0;b=1;
	x0=0;y0=0;
	h=0.1;
    DevEuler(a,b,x0,y0,h,2);
	Runge_Kutta(a,b,x0,y0,h*2,2);
}