runge_kutta.cpp

来自「Runge_Kutta算法 在计算方法中对微分方程进行求解」· C++ 代码 · 共 57 行

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

const double pi = 3.141592653585;
double f1(double t, double u1, double u2)
{
	return u2;
}

double f2(double t, double u1, double u2)
{
	return ((-16.085) * sin(u1));
}

void main()
{
	double t, w1, w2;
	double a, b;
	int N;
	double h;
	double k[5][3];
	a = 0;
	b = 10;
	
	N = 100;

	h = (b-a) / N;

	t = 0;
	w1 = pi / 6;
	w2 = 0;

	cout << t << " " << w1 << " " << w2 << endl;

	for(int i = 1; i <= N; i ++){
		k[1][1] = h * f1(t, w1, w2);
		k[1][2] = h * f2(t, w1, w2);
		k[2][1] = h * f1(t+h/2, w1+k[1][1]/2, w2+k[1][2]/2);
		k[2][2] = h * f2(t+h/2, w1+k[1][1]/2, w2+k[1][2]/2);
		k[3][1] = h * f1(t+h/2, w1+k[2][1]/2, w2+k[2][2]/2);
		k[3][2] = h * f2(t+h/2, w1+k[2][1]/2, w2+k[2][2]/2);
		k[4][1] = h * f1(t+h, w1+k[3][1], w2+k[3][2]);
		k[4][2] = h * f2(t+h, w1+k[3][1], w2+k[3][2]);

		w1 = w1 + (k[1][1] + 2*k[2][1] + 2*k[3][1] + k[4][1])/6;
		w2 = w2 + (k[1][2] + 2*k[2][2] + 2*k[3][2] + k[4][2])/6;

		t = a + i*h;

		cout << t << " " << w1 << " " << w2 << endl;
	}
}