微分方程式（ルンゲ・クッタ）

/**********************************/
/* 微分方程式（ルンゲ・クッタ法） */
/*      例：d2x/dt+3dx/dt+2x=1    */
/*      coded by Y.Suganuma       */
/**********************************/
#include <stdio.h>

void snx(double, double *, double *);
double rkg(double, double, double *, double *, double **, int,
           void (*)(double, double *, double *));
int main()
{
	double time, h, *x, *dx, **g;
	int i1, n = 2;
/*
          初期設定
*/
	x  = new double [n];
	dx = new double [n];
	g  = new double * [4];
	for (i1 = 0; i1 < 4; i1++)
		g[i1] = new double [n];
	time = 0.0;
	h    = 0.01;
	x[0] = 0.0;
	x[1] = 0.0;
/*
          計算と出力
*/
	for (i1 = 0; i1 < 101; i1++) {
		printf("time = %f, x = %f\n", time, x[0]);
		time = rkg(time, h, x, dx, g, n, snx);
	}

	return 0;
}

/****************/
/* 微係数の計算 */
/****************/
void snx(double time, double *x, double *dx)
{
	dx[0] = x[1];
	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
}

/*******************************************/
/* ルンゲ・クッタ法  dx/dt=f(t,x)          */
/*      time : 現在の時間                  */
/*      h : 時間刻み幅                     */
/*      x : 現在の状態                     */
/*      dx : 微係数（f(t,x):snxで計算する）*/
/*      g : 作業域（g[4][n]）              */
/*      n : 微分方程式の次数               */
/*      snx : 微係数を計算する関数の名前   */
/*      return : time+h                    */
/*******************************************/
double rkg(double time, double h, double *x, double *dx, double **g,
           int n, void (*sub)(double, double *, double *))
{
	int i1;
	double h2;

	h2 = 0.5 * h;

	(*sub)(time, x, dx);
	for (i1 = 0; i1 < n; i1++)
		g[0][i1] = h * dx[i1];

	time += h2;
	for (i1 = 0; i1 < n; i1++)
		g[1][i1] = x[i1] + 0.5 * g[0][i1];
	(*sub)(time, g[1], dx);
	for (i1 = 0; i1 < n; i1++)
		g[1][i1] = h * dx[i1];

	for (i1 = 0; i1 < n; i1++)
		g[2][i1] = x[i1] + 0.5 * g[1][i1];
	(*sub)(time, g[2], dx);
	for (i1 = 0; i1 < n; i1++)
		g[2][i1] = h * dx[i1];

	time += h2;
	for (i1 = 0; i1 < n; i1++)
		g[3][i1] = x[i1] + g[2][i1];
	(*sub)(time, g[3], dx);
	for (i1 = 0; i1 < n; i1++)
		g[3][i1] = h * dx[i1];

	for (i1 = 0; i1 < n; i1++)
		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;

	return time;
}