微分方程式を解くルンゲクッタ法を試してみます。
今回は
という微分方程式を解いて、空気抵抗がある自由落下の速度を求めてみます。g=9.8、k=0.01としてます。
空気抵抗により、速度が一定になることがわかります。
プログラム中、funcメソッドが求める微分方程式の右辺になります。
public class RungeKutta { /** 微分方程式 dy/dx = f(x,y) の f(x,y)*/ public static double func(double x, double y){ //重力加速度9.8、空気抵抗係数0.01のときの落下速度 return 9.8 - .01 * y * y; } public static void main(String[] args){ //描画準備 BufferedImage img = new BufferedImage( 400, 300, BufferedImage.TYPE_INT_RGB); Graphics g = img.getGraphics(); g.setColor(Color.WHITE); g.fillRect(0, 0, 400, 300); g.setColor(Color.BLUE); //計算範囲 double x0 = 0; double xm = 20; //初期条件の設定 double x = x0; double y = 0; //刻み幅 double h = .1; for(; x < xm; x += h){ //点を打つ int px = (int) ((x - x0) * 300 / (xm - x0)); int py = (int)(y / 30 * 100); g.fillOval(10 + px, 200 - py, 2, 2); //ルンゲクッタ法で次のyを求める double k1 = h * func(x, y); double k2 = h * func(x + h / 2, y + k1 / 2); double k3 = h * func(x + h / 2, y + k2 / 2); double k4 = h * func(x + h, y + k3); double k = (k1 + 2 * k2 + 2 * k3 + k4) / 6; y += k; } //表示 JFrame f = new JFrame("ルンゲクッタで1階微分方程式"); f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); f.setSize(400, 300); f.add(new JLabel(new ImageIcon(img))); f.setVisible(true); } }