自由落下をオイラー法とルンゲクッタ法で解く

1. 自由落下の運動方程式
2. 解析的手法
3. 解析的手法（跳ね返り）
4. オイラー法
5. オイラー法（跳ね返り）
6. ルンゲクッタ法
7. ルンゲクッタ法（跳ね返り）
8. C言語による実装（跳ね返り無し）
9. Javaによる実装（跳ね返りあり）

以前こちらで雑に考察した自由落下。

かれこれ3年以上前ということで感慨深いものがありましたが、

中途半端に終わっているので決着をつけようと思います。

8. C言語による実装（跳ね返り無し）

9. Javaによる実装（跳ね返りあり）

1. 自由落下の運動方程式

真空中を自由落下する物体には重力のみがはたらくので、次の運動方程式が成立します。

$\displaystyle \frac{d^2}{dt^2}y = -g \ \ \ \ (1)$

　　 $t$

：　時間

　　 $y$

：　位置

　　 $g$

：　重力加速度

この式をオイラー法、ルンゲクッタ法で位置 $y$ について解き、解析的手法との比較を行います。

2. 解析的手法

式(1)を1回積分、2回積分することで、次式が得られます。

$\displaystyle v = v_0 - gt \ \ \ \ (2)$
$\displaystyle y = y_0 + v_0t - \frac{1}{2}gt^2 \ \ \ \ (3)$

ここで $v$ は速度、 $v_0$ は初速度、 $y_0$ は初期位置を表します。

3. 解析的手法（跳ね返り）

式(2)、(3)によってボールの位置が解析的に明らかになりました。

少し複雑になりますが、プログラム上で数値解法と比較するため、

時間ステップ $H$ で進めた際の、解析的な跳ね返りについて考えます。

一先ず落下のみ考えたいという方は、跳ね返りの部分は読み飛ばして下さい。

解析的な跳ね返りは次のような考えで実装しました。

位置がマイナスになったとき、次の処理を行います。

初期位置 $y_0$ 、初速度 $v_0$ から、地面に接する（位置が $0$ になる）時刻 $t_c$ を求める。
時刻 $t_c$ における速度を求め、跳ね返り係数を乗じて初速度 $v_0$ の値を更新する。
初期位置 $y_0$ の値を $0$ に更新する。
地面接触までの累計経過時刻 $t_n$ (初期値は $0$ )に、 $t_c$ を加算する。
現在時刻 $t$ から $t_n$ を減算した時間（ $t - t_n$ ＜ $H$ ）だけ進んだ位置を求め、 $y$ を更新する。
次回ステップ以降の計算においても、跳ね返りの瞬間時刻を起点とするため、 $t_n$ を減算した時間を用いて位置を算出する。

若干ややこしいですが、この考えによって、数値計算と解析的手法の比較が可能になります。

地面に接する時刻 $t_c$ は、式(3)の左辺に $0$ を代入して $t$ について解くことで、次のように求まります。

$\displaystyle t_c = v_0 + \frac{\sqrt{v_0^2+2gy_0}}{g}$

4. オイラー法

次はオイラー法による手法をまとめていきます。

式(1)の二階常微分方程式は、速度 $v$ を用いて次のように一階常微分方程式に分解できます。

$\displaystyle\frac{d}{dt}y = v \ \ \ \ (4)$

$\displaystyle\frac{d}{dt}v = -g \ \ \ \ (5)$

オイラー法では、微分方程式を差分方程式に近似することで、数値計算可能な形にします。

式(4)、(5)は、次のような差分方程式に近似されます。

$\displaystyle y_{n+1} = y_n + v_n \times H \ \ \ \ (6)$

$\displaystyle v_{n+1} = v_n - g \times H \ \ \ \ (7)$

ここで $H$ は時間ステップを表します。

式(6)、(7)の $n$ の値を増やしていくことで、数値的に自由落下の位置を計算できます。

実際のプログラム上ではメモリ節約のため、直前の $y_n$ 、 $v_n$ だけ記憶させておき、

値の更新をすることによって変数を再利用するのが一般的です（多分）。

また、式(6)、(7)では、自由落下（ $v_0 = 0$ ）の際には1番初めのステップを進めても変化がありません。

それに対して、解析解でははじめから少しずつ落下します。ここで時間 $H$ のズレが発生します。

これを解消するために、まずは $v_{n+1}$ の値を更新し、最新の $v_{n+1}$ の値を、 $y_{n+1}$ を求める際に利用します。

数式は以下のようになります。

$\displaystyle v_{n+1} = v_n - g \times H$

$\displaystyle y_{n+1} = y_n + v_{n+1} \times H \ \ \ \ (8)$

5. オイラー法（跳ね返り）

私は昔、位置がマイナスになったら位置を $0$ とし、

速度に $-e$ （反発係数）を乗じるだけ、という方法をとっていましたが、

さすがにお粗末なので、今回はちゃんとした跳ね返りを考えます。

と言っても実装上は解析的手法の跳ね返りよりも簡単で、

位置がマイナスになったとき、次のような手順を行います。

直前の位置 $y_n$ と速度 $v_n$ を用いて、地面に接する（位置が $0$ になる）時刻 $t_c$ を求める。
時刻 $t_c$ を用いて、地面に衝突した瞬間の速度 $v_c$ を計算する。
位置を $0$ 、跳ね返り直後の速度を $-v_c\times e$ とする。
次のステップまでの残り時間 $H-t_c$ だけ時間を進め、 $y_{n+1}$ 、 $v_{n+1}$ を算出する。

以上となります。

（時間を巻き戻して）地面に接する瞬間の時刻、次のステップまでに跳ね返っている位置を求める点で、

解析的手法と同様ではありますが、時間の起点を考慮しなくて良い分考えやすくなっています。

地面に接する時刻 $t_c$ は、式(8)の左辺に $0$ を代入して $t$ について解くことで、次のように求まります。

$\displaystyle t_c = \frac{y_n}{v_{n+1}}$

6. ルンゲクッタ法

そしてルンゲクッタ法。

ルンゲクッタ法は複数勾配の重み付け平均を利用することで、高精度な近似を行う手法です。

今回は一般的な4次のものを利用します。

詳細は割愛して、今回のケースで利用する式をどんどん書いていきます。

二階微分方程式のルンゲクッタ法ということで混乱しそうになるかもしれませんが、

分解した各一階微分方程式の勾配が、計算過程で相互に影響し合うように書きます。

