runge_kutta.cpp

runge-kutta算法

#include<iostream>
using namespace std;
#include<math.h>
double f(double t,double y)
{
	return -0.6*0.1*0.1*sqrt(64.2)*pow(y,-1.5);
}
double rungekutta(double a,double b,int N,double alpha)
{
	double h=(b-a)/N;
	double t=a;
	double w=alpha;
	//cout<<t<<" "<<w<<endl;
	for(int i=1;i<=N;i++)
	{
		double k1=h*f(t,w);
		double k2=h*f(t+h/2,w+k1/2);
		double k3=h*f(t+h/2,w+k2/2);
		double k4=h*f(t+h,w+k3);
		w=w+(k1+2*k2+2*k3+k4)/6;
		t=a+i*h;
		cout<<t<<" "<<w<<endl;
	}
	cout<<t<<" "<<w<<endl;
	return w;
}
int main()
{
	rungekutta(0,600,30,8);
	return 0;
}