1,常備分方程式 (17KB) | Euler法 Runge−Kutta法 レポート |
2,偏微分方程式 | Euler法 Crank-Nicholson法 レポート |
3,確率論 (10KB) | モンテカルロ法 レポート |
偏微分方程式は、常備分方程式の解法と比べて解きにくいものになっています。
それは、それぞれの偏微分方程式にあわせて、個々に式を作る必要があるからです。
ここでは、次の拡散方程式を解くことにします。
∂u(r,t) = D ∇2u ∂t
ここで、uは濃度、Dは拡散係数です。
今の方程式は次のように変形できます。
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で評価しているので
誤差が大きくなってしまいます。
Euler法の誤差の原因である時間のずれを補正したのがCrank-Nicholson法です。
これで精度が上がったのですが、その代償として式が複雑になっています。
Crank-Nicholson法の基本式が次の通りです。
u(x,t+Δt)−u(x,t) = 1{∂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で初期状態に)
以下の解答は、先ほどのプログラムを改造することによって制作できます。
先程の関数は、
for (int i=0; i<=N+1; i++) u[i]=0.0; u[20]=1.0/delx;
でした、この部分を変えることでいろいろな関数についてシミュレーションできます。
吸収壁はu[0]=u[N+1]=0 とすることによって実現できます。
関数「boundary」wp次のようにします。
public void boundary(){ u[0]=u[N+1]=0;}
これで終了です。
湧きだしの方程式は次のようになります。
∂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行目で湧き出す場所だけ湧き出しの項を加えます。
場所によって拡散係数が違うので、Dを配列にします。
また、Dの関数であるdel変数も配列にする必要があります。
Dが変数の場合の運動方程式は次のようになります。
∂u = ∇{D(x) ∇u(x,t)} ∂t
以下が実行結果です。(Clickで動作/停止 DoubleClickで初期状態に)
バグ?のためにDの値が変わる境界地で値が不安定になっています。
・・省略・・