タイトル

計算物理学

 偏微分方程式

タイトル

1,常備分方程式 (17KB) Euler法 Runge−Kutta法 レポート
2,偏微分方程式 Euler法 Crank-Nicholson法 レポート
3,確率論 (10KB) モンテカルロ法 レポート

葉

2,偏微分方程式の解法

偏微分方程式は、常備分方程式の解法と比べて解きにくいものになっています。
それは、それぞれの偏微分方程式にあわせて、個々に式を作る必要があるからです。
ここでは、次の拡散方程式を解くことにします。


∂u(r,t) = D ∇2u
∂t

ここで、uは濃度、Dは拡散係数です。

ボーダー

Euler法

今の方程式は次のように変形できます。


u(x,t+Δt)−u(x,t) = ∂u2(x,t)
    Δt

u(x,t+Δt) = u(x,t)+D∂2u(x,t)Δt

さらに

∂2u(x,t) = u(x−Δ) + u(x+Δ) − 2u(x)
              Δ2
なので
u(x,t+Δt)=u(x,t) + DΔt u(x−Δ) + u(x+Δ) − 2u(x)
                     Δ2

このように比較的簡単に書き表すことができます。
これを「Euler法」といいます。この式は、未知量t+Δtが右辺だけで解ける
形になっているので「陽的差分(explicit)」といいます。

しかし、Euler法は空間微分をtで、時間微分をt+Δt/2で評価しているので
誤差が大きくなってしまいます。

ボーダー

Crank-Nicholson法

Euler法の誤差の原因である時間のずれを補正したのがCrank-Nicholson法です。
これで精度が上がったのですが、その代償として式が複雑になっています。
Crank-Nicholson法の基本式が次の通りです。


u(x,t+Δt)−u(x,t) = {∂u2(x,t) + ∂u2(x,t+Δt)}
    Δt       2

Euler法と違うのは右辺がtでの微分ではなく
t+Δt/2での微分になっているということです。

右辺と左辺の時間は一致しているため精度が上がっています。
しかし、Euler法のように左辺に未知量数について解けた式ではなく
連立方程式を使って解くことになります。
このようなスキームを「陰的差分(implicit)」といいます。


問題
Crank-Nicholson法を用いて拡散方程式をシミュレーションせよ

解答
以下のプログラムで実行できます。

なお、重要でない部分は省略しています。

import java.applet.*;
import java.awt.*;

public class j30b1 extends Applet implements Runnable
{
    private Thread m_j30b1 = null;
    boolean move=false;
    int N=64;
    double []u=new double[N+2];
    double [][]a=new double[N+1][N+1];
    double []x=new double[N+1];
    double []b=new double[N+1];
    double D,del,delx,delt;

    public j30b1(){}

    public boolean gauss() 
    {
        double [][]work=new double[N+1][N+2];
        int i,j,k,l;
        for (i=1; i<=N; i++){
            for (j=1; j<=N; j++)
                work[i][j]=a[i][j];
            work[i][N+1]=b[i];
        }
        //前進部分
        for (k=1; k<=N-1; k++){
            if (Math.abs(work[k][k])<1.0e-6){return false;}
            for(j=k+1; j<=N+1; j++){
                for(i=k+1; i<=N; i++)
                    work[i][j]-=(work[i][k]/work[k][k])*work[k][j];
            }
        }
        //後進部分
        for (k=N; k>=1; k--){
            x[k]=work[k][N+1];
            for (l=N; l>=k+1; l--)
                x[k]=x[k]-x[l]*work[k][l];
            x[k]=x[k]/work[k][k];
            }
        return true;
    }

    public void boundary() 
    {
        u[1]=u[N]=0.0;
    }

