目次
以前こちらで雑に考察した自由落下。
かれこれ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」に書き出すプログラムとなっています。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 |
#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コンポーネント書きました。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 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法の実装には細かいミス多発で結構手こずりました。
また時間があったら、もっと複雑な微分方程式を解いてみて、性能検証してみようと思います。疲れた。