euler.cpp

数值分析的实验报告

#inculde<iostream>
#include<cmath>
using namespace std;
#define N 9											//step length
#define E 0.00001									//presicion 10-5
double g(double x)							//y=g(x)
{return pow(((x*x+2)/3),3.0/2);}
double f(double,double);
double g(double);
void Euler(double,double,double *,double *);
void Implicit_Euler(double,double,double *,double *);
void R_K(double,double,double *,double *);
double f(double x,double y)
{return fabs(sin(2*x));}
double h=0.1;
void Euler(double x0,double y0,double *X,double *Y)
{
	X[0]=x0+h;
	Y[0]=y0+h*f(x0,y0);
	for(int i=1;i<N;i++)
	{
		X[i]=X[0]+i*h;
		Y[i]=Y[i-1]+h*f(X[i-1],Y[i-1]);
	}
}
void Implicit_Euler(double x0,double y0,double *X,double *Y)
{
	double e=0;
	X[0]=x0+h;
	Y[0]=y0+h*f(x0,y0);
	for(int j=0;;j++)
	{
		e=Y[0];
		Y[0]=y0+h*f(X[0],Y[0]);
		if(fabs(e-Y[0])<E) break;
	}
	for(int i=1;i<N;i++)
	{
		X[i]=X[0]+i*h;
		Y[i]=Y[i-1]+h*f(X[i-1],Y[i-1]);
		for(int j=0;;j++)
		{
			e=Y[i];
			Y[i]=Y[i-1]+h*f(X[i],Y[i]);
			if(fabs(e-Y[i])<E) break;
		}
	}
}
int main()
{
	double x0=0,y0=1;
//	double x0=1,y0=1;
//	double x0=0,y0=100;
	double X[N],Y[N];
	int i;
//	Euler(x0,y0,X,Y);
//	cout<<x0<<"  "<<y0<<endl;
//	for(i=0;i<N;i++)
//		cout<<X[i]<<"  "<<Y[i]<<endl;
//	cout<<endl;

//	for(i=0;i<N;i++)
//		cout<<X[i]<<"  "<<g(X[i])<<endl;

	Implicit_Euler(x0,y0,X,Y);
	for(i=0;i<N;i++)
		cout<<X[i]<<"  "<<Y[i]<<endl;
	cout<<endl;
	return 0;
}