移流方程式の計算で数値計算の問題点をみる

数値計算でどんな微分方程式でもいい感じに計算できるかというと、そういうわけはもちろんなく。

という移流方程式をいくつかの計算をやってみて、理論どおりにならないことを見てみました。

左辺の第２項の偏微分について、いくつかのやりかたをしています。
一番上が理論値です。形は変わらずにただ移動するだけになるはずです。
２番目は後退差分。どんどん形がつぶれていくのがわかります。
３番目は、後退差分より精度がいいはずの中心差分。なんか、ひどいノイズが出てます。
４番目は、ラックス-ベンドロフ法という方法で中心差分に補正をいれています。反射してきたノイズがちょっと抑えられています。

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

public class Advect {
    public static void main(String[] args){
        final BufferedImage img = new BufferedImage(400, 300, BufferedImage.TYPE_INT_RGB);
        final JLabel label = new JLabel(new ImageIcon(img));
        
        JFrame f = new JFrame("移流");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.add(label);
        f.setSize(400, 350);
        f.setVisible(true);
        
        final double[] data = new double[50];
        final double dt = 0.01;
        final double dx = 1;
        final double c = 1;
        final String[] names = {"理論値", "後退差分", "中心差分", "ラックス-ベンドロフ"};
        //初期値
        for(int i = 1; i < 5; ++i) data[i] = 1;
        
        new Thread(){
            @Override
            public void run() {
                double t = 0;
                double d[][] = new double[4][];
                double[][] next = new double[4][data.length];
                for(int i = 0; i < d.length; ++i) d[i] = data.clone();
        
                for(;;){
                    //計算
                    for(int n = 0; n < 30; ++n){
                        t += dt;
                        for(int x = 1; x < data.length - 1; ++x){
                            //理論値
                            int ii = (int) (x - t);
                            if(ii >= 0 && ii < data.length){
                                next[0][x] = data[ii];
                            }
                            //後退差分
                            next[1][x] = d[1][x] - c * dt / dx * (d[1][x] - d[1][x - 1]);
                            //中心差分
                            next[2][x] = d[2][x] - c * dt / dx / 2 * (d[2][x + 1] - d[2][x - 1]);
                            //ラックス-ベンドロフ法による補正
                            next[3][x] = d[3][x] - c * dt / dx / 2 * (d[3][x + 1] - d[3][x - 1]) + 
                                    c * c * dt * dt / dx / dx / 2 * (d[3][x + 1] - 2 * d[3][x] + d[3][x - 1]);
                        }
                        d = next.clone();
                    }
                    
                    //表示
                    Graphics2D g = (Graphics2D) img.getGraphics();
                    g.setColor(Color.WHITE);
                    g.fillRect(0, 0, 400, 300);
                    for(int n = 0; n < d.length; ++n){
                        GeneralPath gp = new GeneralPath();
                        for(int i = 0; i < d[n].length; ++i){
                            double x = i * 7 + 20;
                            double y = 40 + n * 75 - d[n][i] * 35;
                            if(i == 0){
                                gp.moveTo(x, y);
                            }else{
                                gp.lineTo(x, y);
                            }
                        }
                        g.setColor(Color.BLACK);
                        g.draw(gp);
                        g.setColor(Color.GRAY);
                        g.drawString(names[n], 10, 10 + n * 75);
                    }
                    g.dispose();
                    
                    label.repaint();
                    try{
                        Thread.sleep(10);
                    }catch(InterruptedException e){}
                }
            }
        }.start();
                
    }
}