目次
以前こちらで雑に考察した自由落下。
かれこれ3年以上前ということで感慨深いものがありましたが、
中途半端に終わっているので決着をつけようと思います。
1. 自由落下の運動方程式
真空中を自由落下する物体には重力のみがはたらくので、次の運動方程式が成立します。
$$\displaystyle \frac{d^2}{dt^2}y = -g \ \ \ \ (1)$$
この式をオイラー法、ルンゲクッタ法で位置 $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;
}
自由落下だけということもあり、割と短いソースコードになりました。
※コードをコピーする際は、タブ文字がスペースに置換されない単純ビュー（画面ダブルクリック）にすることをお勧めします。
書き出された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);
}
}
シミュレーション結果は次のようになりました。
尚、反発係数は0.9、時間ステップは0.05、0.5としています。
時間ステップを増加させると、Euler法では誤差が大きくなっていくことがわかります。
一方で、Runge-Kutta法では全くといっていいほど誤差が見られません。
数学的に同一の処理をしているので当然と言えば当然なのですが、Runge-Kutta法の実装には細かいミス多発で結構手こずりました。
また時間があったら、もっと複雑な微分方程式を解いてみて、性能検証してみようと思います。疲れた。
こんにちは。質問なのですが、
跳ね返りなしのルンゲクッタ法の公式が初めて見る形式で、理論的な部分を知りたいのですが参考になる記事や書籍などございますか？