３−４．偏微分方程式のまとめと演習

Elliptic （楕円型） Equations

Laplace Equation:  ∂²u/∂x² + ∂²u/∂y² = 0

Poisson Equation:  ∂²u/∂x² + ∂²u/∂y²  + f (x, y) = 0

• 線形の多元連立方程式に還元されるので、４章の行列演算で取り扱う。

　

Hyperbolic （双曲型） Equations

Wave Equation:  ∂²y/∂t² - c² ∂²y/∂x² = 0

• Simple method
  
  u_jⁿ⁺¹ = u_jⁿ - cdt/2dx (u_j+1ⁿ - u_j-1ⁿ)
  y_jⁿ⁺¹ = y_jⁿ + cdt/2dx (y_j+1ⁿ - y_j-1ⁿ + 2u_jⁿ⁺¹)

• the Lax method
  
  u_jⁿ⁺¹ = (u_j+1ⁿ + u_j-1ⁿ )/2 - cdt/2dx (u_j+1ⁿ - u_j-1ⁿ)
  y_jⁿ⁺¹ = (y_j+1ⁿ + y_j-1ⁿ )/2  + cdt/2dx (y_j+1ⁿ - y_j-1ⁿ + 2u_jⁿ⁺¹)

　

Parabolic （放物型） Equations

Diffusion Equation:  ∂u/∂t - k ∂²u/∂x² = 0

• Simple method
  
  u_jⁿ⁺¹ = u_jⁿ - kdt/dx² (u_j+1ⁿ - 2u_jⁿ + u_j-1ⁿ)

• the Crank-Nicholson method
  
  u_jⁿ⁺¹ = u_jⁿ - kdt/2dx² (u_j+1ⁿ⁺¹ - 2u_jⁿ⁺¹ + u_j-1ⁿ⁺¹ + u_j+1ⁿ - 2u_jⁿ + u_j-1ⁿ)
  
  （これは陰解法。消去法で連立方程式を解くか、反復法で解を求める。）

　

        《演習》

1. 初期温度分布  u(x,0) = 2x  [0 < x < 0.5]、2(1-x)  [0.5 < x < 1]、
   境界条件  u(0,t) = u(1,t) = 0 [for any t]、として時刻t=0.1における温度分布を計算する。
   dx=0.1を固定し、dtを変化させてSimple methodの安定性を議論せよ。a=kdt/dx²がどのあたりで不安定化するか？
   （Diffusion Equationで熱拡散係数k=1とする。）
   
   サンプルプログラムを参考にしてもよい。[fortran, C]
   f77 -O -o simple simple.f
   ./simple
   
   結果はgnuplotを使ってグラフ化せよ。
   gnuplot
   .....
   gnuplot> plot 'fort.10' w lines, 'fort.11' w points     ←解析解との比較
   
   解析解：　u(x,t) = 8/p²��_i=n^∞(1/n²)sin(np/2)sin(npx)exp(-n²p²t)  (data at t=0.1[fort.10])
   
   動画[a=0.1, 0.5, 0.6]
   　
2. 同じ問題でGaussの消去法を使ったCrank-Nicholson methodを試せ。[fortran, C]
   　
3. さらにGauss-Seidel反復法を使ったCrank-Nicholson methodを試せ。[fortran, C]
   　
4. （アドバンストな自習問題）臨界点近傍にある系を現象論的に記述する方程式としてTDGL(Time Dependent Ginzburg-Landau)式がある。磁性など秩序変数(j)が保存量でない場合には、系の時間発展は以下の非保存型TDGL式で記述できる。この放物型の偏微分方程式をSimple methodで数値シミュレーションするプログラムを作り、bの正負で系のダイナミクスがどう変化するかを議論せよ。
   
   GL Free Energy:  F=∫dxdy ( j⁴/4 + bj²/2 + c(∇j)²)
   TDGL方程式（非保存）:  ∂j/∂t = -L dF/dj = -L(j³ + bj - c ∇²j) = 0
   
   簡単のため、b=±0.5、c = 0.1、L = 1、dx = 1、dt = 0.1とし、空間の次元は２次元とする。システムサイズは128x128で周期的境界条件、初期条件はj=±0.025の範囲でランダムに与えよ。ダイナミクスの議論には、例えば界面積（２次元では界面長）が時間に対してどのように変化するのかを検証せよ。
   
   Simple methodで書いたサンプルプログラムはここ [fortran, C]
   f77 -O -o tdgl_simple tdgl_simple.f
   ./tdgl_simple
   
   gnuplot
   .....
   gnuplot> set contour
   gnuplot> set data style lines
   gnuplot> splot 'fort.10'　　　　　　　←秩序変数(j)の空間分布
   gnuplot> set logscale xy
   gnuplot> plot 'fort.11' u 1:2　　　　 ←時間-界面エネルギーの両対数プロット
   gnuplot> plot 'fort.11' u 1:4　　　 　←時間-界面長の両対数プロット
   
   動画[mpeg]
   　
5. （アドバンストな自習問題）気液相転移など秩序変数(j)が保存量である場合には、系の時間発展は以下の保存型TDGL式で記述できる。Simple methodで数値シミュレーションする２次元 のプログラムを作り、系のダイナミクスを非保存の場合と比較せよ。
   
   TDGL方程式（保存）:  ∂j/∂t = L ∇²dF/dj = L∇²(j³ + bj - c ∇²j) = 0
   
   Simple methodで書いたサンプルプログラムはここ [fortran]
   
   動画[mpeg]
   　
6. （さらにアドバンストな自習問題）上の４、５についてCrank-Nicholson法への改良を試みよ。