波動方程式を解いてみるよ

機械学習は飽きてきたので、今日は波動方程式を解いてみた。
弦をひっぱって手を離したらどういう運動をするかというシミュレーション。右下方向が時間のたつ向き

つっても、ほぼ、本に載ってる通りなんだけど。初期値が違うくらい。
ちなみに数値計算の参考にしてるのはこの本。

数値計算のはなし―パソコンを使って解こう

ほかの本は、行列の計算とかフーリエ解析とかまでやってるけど、この本は明確に偏微分方程式を解くための本で、その過程として連立1次方程式解いて積分して、という感じになっているので、シミュレーションしたいときには最短距離になる感じ。説明も小難しくない。
サンプルがBASICだったりして最初は萎えたけど、Cで書かれてるよりも簡単だし、Javaでサンプル作ってる本だと変に「オブジェクト指向」に凝っててプログラムが見にくかったりするので、実はとてもよいです。

なので、これはほとんど、その本のサンプルのJava版という感じ。

//wave.java
import java.awt.*;
import java.awt.geom.Path2D;
import java.awt.image.BufferedImage;
import javax.swing.*;

public class Wave {
    public static void main(String[] args){
        final int TDIV = 120;//時間の分割数
        final int XDIV = 20;//座標の分割数
        final double TLIM = 3;//時間
        double[][] data = new double[TDIV][XDIV];
        double a = TLIM * XDIV / TDIV;
        a = a * a;

        //初期値
        //最初の形
        for(int x = 0; x < XDIV / 2; ++x){
            data[0][x] = 2. * x / XDIV;
            data[0][XDIV - x - 1] = 2. * x / XDIV;
        }
        for(int x = 1; x < XDIV - 1; ++x){
            data[1][x] = data[0][x] + 
                    a * (data[0][x + 1] - 2 * data[0][x] + data[0][x - 1]) / 2;
        }
        //偏微分方程式を求める
        for(int i = 0; i < 1000; ++i){
            double[][] pre = new double[TDIV][XDIV];
            //値の更新
            for(int t = 1; t < TDIV - 1; ++t){
                for(int x = 1; x < XDIV - 1; ++x){
                    pre[t + 1][x] = data[t + 1][x];
                    data[t + 1][x] = 2 * data[t][x] - data[t - 1][x] + 
                            a * (data[t][x + 1] - 2 * data[t][x] + data[t][x - 1]);
                }
            }
            //変位の計算
            double p = 0;
            double sum = 0;
            for(int t = 1; t < TDIV - 1; ++t){
                for(int x = 1; x < XDIV - 1; ++x){
                    p += Math.abs(pre[t][x] - data[t][x]);
                    sum += Math.abs(data[t][x]);
                }
            }
            //変位が少なければ収束したことにする
            if(p < sum * 0.0001) break;
        }
        //グラフ作成
        Image img = new BufferedImage(400, 300, BufferedImage.TYPE_INT_RGB);
        Graphics2D g = (Graphics2D) img.getGraphics();
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, 400, 300);
        g.setColor(Color.BLACK);
        for(int t = 0; t < TDIV; ++t){
            int tx = t * 200 / TDIV;
            Path2D p = new Path2D.Double();
            p.moveTo(tx + 30, tx / 2- data[t][0] + 160);
            for(int x = 1; x < XDIV; ++x){
                p.lineTo(x * 150 / XDIV + tx  + 30, 
                        tx / 2 - x * 150 / XDIV  + 160 - data[t][x] * 50);
            }
            g.draw(p);
        }
        g.dispose();
        //画面表示
        JFrame f = new JFrame("波動方程式");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(400, 300);
        f.add(new JLabel(new ImageIcon(img)));
        f.setVisible(true);
    }
}