$\displaystyle y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 +2 k_3 + k_4) \ \ \ \ (9)$

$\displaystyle v_{n+1} = v_n + \frac{1}{6}(l_1 + 2l_2 +2 l_3 + l_4) \ \ \ \ (10)$

と定義したとき、各勾配は次のように計算していきます。

$\displaystyle k_1 = H\times v_n$

$\displaystyle l_1 = H\times (-g)$

$\displaystyle k_2 = H\times \left( v_n +\frac{ l_1}{2} \right)$

$\displaystyle l_2 = H\times (-g)$

$\displaystyle k_3 = H\times \left( v_n + \frac{l_2}{2} \right)$

$\displaystyle l_3 = H\times (-g)$

$\displaystyle k_4 = H\times (v_n + l_3)$

$\displaystyle l_4 = H\times (-g)$

これら式を用いることで、ルンゲクッタ法による自由落下を計算することができます。

7. ルンゲクッタ法（跳ね返り）

ルンゲクッタ法による跳ね返りは、オイラー法の時と同様の流れになります。

しかし、地面に接する時刻 $t_c$ を求めるのは若干ややこしくなります。

勾配の式を整理しながら式(9)に代入していくと、次のようになります。

$\displaystyle y_{n+1} = \frac{v_n H + 4\left( v_n - \frac{gH}{2} \right) + H(v_n - gH)}{6}$

$\displaystyle y_{n+1} = y_n + v_n H - \frac{1}{2}gH^2 \ \ \ \ (11)$

式(11)の左辺を $0$ とすると、その時の $H$ が $t_c$ であるため、この連立方程式を解くと

$\displaystyle t_c = \frac{v_n + \sqrt{v_n^2 - 2gy_n}}{g}$

となります。

・・・ここまで計算して、やってしまった感が出てきました。

4次のルンゲクッタの式を変形した結果、 $t_c$ が解析解のものと一致しています。

自由落下の問題において、4次のルンゲクッタによる数値解は、解析解に等しいということに帰着してしまいそうです。

そもそも式(11)の段階でそのことが半分証明されていました。もう証明してしまいましょう。

勾配の式を整理しながら式(10)に代入していくと、次のようになります。

$\displaystyle v_{n+1} = v_n - gH \ \ \ \ (12)$

式(12)は公差 $-gH$ の等差数列なので、初項を $v_0$ とすると

$v_n = v_0 - gH(n-1)$

となります。

また、式(11)の $n$ の値をずらしながら方程式を立て、加算していくことで次式が得られます。

$\displaystyle y_{n+1} = y_0 + \sum_{i=0}^n v_n H - \frac{1}{2}gH^2n$

$\displaystyle y_{n+1} = y_0 + \sum_{i=0}^n \{v_0-gH(n-1) \} H - \frac{1}{2}gH^2n$

$\displaystyle y_{n+1} = y_0 + (v_0 + gH) Hn - \sum_{i=0}^n gH^2n - \frac{1}{2}gH^2n$

$\displaystyle y_{n+1} = y_0 + v_0 Hn + gH^2n - gH^2\frac{(n+1)n}{2} - \frac{1}{2}gH^2n$

$\displaystyle y_{n+1} = y_0 + v_0 Hn - \frac{1}{2}g(Hn)^2$

ここで、時刻 $t$ は、ステップ $H$ を $n$ 回進めた時間 $Hn$ に等しいため、

$\displaystyle y_{n+1} = y_0 + v_0 t - \frac{1}{2}gt^2$

以上により、自由落下の解析解と、ルンゲクッタによる数値解が等しいことが証明されました。

証明されてしまいましたが、見なかったことにして実装していきます。

8. C言語による実装（跳ね返り無し）

跳ね返りを考えない自由落下を、C言語で実装すると次のようになります。

計算結果を「freefall.csv」に書き出すプログラムとなっています。

#include <stdio.h>
#define G 9.80665 		//Gravitational acceleration
#define H 0.05			//Time step

typedef struct ball{
	double y, v, y0, v0;
}BALL;

/*	Prototype declaration	*/
void calc_real(BALL *ball, double t);
void calc_euler(BALL *ball);
void calc_runge(BALL *ball);

int main(void){
	FILE *fp;
	fp= fopen("freefall.csv", "w");

	double t = 0; 					//Time
	BALL real = {100, 0, 100, 0}; 	//Analytical method
	BALL euler = {100, 0, 100, 0};	//Euler method
	BALL runge = {100, 0, 100, 0};	//Runge-Kutta metod

	/*	Write each initial values	*/
	fprintf(fp, "time,Analysis,Euler,Runge-Kutta\n");
	fprintf(fp, "%.3lf,%.14lf,%.14lf,%.14lf\n", t, real.y, euler.y, runge.y);

	while (t < 5){
		/*	Advance time step	*/
		t += H;

		calc_real(&real, t);   	//Analytical method
		calc_euler(&euler);		//Euler method
		calc_runge(&runge);		//Runge-Kutta method

		/*	Write each results	*/
		fprintf(fp, "%.3lf,%.14lf,%.14lf,%.14lf\n", t, real.y, euler.y, runge.y);
	}

	fclose(fp);
	return 0;
}

/*	Analytical method	*/
void calc_real(BALL *ball, double t){
	ball->y = ball->y0 + ball->v0 * t - 0.5 * G * t * t;
}

/*	Euler method	*/
void calc_euler(BALL *ball){
	ball->v = ball->v0 - G * H;
	ball->y = ball->y0 + ball->v * H;

	/*	Update values	*/
	ball->y0 = ball->y;
	ball->v0 = ball->v;
}

