ルンゲクッタ

微分方程式の授業で Runge-Kutta の２次元ベクトルを扱うやつをやっているんだが、なかなかできなくて困る。解が爆発してしまうのだ。ベクトルを扱うのだから、C言語でふつうの構造体を使ってもできるのだが、やはりオブジェクト指向の練習と言うことで Java で組んでみた。

この課題はレポート形式で２つの部分に分かれていて、ひとつめは２種類間の動物で強制し同じ資源を得ようとして競争する場合である。これは単なる２元連立のロジスティック方程式で解く。もうひとつは一方が寄生動物（肉食動物）で、他方が宿主動物（草食動物）の場合である。半分は終わったのだが、後半部分はけっこうむずかしい。後半はロトカ・ボルテラ方程式と呼ばれているらしい。

微分方程式

ものごとを変化する様子を数学的に表現し記述する方法のひとつ。微分の階数により１階、２階などと分類される。

ルンゲクッタ法

微分方程式は一般には解けないので、数値計算的に具体的な値をごりごり計算するよいアルゴリズムである。

ロトカ・ボルテラ方程式

寄生する側とされる側の個体の増減を微分方程式にするときによく用いられる方程式らしい。

参考までにソースコードを乗せておく。ライセンスは GPL です。使う人がいれば、の話。

public class Calc {

	public static void main(String[] args) {

		double x0 = 0.0; // Start
		double xn = 100.0; // End
		int step = 99; // Times

		Vect y = new Vect();
		
		// Initial Value
		y.setVect(24, 84);

		Vect k1 = new Vect();
		Vect k2 = new Vect();
		Vect k3 = new Vect();
		Vect k4 = new Vect();

		double h = (xn - x0) / step; // Steps

		for (double x = x0; x < xn; x += h) {
			k1 = Vect.firstCalc(x, y, h);
			k2 = Vect.secondCalc(x, y, h, k1);
			k3 = Vect.thirdCalc(x, y, h, k2);
			k4 = Vect.fourthCalc(x, y, h, k3);

			System.out.println(x + " " + y.getn1() + " " + y.getn2());
			y = Vect.nextY(y, k1, k2, k3, k4);
		}

	}
}

class Vect {

	private double n1;

	private double n2;

	// Constractor
	Vect() {
		this.n1 = 0;
		this.n2 = 0;
	}

	Vect(int p, int q) {
		this.n1 = p;
		this.n2 = q;
	}

	Vect(double p, double q) {
		this.n1 = p;
		this.n2 = q;
	}

	// Set Variable
	void setVect(double p, double q) {
		this.n1 = p;
		this.n2 = q;
	}

	double getn1() {
		return this.n1;
	}

	double getn2() {
		return this.n2;
	}

	// Ordinary Diffetential Equations
	static Vect diffFunc(double x, Vect y) {
		Vect newy = new Vect();

		double a1 = 0.1;
		double a2 = 0.075;
		double b11 = 0.0007;
		double b12 = 0.001;
		double b21 = 0.0007;
		double b22 = 0.0007;

		double n1 = y.n1 * (a1 - b11 * y.n1 - b12 * y.n2);
		double n2 = y.n2 * (a2 - b21 * y.n1 - b22 * y.n2);

		newy.n1 = n1;
		newy.n2 = n2;

		return newy;
	}

	// k1 = h * f(x, y)
	static Vect firstCalc(double x, Vect y, double h) {
		return diffFunc(x, y).times(h);
	}

	//	k2 = h * f(x + h/2, y + k1/2)
	static Vect secondCalc(double x, Vect y, double h, Vect k1) {
		return diffFunc(x + h / 2.0, y.addHalf(k1)).times(h);
	}

	// k3 = h * f(x + h/2, y + k2/2)
	static Vect thirdCalc(double x, Vect y, double h, Vect k2) {
		return diffFunc(x + h / 2.0, y.addHalf(k2)).times(h);
	}

	// k4 = h * f(x + h, y + k3)
	static Vect fourthCalc(double x, Vect y, double h, Vect k3) {
		return diffFunc(x + h, y.add(k3)).times(h);
	}

	// y[i+1] = y[i] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)
	static Vect nextY(Vect y, Vect k1, Vect k2, Vect k3, Vect k4) {
		y.n1 += (k1.n1 + 2.0 * k2.n1 + 2.0 * k3.n1 + k4.n1) / 6.0;
		y.n2 += (k1.n2 + 2.0 * k2.n2 + 2.0 * k3.n2 + k4.n2) / 6.0;
		return y;
	}

	Vect times(double h) {
		this.n1 *= h;
		this.n2 *= h;
		return this;
	}

	Vect add(Vect h) {
		this.n1 += h.n1;
		this.n2 += h.n2;
		return this;
	}

	Vect addHalf(Vect h) {
		this.n1 += h.n1 / 2.0;
		this.n2 += h.n2 / 2.0;
		return this;
	}

}