| １，常備分方程式 (17KB) | Ｅｕｌｅｒ法 Ｒｕｎｇｅ−Ｋｕｔｔａ法 レポート |
| ２，偏微分方程式 | Ｅｕｌｅｒ法 Crank-Nicholson法 レポート |
| ３，確率論 (10KB) | モンテカルロ法 レポート |
偏微分方程式は、常備分方程式の解法と比べて解きにくいものになっています。
それは、それぞれの偏微分方程式にあわせて、個々に式を作る必要があるからです。
ここでは、次の拡散方程式を解くことにします。
∂ｕ(r,t) ＝ Ｄ ∇2ｕ ∂ｔ
ここで、ｕは濃度、Ｄは拡散係数です。
今の方程式は次のように変形できます。
ｕ(x,t+Δt)−ｕ(x,t) ＝ ∂ｕ2(x,t) Δｔ ｕ(x,t+Δt) ＝ ｕ(x,t)＋Ｄ∂2ｕ(x,t)Δt さらに ∂2ｕ(x,t) ＝ ｕ(ｘ−Δ) ＋ ｕ(ｘ＋Δ) − ２ｕ(ｘ) Δ2 なので ｕ(x,t+Δt)＝ｕ(x,t) ＋ ＤΔt ｕ(ｘ−Δ) ＋ ｕ(ｘ＋Δ) − ２ｕ(ｘ) Δ2
このように比較的簡単に書き表すことができます。
これを「Ｅｕｌｅｒ法」といいます。この式は、未知量t+Δtが右辺だけで解ける
形になっているので「陽的差分(explicit)」といいます。
しかし、Euler法は空間微分をｔで、時間微分をｔ＋Δt/2で評価しているので
誤差が大きくなってしまいます。
Euler法の誤差の原因である時間のずれを補正したのがCrank-Nicholson法です。
これで精度が上がったのですが、その代償として式が複雑になっています。
Crank-Nicholson法の基本式が次の通りです。
ｕ(x,t+Δt)−ｕ(x,t) ＝ １｛∂ｕ2(x,t) ＋ ∂ｕ2(x,t+Δt)｝ Δｔ ２
Euler法と違うのは右辺がｔでの微分ではなく
ｔ＋Δｔ／２での微分になっているということです。
右辺と左辺の時間は一致しているため精度が上がっています。
しかし、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;
でした、この部分を変えることでいろいろな関数についてシミュレーションできます。
吸収壁はｕ[０]＝ｕ[Ｎ＋１]＝０ とすることによって実現できます。
関数「boundary」ｗｐ次のようにします。
public void boundary(){
u[0]=u[N+1]=0;}
これで終了です。
湧きだしの方程式は次のようになります。
∂ｕ ＝ Ｄ ∂2ｕ ＋ αδ(x-xo) ∂ｔ ∂ｘ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;
１，２行は今までと変わりません、３行目で湧き出す場所だけ湧き出しの項を加えます。
場所によって拡散係数が違うので、Ｄを配列にします。
また、Ｄの関数であるdel変数も配列にする必要があります。
Ｄが変数の場合の運動方程式は次のようになります。
∂ｕ ＝ ∇｛Ｄ(x) ∇ｕ(x,t)｝ ∂ｔ
以下が実行結果です。(Clickで動作／停止 DoubleClickで初期状態に)
バグ？のためにＤの値が変わる境界地で値が不安定になっています。
・・省略・・
上へ
Topに戻る