/*	Runge-Kutta method	*/
void calc_runge(BALL *ball){
	double k[4], l[4];
	k[0] = H * ball->v0;
	l[0] = H * (-G);
	k[1] = H * (ball->v0 + l[0] / 2.0);
	l[1] = H * (-G);
	k[2] = H * (ball->v0 + l[1] / 2.0);
	l[2] = H * (-G);
	k[3] = H * (ball->v0 + l[2]);
	l[3] = H * (-G);

	ball->y = ball->y0 + (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6.0;
	ball->v = ball->v0 + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;

	/*	Update values	*/
	ball->y0 = ball->y;
	ball->v0 = ball->v;
}

#include <stdio.h>

#define G 9.80665 //Gravitational acceleration

#define H 0.05 //Time step

typedef struct ball{

double y, v, y0, v0;

}BALL;

/* Prototype declaration */

void calc_real(BALL *ball, double t);

void calc_euler(BALL *ball);

void calc_runge(BALL *ball);

int main(void){

FILE *fp;

fp= fopen("freefall.csv", "w");

double t = 0; //Time

BALL real = {100, 0, 100, 0}; //Analytical method

BALL euler = {100, 0, 100, 0}; //Euler method

BALL runge = {100, 0, 100, 0}; //Runge-Kutta metod

/* Write each initial values */

fprintf(fp, "time,Analysis,Euler,Runge-Kutta\n");

fprintf(fp, "%.3lf,%.14lf,%.14lf,%.14lf\n", t, real.y, euler.y, runge.y);

while (t < 5){

/* Advance time step */

t += H;

calc_real(&real, t); //Analytical method

calc_euler(&euler); //Euler method

calc_runge(&runge); //Runge-Kutta method

/* Write each results */

fprintf(fp, "%.3lf,%.14lf,%.14lf,%.14lf\n", t, real.y, euler.y, runge.y);

}

fclose(fp);

return 0;

}

/* Analytical method */

void calc_real(BALL *ball, double t){

ball->y = ball->y0 + ball->v0 * t - 0.5 * G * t * t;

}

/* Euler method */

void calc_euler(BALL *ball){

ball->v = ball->v0 - G * H;

ball->y = ball->y0 + ball->v * H;

/* Update values */

ball->y0 = ball->y;

ball->v0 = ball->v;

}

/* Runge-Kutta method */

void calc_runge(BALL *ball){

double k[4], l[4];

k[0] = H * ball->v0;

l[0] = H * (-G);

k[1] = H * (ball->v0 + l[0] / 2.0);

l[1] = H * (-G);

k[2] = H * (ball->v0 + l[1] / 2.0);

l[2] = H * (-G);

k[3] = H * (ball->v0 + l[2]);

l[3] = H * (-G);

ball->y = ball->y0 + (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6.0;

ball->v = ball->v0 + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;

/* Update values */

ball->y0 = ball->y;

ball->v0 = ball->v;

}

自由落下だけということもあり、割と短いソースコードになりました。

※コードをコピーする際は、タブ文字がスペースに置換されない単純ビュー（画面ダブルクリック）にすることをお勧めします。

書き出されたcsvをExcelで開き、グラフ化してみます。

解析解とオイラー法では僅かに誤差が観察されました。

ルンゲクッタ法とは線分が完全に一致してしまっています。

次に、解析解との差分を計算し、誤差曲線を描いてみました。

オイラー法では誤差は線形に増加していくのがわかります。

同じスケールでは、ルンゲクッタ法の誤差は全くわかりません。

拡大すると次のようになりました。

10のマイナス13乗のスケールでこつこつと誤差が増えていきます。

これは計算順序の違いによって発生したdouble型の精度の問題と考えられます。

9. Javaによる実装（跳ね返りあり）

次は、Javaで跳ね返りありのアニメーション付きで実装してみました。

最近NetBeansのGUIエディタに頼ってばっかりだったので、今回はEclipseでせっせとSwingコンポーネント書きました。

import java.awt.Color;
import java.awt.Dimension;
import java.awt.FlowLayout;
import java.awt.Font;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.RenderingHints;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.text.DecimalFormat;

import javax.swing.JButton;
import javax.swing.JCheckBox;
import javax.swing.JLabel;
import javax.swing.JPanel;
import javax.swing.JSpinner;
import javax.swing.SpinnerNumberModel;
import javax.swing.border.EtchedBorder;
import javax.swing.border.TitledBorder;

public class FreeFall extends javax.swing.JFrame {

	double t = 0.0, tn = 0.0; 		//Time
	int ground = 70;    	  		//Ground heights
	double G = 9.80665; 	  		//Gravitational acceleration
	double E = 1.00;      			//Restitution coefficient
	double H = 0.05;    	  		//Time step
	Ball real, euler, runge;	  	//Ball instance
	boolean start = false;			//Animation flag
	/*	Draw line flag	*/
	boolean drawLineAna = false, drawLineEuler = false, drawLineRunge = false;
	/*	Parameter components	*/
	public static JSpinner spPosition, spVelocity, spBall, spE, spG, spH;

	/*	Ball information class	 */
	class Ball {
		double r = 0, x = 0, y = 0, vy = 0;
		double y0 = 0, vy0 = 0;

		/*	Constructor	*/
		Ball(double x) {
			this.x = x;
			reset();
		}

		/*	Return x drawing position	*/
		int getX() {
			return (int) Math.round(x - r);
		}

		/*	Return y drawing position	*/
		int getY() {
			return (int) Math.round(y - r) + ground;
		}

		/*	Return ball diameter	*/
		int getD() {
			return (int) Math.round(2 * r);
		}

		/*	Initialization	*/
		void reset() {
			r = (Integer) spBall.getValue();
			y = (Integer) spPosition.getValue() + r;
			y0 = y;
			vy = (Double) spVelocity.getValue();
			vy0 = vy;

		}
	}

	/*	Analytical method	*/
	public void calcReal(Ball ball) {
		if ((ball.y0 - ball.r) == 0 && ball.vy0 == 0)	return;

		ball.y = ball.y0 + ball.vy0 * (t - tn) - 0.5 * G * (t - tn) * (t - tn);
		ball.vy = ball.vy0 - G * (t - tn);

		/*	Rebound process	*/
		if (ball.y - ball.r < 0) {
			double tc = (ball.vy0 + Math.sqrt(ball.vy0 * ball.vy0 + 2 * G * (ball.y0 - ball.r))) / G;
			ball.vy0 = -(ball.vy0 - G * tc) * E;
			tn += tc;
			ball.y0 = ball.r;
			ball.y = ball.y0 + ball.vy0 * (t - tn) - 0.5 * G * (t - tn) * (t - tn);
			ball.vy = ball.vy0 - G * (t - tn);
		}
		
		/*	Stationary judgment	*/
		if (ball.y - ball.r < 0) {
			ball.vy = ball.vy0 = 0;
			ball.y = ball.y0 = ball.r;
		}
	}

	/*	Euler method	*/
	public void calcEuler(Ball ball) {
		if ((ball.y0 - ball.r) == 0 && ball.vy0 == 0)	return;

		ball.vy = ball.vy0 - G * H;
		ball.y = ball.y0 + ball.vy * H;

		/*	Rebound process	*/
		if (ball.y - ball.r < 0) {
			double diff = (ball.r - ball.y0) / ball.vy;
			ball.vy = ball.vy - G * diff;

			ball.y = ball.r;
			ball.vy *= -E;

			ball.vy = ball.vy - G * (H - diff);
			ball.y = ball.y + ball.vy * (H - diff);
		}
		
		/*	Update Values	*/
		ball.y0 = ball.y;
		ball.vy0 = ball.vy;

		/*	Stationary judgment	*/
		if (ball.y - ball.r < 0) {
			ball.vy = ball.vy0 = 0;
			ball.y = ball.y0 = ball.r;
		}
	}

	/*	Runge-Kutta method	*/
	public void calcRunge(Ball ball) {
		if ((ball.y0 - ball.r) == 0 && ball.vy0 == 0)	return;

		double[] k = new double[4];
		double[] l = new double[4];
		
		k[0] = H * ball.vy0;
		l[0] = H * (-G);

		k[1] = H * (ball.vy0 + l[0] / 2.0);
		l[1] = H * (-G);

		k[2] = H * (ball.vy0 + l[1] / 2.0);
		l[2] = H * (-G);

		k[3] = H * (ball.vy0 + l[2]);
		l[3] = H * (-G);

		ball.y = ball.y0 + (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6.0;
		ball.vy = ball.vy0 + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;

		/*	Rebound process	*/
		if (ball.y - ball.r < 0) {
			double diff = (ball.vy0 + Math.sqrt(ball.vy0*ball.vy0-2*G*(ball.r-ball.y0))) / G;
			
			k[0] = diff * ball.vy0;
			l[0] = diff * (-G);

			k[1] = diff * (ball.vy0 + l[0] / 2.0);
			l[1] = diff * (-G);

			k[2] = diff * (ball.vy0 + l[1] / 2.0);
			l[2] = diff * (-G);

			k[3] = diff * (ball.vy0 + l[2]);
			l[3] = diff * (-G);

			ball.vy = ball.vy0 + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;

			
			ball.y = ball.r;
			ball.vy *= -E;

			k[0] = (H - diff) * ball.vy;
			l[0] = (H - diff) * (-G);

			k[1] = (H - diff) * (ball.vy + l[0] / 2.0);
			l[1] = (H - diff) * (-G);

			k[2] = (H - diff) * (ball.vy + l[1] / 2.0);
			l[2] = (H - diff) * (-G);

			k[3] = (H - diff) * (ball.vy + l[2]);
			l[3] = (H - diff) * (-G);

			ball.y = ball.y + (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6.0;
			ball.vy = ball.vy + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;
		}

		/*	Update Values	*/
		ball.y0 = ball.y;
		ball.vy0 = ball.vy;

		/*	Stationary judgment	*/
		if (ball.y - ball.r < 0) {
			ball.vy = ball.vy0 = 0;
			ball.y = ball.y0 = ball.r;
		}
	}

	/*	Drawing class	*/
	class DrawComponent extends JPanel {

		@Override
		public void paintComponent(Graphics g) {
			Graphics2D g2 = (Graphics2D) g;
			/*	Antialiasing	*/
			g2.setRenderingHint(RenderingHints.KEY_ANTIALIASING,
					RenderingHints.VALUE_ANTIALIAS_ON);

			/*	Fill the screen	*/
			g2.setColor(Color.white);
			g2.fillRect(0, 0, getWidth(), getHeight());

			/*	Draw information	*/
			g2.setFont(new Font(null, Font.BOLD, 13));
			g2.setColor(Color.black);
			g2.drawString("Analysis", real.getX(), getHeight() - 50);
			g2.drawString("Euler", euler.getX() + 10, getHeight() - 50);
			g2.drawString("Runge-Kutta", runge.getX() - 10, getHeight() - 50);
			g2.drawString("pos", 15, getHeight() - 32);
			g2.drawString("vel", 15, getHeight() - 12);

			DecimalFormat df0 = new DecimalFormat("00.000");
			DecimalFormat df = new DecimalFormat("000.000000");
			g2.setFont(new Font(null, Font.PLAIN, 13));
			g2.drawString("t = " + df0.format(t), 10, 20);
			g2.drawString("" + df.format(real.y - real.r), real.getX() - 5, getHeight() - 30);
			g2.drawString("" + df.format(euler.y - euler.r), euler.getX() - 5, getHeight() - 30);
			g2.drawString("" + df.format(runge.y - runge.r), runge.getX() - 5, getHeight() - 30);

			g2.drawString("" + df.format(real.vy), real.getX() - 5, getHeight() - 10);
			g2.drawString("" + df.format(euler.vy), euler.getX() - 5, getHeight() - 10);
			g2.drawString("" + df.format(runge.vy), runge.getX() - 5, getHeight() - 10);

			/*	Coordinate transformation	*/
			g2.translate(0, getHeight() - 1);
			g2.scale(1, -1);

			/*	Draw the ground	*/
			g2.drawLine(0, ground, getWidth(), ground);

			/*	Draw the balls	*/
			g2.drawOval(real.getX(), real.getY(), real.getD(), real.getD());
			g2.drawOval(euler.getX(), euler.getY(), euler.getD(), euler.getD());
			g2.drawOval(runge.getX(), runge.getY(), runge.getD(), runge.getD());

			/* Draw the ball line	*/
			g2.setColor(Color.gray);
			if (drawLineAna) {
				g2.drawLine(0, real.getY(), getWidth(), real.getY());
			}
			if (drawLineEuler) {
				g2.drawLine(0, euler.getY(), getWidth(), euler.getY());
			}
			if (drawLineRunge) {
				g2.drawLine(0, runge.getY(), getWidth(), runge.getY());
			}
		}
	}

	/*	Animation class	*/
	class Animation extends Thread {
		long delay = 1;
		Animation(){
			delay = (long)(H / 0.01);
			if (delay == 0) delay = 1;
		}
		
		@Override
		public void run() {
			while (start) {
				/*	Advance time step	*/
				t += H;

				/*	Update the ball position	*/
				calcReal(real);
				calcEuler(euler);
				calcRunge(runge);

				repaint();
				try {
					Thread.sleep(delay);
				} catch (InterruptedException e) {
				}
			}
		}

	}

	/*	Constructor to set the GUI	*/
	FreeFall(String title) {
		/*	Set the JFrame	*/
		setTitle(title);
		setDefaultCloseOperation(EXIT_ON_CLOSE);
		setSize(640, 480);
		setLocationRelativeTo(null);
		setLayout(new FlowLayout());

		/*	Add the draw panel on the JFrame	*/
		DrawComponent panelDraw = new DrawComponent();
		panelDraw.setPreferredSize(new Dimension(500, 480));
		panelDraw.setBorder(new EtchedBorder());
		add(panelDraw);

		/*	Add the control panel on the JFrame	*/
		JPanel panelControl = new JPanel();
		panelControl.setPreferredSize(new Dimension(230, 480));
		add(panelControl);

		{
			
			/*	Add the parameter panel on the control panel	*/
			JPanel panelParameter = new JPanel();
			panelParameter.setPreferredSize(new Dimension(230, 270));
			panelParameter.setBorder(new TitledBorder("Parameters"));
			panelParameter.setLayout(new FlowLayout());
			panelControl.add(panelParameter);
			{

				/*	Add the initial position form on the parameter panel		*/
				JLabel labelPositon = new JLabel("Initial position (y0)");
				labelPositon.setPreferredSize(new Dimension(155, 30));
				panelParameter.add(labelPositon);
				spPosition = new JSpinner();
				spPosition.setModel(new SpinnerNumberModel(300, 0, null, 1));
				spPosition.setPreferredSize(new Dimension(50, 20));
				panelParameter.add(spPosition);

				/*	Add the initial velocity form on the parameter panell	*/
				JLabel labelVelocity = new JLabel("Initial velocity (vy0)");
				labelVelocity.setPreferredSize(new Dimension(155, 30));
				panelParameter.add(labelVelocity);
				spVelocity = new JSpinner();
				spVelocity.setModel(new SpinnerNumberModel(0d, null, null, 1d));
				spVelocity.setPreferredSize(new Dimension(50, 20));
				panelParameter.add(spVelocity);

				/*	Add the ball radious form on the parameter panel	*/
				JLabel labelBall = new JLabel("Ball radious (r)");
				labelBall.setPreferredSize(new Dimension(155, 30));
				panelParameter.add(labelBall);
				spBall = new JSpinner();
				spBall.setModel(new SpinnerNumberModel(30, null, null, 1));
				spBall.setPreferredSize(new Dimension(50, 20));
				panelParameter.add(spBall);

				/*	Add the restitution coefficient form on the parameter panel	*/
				JLabel labelE = new JLabel("Restitution coefficient (E)");
				labelE.setPreferredSize(new Dimension(155, 30));
				panelParameter.add(labelE);
				spE = new JSpinner();
				spE.setModel(new SpinnerNumberModel(1.0d, 0d, 1d, 0.1d));
				spE.setPreferredSize(new Dimension(50, 20));
				panelParameter.add(spE);

				/*	Add the gravitational acceleration form on the parameter panel	*/
				JLabel labelG = new JLabel("Gravitational acceleration (G)");
				labelG.setPreferredSize(new Dimension(155, 30));
				panelParameter.add(labelG);
				spG = new JSpinner();
				spG.setModel(new SpinnerNumberModel(9.80665d, null, null, 0.1d));
				spG.setPreferredSize(new Dimension(50, 20));
				panelParameter.add(spG);

				/*	Add the time step form on the parameter panel	*/
				JLabel labelH = new JLabel("Step size (H)");
				labelH.setPreferredSize(new Dimension(155, 30));
				panelParameter.add(labelH);
				spH = new JSpinner();
				spH.setModel(new SpinnerNumberModel(0.05d, 0d, null, 0.005d));
				spH.setPreferredSize(new Dimension(50, 20));
				panelParameter.add(spH);

				JLabel indent = new JLabel();
				indent.setPreferredSize(new Dimension(85, 20));
				panelParameter.add(indent);

				/*	Add the apply button on the parameter panel	*/
				JButton btnApply = new JButton("Apply");
				btnApply.setPreferredSize(new Dimension(120, 30));
				btnApply.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						E = (Double) spE.getValue();
						G = (Double) spG.getValue();
						H = (Double) spH.getValue();
						/*	Initialize the time	*/
						t = 0;
						tn = 0;
						/*	Initialize the ball position	*/
						real.reset();
						euler.reset();
						runge.reset();
						start = false;
						repaint();
					}
				});
				panelParameter.add(btnApply);
			}

			/*	Add the view setting panel on the control panel	*/
			JPanel panelView = new JPanel();
			panelView.setPreferredSize(new Dimension(230, 100));
			panelView.setBorder(new TitledBorder("View setting"));
			panelView.setLayout(new FlowLayout());
			panelControl.add(panelView);
			{
				JCheckBox checkAna = new JCheckBox("Draw the analysis ball line");
				checkAna.setPreferredSize(new Dimension(200, 20));
				checkAna.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						drawLineAna = !drawLineAna;
						repaint();
					}
				});
				panelView.add(checkAna);

				JCheckBox checkEuler = new JCheckBox("Draw the Euler ball line");
				checkEuler.setPreferredSize(new Dimension(200, 20));
				checkEuler.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						drawLineEuler = !drawLineEuler;
						repaint();
					}
				});
				panelView.add(checkEuler);

				JCheckBox checkRunge = new JCheckBox("Draw the Runge-Kutta ball line");
				checkRunge.setPreferredSize(new Dimension(200, 20));
				checkRunge.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						drawLineRunge = !drawLineRunge;
						repaint();
					}
				});
				panelView.add(checkRunge);
			}

			/*	Add the simulation panel on the control panel	*/
			JPanel panelSimulation = new JPanel();
			panelSimulation.setPreferredSize(new Dimension(230, 95));
			panelSimulation.setBorder(new TitledBorder("Simulation"));
			panelSimulation.setLayout(new FlowLayout());
			panelControl.add(panelSimulation);
			{
				/*	Add the start button on the simulation panel	*/
				JButton btnStart = new JButton("Start");
				btnStart.setPreferredSize(new Dimension(100, 30));
				btnStart.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						if (!start) {
							start = true;
							Animation anime = new Animation();
							anime.start();
						}
					}
				});
				panelSimulation.add(btnStart);

				/*	Add the stop button on the simulation panel	*/
				JButton btnStop = new JButton("Stop");
				btnStop.setPreferredSize(new Dimension(100, 30));
				btnStop.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						start = false;
					}
				});
				panelSimulation.add(btnStop);

				/*	Add the reset button on the simulation panel	*/
				JButton btnReset = new JButton("Reset");
				btnReset.setPreferredSize(new Dimension(100, 30));
				btnReset.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						E = (Double) spE.getValue();
						G = (Double) spG.getValue();
						H = (Double) spH.getValue();
						/*	Initialize the time	*/
						t = 0;
						tn = 0;
						/*	Initialize the ball position	*/
						real.reset();
						euler.reset();
						runge.reset();
						start = false;
						repaint();
					}
				});
				panelSimulation.add(btnReset);

				/*	Add the step by step button on the simulation panel	*/
				JButton btnStep = new JButton("Step by step");
				btnStep.setPreferredSize(new Dimension(100, 30));
				btnStep.addActionListener(new ActionListener() {
					public void actionPerformed(ActionEvent e) {
						/*	Advance time step	*/
						t += H;

						/*	Update the ball position	*/
						calcReal(real);
						calcEuler(euler);
						calcRunge(runge);

						repaint();
					}
				});
				panelSimulation.add(btnStep);
			}
		}
		pack();

		/*	Initialize the balls	*/
		real = new Ball(100);
		euler = new Ball(250);
		runge = new Ball(400);
	}

	public static void main(String[] args) {
		/*	Set the Windows look and feel	*/
		try {
			javax.swing.UIManager
					.setLookAndFeel("com.sun.java.swing.plaf.windows.WindowsLookAndFeel");
		} catch (Exception e) {
		}

		/*	Create and display the form	*/
		FreeFall frame = new FreeFall("FreeFall");
		frame.setVisible(true);
	}
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

