rungekutta.cpp

用C++源程序实现经典龙格库塔方法 解四阶微分方程

/*
 * Runge-Kutta Program
 */



#include<iostream.h>


float f(float x, float y,float z)     
    
   
{
	//输入方程右端函数
    return 2*y*y*y;
   
}
    
float RungeKutta_z(float x,float y,float z,float h)   
{
    float k1,k2,k3,k4;
    k1=f(x,y,z);
    k2=f(x+h/2,y+h*z/2,z+h*k1/2);
    k3=f(x+h/2,y+h*(z+h*k1/2)/2,z+h*k2/2);  
    k4=f(x+h,y+h*(z+h*k2/2),z+h*k3);
    return z+h*(k1+2*k2+2*k3+k4)/6;
}

float RungeKutta_y(float x,float y,float z,float h)
{
	float k1,k2,k3;
    k1=f(x,y,z);
    k2=f(x+h/2,y+h*z/2,z+h*k1/2);
    k3=f(x+h/2,y+h*(z+h*k1/2)/2,z+h*k2/2);
	return y+h*z+h*h*(k1+k2+k3)/6;
}

/************
Main function
Entry Point here
*************/
void main()
{ 
	float x,y,z,h,b;
	float xi,yi,zi;
	int i=0;
	cout<<"用龙格－库塔公式求解微分方程的数值解\n";
	cout<<"输入初始条件x0: ";
	cin>>x;
	cout<<"输入初始条件y0: ";
	cin>>y;
	cout<<"输入初始条件z0: ";
	cin>>z;
	cout<<"输入初始条件,步长h: ";
	cin>>h;
	cout<<"输入初始条件,区间右端点b: ";
	cin>>b;

	xi = x;
	yi = y;
	zi = z;
	cout<<"结果为：\n";
	do
	{ 
		i++;
		if(i==1)
			y = RungeKutta_y(x,y,z,h);
		else
		{		
			z = RungeKutta_z(xi,yi,zi,h);
			yi = y;
			y = RungeKutta_y(x,y,z,h);
			xi = x;		
			zi = z;
		}
		x=x+h;  
		cout<<"x"<<i<<" = "<<x<<"	"<<"y"<<i<<" = "<<y<<"\n";
	}while(x<b);

	//getch();
}