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

微分方程式を解くルンゲクッタ法を試してみます。
今回は
\frac{dv}{dt}=g-kv^2
という微分方程式を解いて、空気抵抗がある自由落下の速度を求めてみます。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);
    }
}