import java.awt.Color;

import java.awt.Dimension;

import java.awt.FlowLayout;

import java.awt.Font;

import java.awt.Graphics;

import java.awt.Graphics2D;

import java.awt.RenderingHints;

import java.awt.event.ActionEvent;

import java.awt.event.ActionListener;

import java.text.DecimalFormat;

import javax.swing.JButton;

import javax.swing.JCheckBox;

import javax.swing.JLabel;

import javax.swing.JPanel;

import javax.swing.JSpinner;

import javax.swing.SpinnerNumberModel;

import javax.swing.border.EtchedBorder;

import javax.swing.border.TitledBorder;

public class FreeFall extends javax.swing.JFrame {

double t = 0.0, tn = 0.0; //Time

int ground = 70; //Ground heights

double G = 9.80665; //Gravitational acceleration

double E = 1.00; //Restitution coefficient

double H = 0.05; //Time step

Ball real, euler, runge; //Ball instance

boolean start = false; //Animation flag

/* Draw line flag */

boolean drawLineAna = false, drawLineEuler = false, drawLineRunge = false;

/* Parameter components */

public static JSpinner spPosition, spVelocity, spBall, spE, spG, spH;

/* Ball information class */

class Ball {

double r = 0, x = 0, y = 0, vy = 0;

double y0 = 0, vy0 = 0;

/* Constructor */

Ball(double x) {

this.x = x;

reset();

}

/* Return x drawing position */

int getX() {

return (int) Math.round(x - r);

}

/* Return y drawing position */

int getY() {

return (int) Math.round(y - r) + ground;

}

/* Return ball diameter */

int getD() {

return (int) Math.round(2 * r);

}

/* Initialization */

void reset() {

r = (Integer) spBall.getValue();

y = (Integer) spPosition.getValue() + r;

y0 = y;

vy = (Double) spVelocity.getValue();

vy0 = vy;

}

/* Analytical method */

public void calcReal(Ball ball) {

if ((ball.y0 - ball.r) == 0 && ball.vy0 == 0) return;

ball.y = ball.y0 + ball.vy0 * (t - tn) - 0.5 * G * (t - tn) * (t - tn);

ball.vy = ball.vy0 - G * (t - tn);

/* Rebound process */

if (ball.y - ball.r < 0) {

double tc = (ball.vy0 + Math.sqrt(ball.vy0 * ball.vy0 + 2 * G * (ball.y0 - ball.r))) / G;

ball.vy0 = -(ball.vy0 - G * tc) * E;

tn += tc;

ball.y0 = ball.r;

ball.y = ball.y0 + ball.vy0 * (t - tn) - 0.5 * G * (t - tn) * (t - tn);

ball.vy = ball.vy0 - G * (t - tn);

}

/* Stationary judgment */

if (ball.y - ball.r < 0) {

ball.vy = ball.vy0 = 0;

ball.y = ball.y0 = ball.r;

}

/* Euler method */

public void calcEuler(Ball ball) {

if ((ball.y0 - ball.r) == 0 && ball.vy0 == 0) return;

ball.vy = ball.vy0 - G * H;

ball.y = ball.y0 + ball.vy * H;

/* Rebound process */

if (ball.y - ball.r < 0) {

double diff = (ball.r - ball.y0) / ball.vy;

ball.vy = ball.vy - G * diff;

ball.y = ball.r;

ball.vy *= -E;

ball.vy = ball.vy - G * (H - diff);

ball.y = ball.y + ball.vy * (H - diff);

}

/* Update Values */

ball.y0 = ball.y;

ball.vy0 = ball.vy;

/* Stationary judgment */

if (ball.y - ball.r < 0) {

ball.vy = ball.vy0 = 0;

ball.y = ball.y0 = ball.r;

}

/* Runge-Kutta method */

public void calcRunge(Ball ball) {

if ((ball.y0 - ball.r) == 0 && ball.vy0 == 0) return;

double[] k = new double[4];

double[] l = new double[4];

k[0] = H * ball.vy0;

l[0] = H * (-G);

k[1] = H * (ball.vy0 + l[0] / 2.0);

l[1] = H * (-G);

k[2] = H * (ball.vy0 + l[1] / 2.0);

l[2] = H * (-G);

k[3] = H * (ball.vy0 + l[2]);

l[3] = H * (-G);

ball.y = ball.y0 + (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6.0;

ball.vy = ball.vy0 + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;

/* Rebound process */

if (ball.y - ball.r < 0) {

double diff = (ball.vy0 + Math.sqrt(ball.vy0*ball.vy0-2*G*(ball.r-ball.y0))) / G;

k[0] = diff * ball.vy0;

l[0] = diff * (-G);

k[1] = diff * (ball.vy0 + l[0] / 2.0);

l[1] = diff * (-G);

k[2] = diff * (ball.vy0 + l[1] / 2.0);

l[2] = diff * (-G);

k[3] = diff * (ball.vy0 + l[2]);

l[3] = diff * (-G);

ball.vy = ball.vy0 + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;

ball.y = ball.r;

ball.vy *= -E;

k[0] = (H - diff) * ball.vy;

l[0] = (H - diff) * (-G);

k[1] = (H - diff) * (ball.vy + l[0] / 2.0);

l[1] = (H - diff) * (-G);

k[2] = (H - diff) * (ball.vy + l[1] / 2.0);

l[2] = (H - diff) * (-G);

k[3] = (H - diff) * (ball.vy + l[2]);

l[3] = (H - diff) * (-G);

ball.y = ball.y + (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6.0;

ball.vy = ball.vy + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;

}

/* Update Values */

ball.y0 = ball.y;

ball.vy0 = ball.vy;

/* Stationary judgment */

if (ball.y - ball.r < 0) {

ball.vy = ball.vy0 = 0;

ball.y = ball.y0 = ball.r;

}

/* Drawing class */

class DrawComponent extends JPanel {

@Override

public void paintComponent(Graphics g) {

Graphics2D g2 = (Graphics2D) g;

/* Antialiasing */

g2.setRenderingHint(RenderingHints.KEY_ANTIALIASING,

RenderingHints.VALUE_ANTIALIAS_ON);

/* Fill the screen */

g2.setColor(Color.white);

g2.fillRect(0, 0, getWidth(), getHeight());

/* Draw information */

g2.setFont(new Font(null, Font.BOLD, 13));

g2.setColor(Color.black);

g2.drawString("Analysis", real.getX(), getHeight() - 50);

g2.drawString("Euler", euler.getX() + 10, getHeight() - 50);

g2.drawString("Runge-Kutta", runge.getX() - 10, getHeight() - 50);

g2.drawString("pos", 15, getHeight() - 32);

g2.drawString("vel", 15, getHeight() - 12);

DecimalFormat df0 = new DecimalFormat("00.000");

DecimalFormat df = new DecimalFormat("000.000000");

g2.setFont(new Font(null, Font.PLAIN, 13));

g2.drawString("t = " + df0.format(t), 10, 20);

g2.drawString("" + df.format(real.y - real.r), real.getX() - 5, getHeight() - 30);

g2.drawString("" + df.format(euler.y - euler.r), euler.getX() - 5, getHeight() - 30);

g2.drawString("" + df.format(runge.y - runge.r), runge.getX() - 5, getHeight() - 30);

g2.drawString("" + df.format(real.vy), real.getX() - 5, getHeight() - 10);

g2.drawString("" + df.format(euler.vy), euler.getX() - 5, getHeight() - 10);

g2.drawString("" + df.format(runge.vy), runge.getX() - 5, getHeight() - 10);

/* Coordinate transformation */

g2.translate(0, getHeight() - 1);

g2.scale(1, -1);

/* Draw the ground */

g2.drawLine(0, ground, getWidth(), ground);

/* Draw the balls */

g2.drawOval(real.getX(), real.getY(), real.getD(), real.getD());

g2.drawOval(euler.getX(), euler.getY(), euler.getD(), euler.getD());

g2.drawOval(runge.getX(), runge.getY(), runge.getD(), runge.getD());

/* Draw the ball line */

g2.setColor(Color.gray);

if (drawLineAna) {

g2.drawLine(0, real.getY(), getWidth(), real.getY());

}

if (drawLineEuler) {

g2.drawLine(0, euler.getY(), getWidth(), euler.getY());

}

if (drawLineRunge) {

g2.drawLine(0, runge.getY(), getWidth(), runge.getY());

}

/* Animation class */

class Animation extends Thread {

long delay = 1;

Animation(){

delay = (long)(H / 0.01);

if (delay == 0) delay = 1;

}

@Override

public void run() {

while (start) {

/* Advance time step */

t += H;

/* Update the ball position */

calcReal(real);

calcEuler(euler);

calcRunge(runge);

repaint();

try {

Thread.sleep(delay);

} catch (InterruptedException e) {

}

/* Constructor to set the GUI */

FreeFall(String title) {

/* Set the JFrame */

setTitle(title);

setDefaultCloseOperation(EXIT_ON_CLOSE);

setSize(640, 480);

setLocationRelativeTo(null);

setLayout(new FlowLayout());

/* Add the draw panel on the JFrame */

DrawComponent panelDraw = new DrawComponent();

panelDraw.setPreferredSize(new Dimension(500, 480));

panelDraw.setBorder(new EtchedBorder());

add(panelDraw);

/* Add the control panel on the JFrame */

JPanel panelControl = new JPanel();

panelControl.setPreferredSize(new Dimension(230, 480));

add(panelControl);

{

/* Add the parameter panel on the control panel */

JPanel panelParameter = new JPanel();

panelParameter.setPreferredSize(new Dimension(230, 270));

panelParameter.setBorder(new TitledBorder("Parameters"));

panelParameter.setLayout(new FlowLayout());

panelControl.add(panelParameter);

{

/* Add the initial position form on the parameter panel */

JLabel labelPositon = new JLabel("Initial position (y0)");

labelPositon.setPreferredSize(new Dimension(155, 30));

panelParameter.add(labelPositon);

spPosition = new JSpinner();

spPosition.setModel(new SpinnerNumberModel(300, 0, null, 1));

spPosition.setPreferredSize(new Dimension(50, 20));

panelParameter.add(spPosition);

/* Add the initial velocity form on the parameter panell */

JLabel labelVelocity = new JLabel("Initial velocity (vy0)");

labelVelocity.setPreferredSize(new Dimension(155, 30));

panelParameter.add(labelVelocity);

spVelocity = new JSpinner();

spVelocity.setModel(new SpinnerNumberModel(0d, null, null, 1d));

spVelocity.setPreferredSize(new Dimension(50, 20));

panelParameter.add(spVelocity);

/* Add the ball radious form on the parameter panel */

JLabel labelBall = new JLabel("Ball radious (r)");

labelBall.setPreferredSize(new Dimension(155, 30));

panelParameter.add(labelBall);

spBall = new JSpinner();

spBall.setModel(new SpinnerNumberModel(30, null, null, 1));

spBall.setPreferredSize(new Dimension(50, 20));

panelParameter.add(spBall);

/* Add the restitution coefficient form on the parameter panel */

JLabel labelE = new JLabel("Restitution coefficient (E)");

labelE.setPreferredSize(new Dimension(155, 30));

panelParameter.add(labelE);

spE = new JSpinner();

spE.setModel(new SpinnerNumberModel(1.0d, 0d, 1d, 0.1d));

spE.setPreferredSize(new Dimension(50, 20));

panelParameter.add(spE);

/* Add the gravitational acceleration form on the parameter panel */

JLabel labelG = new JLabel("Gravitational acceleration (G)");

labelG.setPreferredSize(new Dimension(155, 30));

panelParameter.add(labelG);

spG = new JSpinner();

spG.setModel(new SpinnerNumberModel(9.80665d, null, null, 0.1d));

spG.setPreferredSize(new Dimension(50, 20));

panelParameter.add(spG);

/* Add the time step form on the parameter panel */

JLabel labelH = new JLabel("Step size (H)");

labelH.setPreferredSize(new Dimension(155, 30));

panelParameter.add(labelH);

spH = new JSpinner();

spH.setModel(new SpinnerNumberModel(0.05d, 0d, null, 0.005d));

spH.setPreferredSize(new Dimension(50, 20));

panelParameter.add(spH);

JLabel indent = new JLabel();

indent.setPreferredSize(new Dimension(85, 20));

panelParameter.add(indent);

/* Add the apply button on the parameter panel */

JButton btnApply = new JButton("Apply");

btnApply.setPreferredSize(new Dimension(120, 30));

btnApply.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

E = (Double) spE.getValue();

G = (Double) spG.getValue();

H = (Double) spH.getValue();

/* Initialize the time */

t = 0;

tn = 0;

/* Initialize the ball position */

real.reset();

euler.reset();

runge.reset();

start = false;

repaint();

}

});

panelParameter.add(btnApply);

}

/* Add the view setting panel on the control panel */

JPanel panelView = new JPanel();

panelView.setPreferredSize(new Dimension(230, 100));

panelView.setBorder(new TitledBorder("View setting"));

panelView.setLayout(new FlowLayout());

panelControl.add(panelView);

{

JCheckBox checkAna = new JCheckBox("Draw the analysis ball line");

checkAna.setPreferredSize(new Dimension(200, 20));

checkAna.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

drawLineAna = !drawLineAna;

repaint();

}

});

panelView.add(checkAna);

JCheckBox checkEuler = new JCheckBox("Draw the Euler ball line");

checkEuler.setPreferredSize(new Dimension(200, 20));

checkEuler.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

drawLineEuler = !drawLineEuler;

repaint();

}

});

panelView.add(checkEuler);

JCheckBox checkRunge = new JCheckBox("Draw the Runge-Kutta ball line");

checkRunge.setPreferredSize(new Dimension(200, 20));

checkRunge.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

drawLineRunge = !drawLineRunge;

repaint();

}

});

panelView.add(checkRunge);

}

/* Add the simulation panel on the control panel */

JPanel panelSimulation = new JPanel();

panelSimulation.setPreferredSize(new Dimension(230, 95));

panelSimulation.setBorder(new TitledBorder("Simulation"));

panelSimulation.setLayout(new FlowLayout());

panelControl.add(panelSimulation);

{

/* Add the start button on the simulation panel */

JButton btnStart = new JButton("Start");

btnStart.setPreferredSize(new Dimension(100, 30));

btnStart.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

if (!start) {

start = true;

Animation anime = new Animation();

anime.start();

}

});

panelSimulation.add(btnStart);

/* Add the stop button on the simulation panel */

JButton btnStop = new JButton("Stop");

btnStop.setPreferredSize(new Dimension(100, 30));

btnStop.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

start = false;

}

});

