ルンゲクッタ法で2階連立微分方程式を解く

java 数値計算

今回は、空気抵抗があるときの物体の軌跡を描いてみました。下のプログラムでは、実際にはアニメーションさせてます。

※追記 2023/4/30 動画を投稿しました

2階以上の微分方程式は、
\frac{d^2y}{dx^2}=f(x, y, \frac{dy}{dx})
を、
w=\frac{dy}{dx}
と置くことで
\left\{w=\frac{dy}{dx}\atop\frac{d^2y}{dx^2}=f(x, y, w)\right.
という1階連立微分方程式に分解して、ルンゲクッタ法を使います。

このとき物体の位置x,yの微分方程式は、空気抵抗の係数をT、重力加速度をgとすると
\left\{ {\frac{d^2x}{dt^2}=-T\frac{dx}{dt}\sqrt{(\frac{dx}{dt})^2+(\frac{dy}{dt})^2}} \atop {\frac{d^2y}{dt^2}=-g-T\frac{dy}{dt}\sqrt{(\frac{dx}{dt})^2+(\frac{dy}{dt})^2}}\right.
となります。
ここで
w=\frac{dx}{dt}, v=\frac{dy}{dt}
と置くことで
\left\{ \frac{dw}{dt}=-Tw\sqrt{w^2+v^2} \atop \frac{dv}{dt}=-g-Tv\sqrt{w^2+v^2}\right.
のような１階微分方程式に変形して計算を行います。
プログラム中では、funcw・funcvメソッドが変形後の式にあたります。

こうやって2階微分方程式が解けるようになると、いろいろな物理現象がシミュレートできるようになって楽しいです。

import java.awt.*;
import java.awt.image.BufferedImage;
import javax.swing.*;

public class RungeKutta2 {
  static double gt = 9.8;//重力加速度
  static double tt = 0.008;//空気抵抗の係数
  
  private static double funcv(double t, double v, double w, double x, double y) {
    return -gt - tt * v * Math.sqrt(w * w + v * v);
  }
  private static double funcw(double t, double v, double w, double x, double y) {
    return -tt * w * Math.sqrt(w * w + v * v);
  }
  private static double funcx(double t, double v, double w, double x, double y) {
    return w;
  }
  private static double funcy(double t, double v, double w, double x, double y) {
    return v;
  }

  public static void main(String[] args) {
    //表示準備
    BufferedImage bi = new BufferedImage(400, 300, BufferedImage.TYPE_INT_RGB);
    final Graphics g = bi.getGraphics();

    //画面に表示
    JFrame f = new JFrame("ルンゲクッタ２階連立微分方程式");
    f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
    f.setSize(400, 300);
    final JLabel label = new JLabel(new ImageIcon(bi));
    f.add(label);
    f.setVisible(true);

    new Thread() {
      @Override public void run() {
        for(;;){
          try { Thread.sleep(300);
          } catch (InterruptedException e) { }
          //初期条件
          double h = .1;//刻み幅
          double x = 0;//初期位置
          double y = 0;
          double sp = Math.random() * 100 + 20;//初速
          double th = (Math.random() * 3 / 4 + 1. / 8) * Math.PI;//角度
          double v = sp * Math.sin(th);
          double w = sp * Math.cos(th);

          for (double t = 0; t < 20; t += h) {
            g.setColor(Color.WHITE);
            g.fillRect(0, 0, 400, 300);
            g.setColor(Color.RED);
            g.fillOval((int)(10 + x * 1.5), (int) (250 - y * 1.5), 7, 7);

            double vk1 = h * funcv(t, v, w, x, y);
            double wk1 = h * funcw(t, v, w, x, y);
            double xk1 = h * funcx(t, v, w, x, y);
            double yk1 = h * funcy(t, v, w, x, y);
            double vk2 = h * funcv(t + h/2, v + vk1/2, w + wk1/2, x + xk1/2, y + yk1/2);
            double wk2 = h * funcw(t + h/2, v + vk1/2, w + wk1/2, x + xk1/2, y + yk1/2);
            double xk2 = h * funcx(t + h/2, v + vk1/2, w + wk1/2, x + xk1/2, y + yk1/2);
            double yk2 = h * funcy(t + h/2, v + vk1/2, w + wk1/2, x + xk1/2, y + yk1/2);
            double vk3 = h * funcv(t + h/2, v + vk2/2, w + wk2/2, x + xk2/2, y + yk2/2);
            double wk3 = h * funcw(t + h/2, v + vk2/2, w + wk2/2, x + xk2/2, y + yk2/2);
            double xk3 = h * funcx(t + h/2, v + vk2/2, w + wk2/2, x + xk2/2, y + yk2/2);
            double yk3 = h * funcy(t + h/2, v + vk2/2, w + wk2/2, x + xk2/2, y + yk2/2);
            double vk4 = h * funcv(t + h, v + vk3, w + wk3, x + xk3, y + yk3);
            double wk4 = h * funcw(t + h, v + vk3, w + wk3, x + xk3, y + yk3);
            double xk4 = h * funcx(t + h, v + vk3, w + wk3, x + xk3, y + yk3);
            double yk4 = h * funcy(t + h, v + vk3, w + wk3, x + xk3, y + yk3);
            double vk = (vk1 + 2 * vk2 + 2 * vk3 + vk4) / 6;
            double wk = (wk1 + 2 * wk2 + 2 * wk3 + wk4) / 6;
            double xk = (xk1 + 2 * xk2 + 2 * xk3 + xk4) / 6;
            double yk = (yk1 + 2 * yk2 + 2 * yk3 + yk4) / 6;
            v += vk;
            w += wk;
            x += xk;
            y += yk;
            if (y < 0) v = -v; //反射させてみる
            label.repaint();
            try { Thread.sleep(20);
            } catch (InterruptedException e) { }
          }
        }
      }
    }.start();
  }
}