ラプラシアンの計算をしてラプラス方程式を解く

温度の拡散などをあらわすラプラス方程式
$\nabla^2f=0$
を解いてみました。
ここでは、左右の辺が1℃、下の辺が0℃になるとき、板の温度分布がどうなるかという計算になっています。

ラプラシアンというのは $\nabla^2$ のことです。
普通の偏微分であらわすと
$\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}$
となります。
この形の偏微分方程式を楕円形方程式といいます。

これを、プログラムで計算するときは、微分を差分にして引き算することになります。で、2階微分は次のように差分にできるんでした。
$\frac{d^2f}{dx^2}=\frac{f[x-1]-2f[x]+f[x+1]}{\Delta x^2}$
そうすると、上のラプラシアンは結局、次のように差分に変換できることになります。
$\nabla^2f=\frac{f[x-1][y]-2f[x][y]+f[x+1][y]}{\Delta x^2}+\frac{f[x][y-1]-2f[x][y]+f[x][y+1]}{\Delta y^2}$
ここで、ΔxとΔyはそれぞれx・yの刻み幅ですが、ΔxとΔyが等しいとすると
$\nabla^2f=\frac{f[x-1][y]-2f[x][y]+f[x+1][y] +f[x][y-1]-2f[x][y]+f[x][y+1]}{\Delta x^2}$
$=\frac{f[x-1][y]+f[x+1][y] +f[x][y-1]+f[x][y+1]-4f[x][y]}{\Delta x^2}$
となります。

ラプラス方程式は、次のような形で、熱の拡散などを表します。
$\nabla^2f=0$
これを差分式であらわすと
$\frac{f[x-1][y]+f[x+1][y] +f[x][y-1]+f[x][y+1]-4f[x][y]}{\Delta x}=0$
ということになって、f[x][y]の値は次のようになります。
$f[x][y]=\frac{f[x-1][y]+f[x+1][y]+f[x][y-1]+f[x][y+1]}{4}$
要するに、その場の温度は、上下左右の温度を足して4で割ったものという、簡単な式になるわけです。
すべての点についてこの方程式がみたされるようにするわけで、連立方程式を解くということになるんですが、実際のところは繰り返し計算し続ければいいということになります。

プログラムはつぎのようになりました。ここで連立方程式を解くとき、値の変移を加速させてます。これをSOR法といいます。

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

public class Laplas {
    public static void main(String[] args){
        //計算
        final int w = 20;
        final int h = 20;
        double[][] p = new double[w][h];
        //境界条件
        for(int y = 0; y < h; ++y){
            p[0][y] = 1;
            p[w - 1][y] = 1;
        }
        
        //連立方程式を解く
        for(int n = 0; n < 1000; ++n){
            double sum = 0;
            for(int x = 1; x < w - 1; ++x){
                for(int y = 1; y < h - 1; ++y){
                    double old = p[x][y];
                    double ne = (p[x - 1][y] + p[x + 1][y] + p[x][y - 1] + p[x][y + 1]) / 4;
                    double dif = (ne - old) * 1.7;//変移を加速させる
                    p[x][y] += dif;
                    sum += dif;
                }
            }
            if(sum < 0.0001) {
                System.out.println(n);
                break;
            }
        }
        //描画
        final int KX = 10;
        final int KY = 5;
        final int PD = 100;//データの長さ
        final int OX = 10;//表示開始位置
        final int OY = 200;
        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 x = 0; x < w; ++x){
            for(int y = 0; y < h; ++y){
                int px = x * KX + y * KX + OX;
                int py = -x * KY + y * KY - (int)(p[x][y] * PD) + OY;
                //横線
                if(x < w - 1){
                    int pxx = (x + 1) * KX + y * KX + OX;
                    int pxy = -(x + 1) * KY + y * KY - (int)(p[x + 1][y] * PD) + OY;
                    g.drawLine(px, py, pxx, pxy);
                }
                //縦線
                if(y < h - 1){
                    int pyx = x * KX + (y + 1) * KX + OX;
                    int pyy = -x * KY + (y + 1) * KY - (int)(p[x][y + 1] * PD) + OY;
                    g.drawLine(px, py, pyx, pyy);
                }
            }
        }
        
        //表示
        JFrame f = new JFrame("ラプラス方程式");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(430, 350);
        JLabel label = new JLabel(new ImageIcon(img));
        f.add(label);
        f.setVisible(true);
        
        
    }
}