panelSimulation.add(btnStop);

/* Add the reset button on the simulation panel */

JButton btnReset = new JButton("Reset");

btnReset.setPreferredSize(new Dimension(100, 30));

btnReset.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

E = (Double) spE.getValue();

G = (Double) spG.getValue();

H = (Double) spH.getValue();

/* Initialize the time */

t = 0;

tn = 0;

/* Initialize the ball position */

real.reset();

euler.reset();

runge.reset();

start = false;

repaint();

}

});

panelSimulation.add(btnReset);

/* Add the step by step button on the simulation panel */

JButton btnStep = new JButton("Step by step");

btnStep.setPreferredSize(new Dimension(100, 30));

btnStep.addActionListener(new ActionListener() {

public void actionPerformed(ActionEvent e) {

/* Advance time step */

t += H;

/* Update the ball position */

calcReal(real);

calcEuler(euler);

calcRunge(runge);

repaint();

}

});

panelSimulation.add(btnStep);

}

pack();

/* Initialize the balls */

real = new Ball(100);

euler = new Ball(250);

runge = new Ball(400);

}

public static void main(String[] args) {

/* Set the Windows look and feel */

try {

javax.swing.UIManager

.setLookAndFeel("com.sun.java.swing.plaf.windows.WindowsLookAndFeel");

} catch (Exception e) {

}

/* Create and display the form */

FreeFall frame = new FreeFall("FreeFall");

frame.setVisible(true);

}

