// work.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "iostream.h"
#include "stdio.h"
#define f(x,y) (y*y)
#define g(x,y) (x/y)
//f(x,y),g(x,y)為所定義的函數(shù)表達(dá)式
#define N 20

int m;
double h;
double x[N],y[N],yb[N];
int main(int argc, char* argv[])
{	
	int i;
	int choice;
	double a,b,y0;
	void euler(double x[],double y[]);
	void rk(double x[],double y[]);
	cout<<"輸入x的范圍[a,b]，上下界a,b分別為："<<endl;
	cin>>a>>b;
	cout<<"y("<<a<<")=";
	cin>>y0;
	cout<<"輸入將區(qū)間劃分的快數(shù)(0<m<20)m=";
	cin>>m;
	while((m>N)||(m<0))
	{
		cout<<"所定義最大區(qū)間為20，請(qǐng)輸入20以內(nèi)整數(shù)m=";
		cin>>m;
	}
	h=(b-a)/m;
	x[0]=a;
	y[0]=y0;
	
	cout<<"請(qǐng)選擇( 1 / 2 ) :1.改進(jìn)的euler法,2.Runge-kutta法"<<endl;
	cin>>choice;
	if (choice==1)
	{
    	euler(x,y);
	   	for(i=0;i<=m;i++)
		{
			cout<<"x"<<i<<"="<<x[i]<<" , ";
		    cout<<"yb"<<i<<"="<<yb[i]<<" , ";
		    cout<<"y"<<i<<"="<<y[i]<<endl;
		}
	}
	else 
	{
		rk(x,y);
		for(i=0;i<=m;i++)
		{
			cout<<"x"<<i<<"="<<x[i]<<" , ";
		    cout<<"y"<<i<<"="<<y[i]<<endl;
		}
	}
	return 0;
}

	
//改進(jìn)的euler法
void euler(double x[],double y[])
{  	   
	 int i;
	 for(i=1;i<=m;i++)
	 {
		x[i]=x[i-1]+h;
		yb[i]=y[i-1]+h*f(x[i-1],y[i-1]);
		y[i]=y[i-1]+h/2*(f(x[i-1],y[i-1])+f(x[i],yb[i]));
	 }
 }
//四階Runge-Kutta法
void rk(double x[],double y[])
{
	int i;
	double k1,k2,k3,k4;
	for(i=1;i<=m;i++)
	{
		k1=h*g(x[i-1],y[i-1]);
		k2=h*g(x[i-1]+h/2,y[i-1]+k1/2);
		k3=h*g(x[i-1]+h/2,y[i-1]+k2/2);
		k4=h*g(x[i-1]+h,y[i-1]+k3);
		x[i]=x[i-1]+h;
		y[i]=y[i-1]+(k1+2*k2+2*k3+k4)/6;	
	}
}