rkf.cpp

ronge-kutta 数值分析问题

/*input a,b,t,x,delta,M
while k<M do
    compute F1,F2,F3,F4,F5,F6
    compute x,y
    compute err
     t=t+h
    output k,t,x,y,err
    h=0.9h[delta/|e|]1/(1+p)
end do
*/

#include <math.h>
#include <stdio.h>
#include <iostream.h>

double sol(double t){
	return(t/(0.5+log(t)));
}


double f( double t, double x){
	double m;
	m=(t*x-pow(x,2))/pow(t,2);
	return(m);
}

double f1( double t, double x ,double h){
	double m;
	m=h*(f(t,x));
	return(m);
}

double f2(double t, double x ,double h){
	double m;
	m=h* f( t+0.25*h, x+0.25*f1(t,x,h) ) ;
	return(m);
}

double f3(double t, double x  ,double h){
	double m;
	m=h* f( t+0.375*h, x+0.09375*f1(t,x,h)+0.28125*f2(t,x,h) ) ;
	return(m);
}

double f4(double t, double x  ,double h){
	double m;
	m=h* f( t+0.923077*h, x+0.879381*f1(t, x, h)-3.2772*f2(t, x, h) +3.32089*f3(t, x, h) );
	return(m);
}

double f5(double t, double x  ,double h){
	double m;
	m=h* f( t+h, x+2.03241*f1(t, x, h)-8.0*f2(t, x, h) +7.17349*f3(t, x, h) -0.205897*f4(t, x, h)) ;
	return(m);
}

double f6(double t, double x  ,double h){
	double m;
	m=h* f( t+0.5*h, x-0.296296*f1(t, x, h)+2.0 *f2(t, x, h) -1.38168*f3(t, x, h)+0.452973*f4(t, x, h)) -0.275*f5(t, x, h);
	return(m);
}


void main(){
    
	double a,b,x,t,h,delta,y,err;
	int M,k;

L:	cout<<"a=";
	cin>>a;
    cout<<"b=";
	cin>>b;
	cout<<"x=";
	cin>>x;

	if((x>b)||(x<a)) {cout<<"a<x<b!"<<endl; goto L;}

	cout<<"t=";
	cin>>t;
	cout<<"M=";
	cin>>M;

	cout<<"h=";
	cin>>h;
	cout<<"delta=";
	cin>>delta;

    printf("k=  0, t=%f, x=%f\n",t,x);
    
	y=x;
    
	
	for(k=1;k<=M;k++){
		err=0.00277778*f1(t,x,h)-0.0299415*f3(t,x,h)-0.0291999*f4(t,x,h)+0.02*f5(t,x,h)+0.0363636*f6(t,x,h);

		x=x+0.115741*f1(t,x,h)+0.548928*f3(t,x,h)+0.535331*f4(t,x,h)-0.2*f5(t,x,h);
	
        y=y+0.118519*f1(t,y,h)+0.518986*f3(t,y,h)+0.506131*f4(t,y,h)-0.18*f5(t,y,h)+0.0363636*f6(t,y,h);
		
		t=t+h;

		//output
		if((fabs(err)<delta)) break;

		printf("k=%3d, t=%f, x=%f, y=%f, sol=%f, err=%e\n", k,t,x,y,sol(t),err);

		h=0.9*h*pow((delta/fabs(err)),0.2);        

	}//for
    
 
}//main