シミュレーション結果は次のようになりました。

尚、反発係数は0.9、時間ステップは0.05、0.5としています。

時間ステップを増加させると、Euler法では誤差が大きくなっていくことがわかります。

一方で、Runge-Kutta法では全くといっていいほど誤差が見られません。

数学的に同一の処理をしているので当然と言えば当然なのですが、Runge-Kutta法の実装には細かいミス多発で結構手こずりました。

また時間があったら、もっと複雑な微分方程式を解いてみて、性能検証してみようと思います。疲れた。

自由落下をオイラー法とルンゲクッタ法で解く

1. 自由落下の運動方程式

2. 解析的手法

3. 解析的手法（跳ね返り）

4. オイラー法

5. オイラー法（跳ね返り）

6. ルンゲクッタ法

7. ルンゲクッタ法（跳ね返り）

8. C言語による実装（跳ね返り無し）

9. Javaによる実装（跳ね返りあり）

コメントを残すコメントをキャンセル

カテゴリー

アーカイブ

人気記事

アクセスカウンター

自由落下をオイラー法とルンゲクッタ法で解く

1. 自由落下の運動方程式

2. 解析的手法

3. 解析的手法（跳ね返り）

4. オイラー法

5. オイラー法（跳ね返り）

6. ルンゲクッタ法

7. ルンゲクッタ法（跳ね返り）

8. C言語による実装（跳ね返り無し）

9. Javaによる実装（跳ね返りあり）

コメントを残す コメントをキャンセル

カテゴリー

アーカイブ

人気記事

アクセスカウンター

コメントを残すコメントをキャンセル