    public void setting() 
    {
        D=1.0;delx=0.5;delt=0.2;
        del=D*delt/(2.0*delx*delx);
        //********* u ************
        for (int i=0; i<=N+1; i++)
            u[i]=0.0;
        u[20]=1.0/delx;

        //********* a ************
        for (int i=2;i<=N-1; i++){
            a[i][i]=1.0+2.0*del;
            a[i][i-1]=-del;
            a[i][i+1]=-del;
        }
        a[1][1]=a[N][N]=1.0+2.0*del;
        a[1][N]=a[1][2]=a[N][N-1]=a[N][1]=-del;
    }

    public void init()
    {
        resize(320, 240);
        setBackground(Color.white);
        setting();
    }

    public void paint(Graphics g)
    {
        g.drawLine(0,239,320,239);
        g.setColor(Color.red);
        for (int i=1; i<N; i++)
            g.drawLine(i*5,239-(int)(u[i]*100),(i+1)*5,239-(int)(u[i+1]*100));
    }

    public void start()
    {
        if (m_j30b1 == null)
        {
            m_j30b1 = new Thread(this);
            m_j30b1.start();
        }
    }

    public void stop()
    {
        if (m_j30b1 != null)
        {
            m_j30b1.stop();
            m_j30b1 = null;
        }
    }

    public void run()
    {
        int i;
        while (true){
            try{Thread.sleep(30);}catch (InterruptedException e){stop();}
            if (move){
                for (i=1; i<=N; i++)
                    b[i]=(1-2*del)*u[i]+del*(u[i-1]+u[i+1]);
                b[60]+=0.1*delt*1.0/delx;
                if (!gauss()){return;}
                for (i=1; i<=N; i++)
                    u[i]=x[i];
                boundary();
                repaint();
            }
        }
    }

    public boolean mouseDown(Event evt, int x, int y)
    {
        move=!move;

        if (evt.clickCount==2){
            setting();
            move=true;}

        if (move){
            m_j30b1.resume();
        }else{
            m_j30b1.suspend();
        }
        return true;
    }
}

以下が実行結果です。(Clickで動作/停止 DoubleClickで初期状態に)

ボーダー

レポート

●課題

  1. 様々な関数でシミュレーションし、グリーン関数の立場から解析せよ
  2. 境界条件が吸収壁の場合、どのような時間発展がみられるか?
  3. 境界条件が吸収壁で湧きだしがあるとき、系が定常状態に近づくことを確かめよ
  4. 拡散係数がD(x)でxの位置によって違う系は、一定の場合とどのように違うか?
  5. その他、いろいろな状況を考えてみよ

以下の解答は、先ほどのプログラムを改造することによって制作できます。

●解答1

先程の関数は、


        for (int i=0; i<=N+1; i++)
            u[i]=0.0;
        u[20]=1.0/delx;

でした、この部分を変えることでいろいろな関数についてシミュレーションできます。

●解答2

吸収壁はu[0]=u[N+1]=0 とすることによって実現できます。
関数「boundary」wp次のようにします。

  public void boundary(){
  u[0]=u[N+1]=0;}

これで終了です。

●解答3

湧きだしの方程式は次のようになります。


∂u = D ∂2u + αδ(x-xo)
∂t     ∂x2

右辺第2項が湧き出しの項です。αはわきだしの強さ、xoは湧き出し地点です。


for (i=1; i<=N; i++)
    b[i]=(1-2*del)*u[i]+del*(u[i-1]+u[i+1]);
b[60]+=0.1*delt*1.0/delx;

1,2行は今までと変わりません、3行目で湧き出す場所だけ湧き出しの項を加えます。

●解答4

場所によって拡散係数が違うので、Dを配列にします。
また、Dの関数であるdel変数も配列にする必要があります。

Dが変数の場合の運動方程式は次のようになります。


∂u = ∇{D(x) ∇u(x,t)}
∂t

以下が実行結果です。(Clickで動作/停止 DoubleClickで初期状態に)

バグ?のためにDの値が変わる境界地で値が不安定になっています。

●解答5

・・省略・・

葉

上上へ

戻るTopに戻る