数理手法Ⅱ のバックアップ(No.15)

講義日程-2007年度冬学期

数理手法Ⅱ †

担当：河村哲也非常勤講師(お茶の水女子大学教授)
1.5単位
- 物工：限定選択
- 数理：限定選択B
- システム：限定選択※
16:30-18:00 工学部六号館 63講義室
教科書
- 数値シミュレーション入門
- 応用編微分方程式

↑

内容 †

偏微分方程式の数値解法とその応用
- 大部分の偏微分方程式は解析的には解けない→数値的に解く。
- 物理法則に現れる偏微分方程式の多くは2階。
  線形なものの方が解きやすい。非線形なものはできる範囲で線形化。

↑

過去問 †

1. 1解波動方程式のLax法での安定条件．
$u_i^j=G^ne^{\sqrt{-1}\xi\Delta x}$ とおき，代入すると，
$|G|=\frac{\Delta t}{\Delta x}$
より&mimetx(\Delta t\le \Delta x);が安定条件．

2. 2次元ラプラス方程式の境界値問題を3×3の格子を用いて解け.
$(x,y)=(ai/3,j/3)$ での解を $u_{i,j}$ と置くと，差分近似して，
$\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{(\Delta x)^2}+\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{(\Delta y)^2}$
また，境界条件より， $u_{0,j}=u(3,j)=0, u_{i,0}=-3, u_{i,3}=3$ なので， $(i,j)=(1,1),(1,2),(2,1),(2,2)$ を入れて連立方程式を解けばよい．

3. 1次元拡散方程式の初期値・境界値問題を変数分離法で解け．
$u=X(x)T(t)$ とおき，
$\frac{1}{X(x)}\frac{\partial^2 X(t)}{\partial x^2}=\frac{1}{T(t)}\frac{\partial^2 T(t)}{\partial t^2}=C^2$
とかおいてごにょごにょ．

4 波動方程式，ラプラス方程式，拡散方程式の解の性質について知っていることを簡潔に述べよ．
- 波動方程式
  
  $\frac{\partial^2 u}{\partial t^2}=k^2\frac{\partial^2 u}{\partial x^2}$
  $u=f(x-kt)+f(x+kt)$ の形で表せる．
- ラプラス方程式
  
  $\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0$
  平均値の定理，最大最小値の定理が成立．
- 拡散方程式の解の性質
  
  $\frac{\partial u}{\partial t}=a^2\frac{\partial^2 u}{\partial x^2}$
  各点で関数の値が平均化する方向に向かう．

5 講義の感想
指がつりました．

↑

授業ノート †

数理物理学に現れる偏微分方程式(PDE,Partial Differential Equation)
⇔常微分方程式(ODE,Ordinal Differential Equation)
- 二つ以上の独立変数に対する偏微分を含んだ方程式。
- ラグランジュの偏微分方程式
  u = u(x,y) は未知
  P,Q,Rは既知として、
  P(x,y,u)∂u/∂x + Q(x,y,u)∂u/∂y = R(x,y,u)~~について、~~補助方程式、dx/P = dy/Q = du/R という連立一階常微分方程式を解き、その解~~f(x,y,z) = a~~g(x,y,u) = b を使って、~~元の方程式の一般解は、任意関数φを使ってφ(f,g) = 0となる。
- ∂u/∂x + c∂u/∂t = 0 を解く。~~dx = dt/c ∧ du = 0 ⇒ u = a ∧ t-cx=bなので、一般解は~~φ(u,t-cx)=0~~uについて解けば u = ψ(t-cx)~~x軸正の方向へ速度1/cで平行移動している。

一次元波動方程式~~{(∂/∂t)^2-(k∂/∂x)^2}u = 0~~微分演算を因数分解。(∂/∂t-k∂/∂x)(∂/∂t-k∂/∂x)u = 0~~ξ=x-kt, η=x+kt と変数変換すると、~~(∂/∂ξ)(∂/∂η) u = 0⇒φ(ξ) + ψ(η)が一般解。

二次元ラプラス方程式~~{(∂/∂x)^2 + (∂/∂y)^2}u = 0 の解は正則関数か、複素共役が正則な関数。

物理現象からの導出
- 弦の微小振動。弦を微小区間に分割して運動方程式を立て、弦の曲がりが小さいとして一時近似。

二階の準線形偏微分方程式の分類
- 拡散方程式。熱伝導方程式。熱流は等温面に垂直に、温度勾配に比例して流れると見なす。~~ある領域Ω内部の熱量変化は、(∂/∂t)∫[Ω]TdV に比例。~~領域表面からの流入量は、温度勾配の面素方向成分に比例すると考えると、∫[∂Ω]gradT dS に比例。~~この二つが比例すると見なすと、(∂/∂t)∫[Ω]TdV = a^2∫[∂Ω]gradTdS = a^2∫[Ω]ΔTdV~~これが任意の領域で成り立つので、∂T/∂t = a^2 ΔT~~比例定数a^2は、質量密度ρ、比熱c、熱伝導率κを使って a^2=κ/(ρc) と書ける。
- ポアソン方程式~~熱源ありの熱伝導方程式の平衡状態の解。~~Δu = -f~~熱源なしの場合はラプラス方程式Δu = 0
- 波動方程式、{(∂/∂t)^2 - Δ}u = 0
- 二階の準線形偏微分方程式：二階の偏微分方程式で、ニ階偏微分について一次の方程式。
  - ２変数の準線形偏微分方程式、A∂[x]^2 u + B∂[x]∂[y]u + C∂[y]^2 u = Dについて、~~ξ = ξ(x,y)~~η=η(x,y)と変数変換すると、~~a = A(∂[x]ξ)^2 + B(∂[x]ξ)(∂[y]ξ) + C(∂[y]ξ)^2~~c = A(∂[x]η)^2 + B(∂[x]η)(∂[y]η) + C(∂[y]η)^2~~b = 2A(∂[x]ξ)(∂[y]η) + B{(∂[x]ξ)(∂[y]η)+(∂[y]ξ)(∂[x]η)} + 2C(∂[y]ξ)(∂[x]η)~~として、a∂[ξ]^2 u + b∂[ξ]∂[η]u + c∂[η]^2 u = dという形に変換される。~~うまくξとηを選べば、a,bを0にすることができる場合がある。~~B-4AC＞0のとき、双曲型。a=c=0となるように選べる。∂[ξ]∂[η]u = fの形、もしくは(∂[ξ]^2-∂[η]^2)u = fの形にまで簡単にできる。~~B-4AC=0のとき、放物型。∂[ξ]^2 u = fの形にまで簡単にできる。~~B-4AC＜0のとき、楕円型。∂[ξ]∂[ξ^*]u = fもしくは(∂[ξ]^2+∂[η]^2)u = f の形にまで簡単にできる。~~標準化した式は、いずれも右辺には一階微分までしか含まない。
- 例){y^2∂[x]^2 - 4x^2 ∂[y]^2}u = 0~~y∂[x]ξ + 2x∂[y]ξ = 0~~y∂[x]η - 2x∂[y]η = 0と変換すると、~~∂[ξ]∂[η]u = (1/(2(ξ^2-η^2)))(η∂[ξ]-ξ∂[η])u

偏微分方程式の解法、 $\LaTeX$ の練習
- 変数分離法
  線形・同次な偏微分方程式で、境界条件も同次なときにはそのまま変数分離可能。
  独立変数ごとの別々な関数の積の形で解を書いてから、境界条件を満たすような重みで重ね合わせる。
  方程式や境界条件が同次でないときには、非同次の特解と同次の一般解を重ね合わせればいい。
  - 有限長の一次元熱伝導
    $u = u(x,t)$
    $\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}$
    境界値 $u(0,t)=u(1,t)=0$
    初期値 $u(x,0)=f(x)$
    与えられた境界値条件(B.C.,Boundary Condition)と初期条件(I.C.,Initial Condition)のもとで
    偏微分方程式の解を求める問題。初期値・境界値問題。
    解の一意性は他で示すとして、発見的に解を見つける。
    変数分離、 $u(x,t) = T(t)X(x)$ という形の解を仮定すると
    T,Xは、 $\frac{dT}{dt}X = T\frac{d^2 X}{dx^2}$ を満たす。
    変形して、 $\frac{1}{T}\frac{dT}{dt} = \frac{1}{X}\frac{d^2 X}{dx^2}$
    左辺は $x$ に依存せず、右辺は $t$ に依存しないので、両辺は定数でなければならない。
    $\frac{1}{T}\frac{dT}{dt} = \frac{1}{X}\frac{d^2 X}{dx^2} = -C^2$
    二つの別々な常微分方程式に帰着できた。
    まず、 $X$ の方を解くと。
    $X = A \sin(Cx) + B \cos(Cx)$
    境界条件を課せば、 $B = 0$ ~、 $A \sin(C) = 0$
    自明でない解は $\sin(C) = 0$ 、 $C=n\pi$ のときのみ存在し、
    $X = A\sin(n\pi x)$
    次に、 $T$ の方を解くと、 $\frac{dT}{dt} = -(n\pi)^2 T$
    $T = \exp(-n^2 \pi^2 t)$
    二つ合わせて、 $TX = C_n \exp(-n^2 \pi^2 t)\sin(n\pi x)$
    これは境界条件を満たす解の一系列になっている。
    次に、元々の偏微分方程式が線形であることに注意して、初期条件を満たすために
    得られた解の線形結合をとる。
    $u = \sum_{n=1}^{\infty}A_n \exp(-n^2 \pi^2 t)\sin(n\pi x)$
    初期条件を満たすように係数 $A_n$ を決めれば最終的な解が求まる。
    $u(x,0) = \sum_{n=1}^{\infty}A_n \sin(n\pi x) = f(x)$
    $u(x,t) = 2\sum_{n=1}^{\infty}\left(\int_{0}^{1}f(y)\sin(n\pi y)dy\right)\exp(-n^2 \pi^2 t)\sin(n\pi x)$
  - 半無限長の一次元熱伝導
    境界条件を $u(x,0)=0$ 、 $x\to\infty$ で $u$ が有界、にする。
    それ以外は同じ。
    変数分離をして $X$ の方を解くと、 $X = A\sin(Cx) + B\cos(Cx)$
    境界条件から $B=0$ 、 $C$ は任意の実数でよい。
    $T$ の方も解いてから重ね合わせると、 $u = \int_0^{\infty}A(C)\exp(-C^2 t)\sin(Cx)dC$
    境界条件を満たすように $A(C)$ を決めると、
    $u(x,t) = \frac{2}{\pi}\int_0^{\infty}\left(\int_0^{\infty}f(y)\sin(Cy)dy\right)\sin(Cx)dC$
  - 熱源のある有限長一次元熱伝導
    熱源の無いときの偏微分方程式に項を追加して、
    $\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + h(x,t)$ とする。他の条件は同じ。
    特解を見つけるために、定数変化法を使う。
    同次方程式の解、 $\sum_{n=1}^{\infty}A_n \exp(-n^2 \pi^2 t)\sin(n\pi x)$ を少し一般化して、
    $\sum_{n=1}^{\infty}A_{n}(x) \exp(-n^2 \pi^2 t)\sin(n\pi x)$ の形で探してみる。
    偏微分法定式に代入してみると。 $A_{n}(x)$ が満たすべき方程式は、
    $\sum_{n=1}^{\infty}\exp(-n^2 \pi^2 t)\{\frac{d^2 A_{n}(x)}{dx^2}\sin(n\pi x) + 2n\pi \frac{dA_{n}(x)}{dx}\cos(n\pi x)\} + h(x) = 0$
    $\TeX$ 化限界により以下略。この方程式を解いて非同次の特解を求め、同次の一般解と重ね合わせれば
    非同次の一般解が得られる。

フーリエ変換による解法
$\hat{f}(k) = \int_a^b K(x,k)f(x)dx$ という形の関数から関数への線形変換を積分変換と呼ぶ。
うまく変換すると方程式本体や境界条件を簡単にすることができる。
$a=-\infty, b=\infty$
$K(k,x)=\frac{1}{\sqrt{2\pi}}\exp(-ikx)$ と選び、~
定義域を適当に設定したものをフーリエ変換と呼ぶ。
- 定義域が急減少超関数なら値域は緩増加関数。
  定義域が緩増加超関数なら値域は急減少関数、などなど。少々うろ覚え。

フーリエ変換
- $\hat{f}(k) = F[f] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(x)e^{-ikx}dx$
  $f(x) = F^{-1}[\hat{f}] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\hat{f}(k)e^{ikx}dk$
- 種々の公式
  $F\left[\frac{d^n f}{dx^n}\right] = (ik)^{n}\hat{f}(k)$
  $F[f*g] = \sqrt{2\pi}F[f]F[g]$ $F[f\cdot g] = \frac{1}{\sqrt{2\pi}}F[f]*F[g]$
  - 畳み込み
    $f*g = \int_{-\infty}^{\infty}f(\xi)g(x-\xi)d\xi$
    畳み込みは線形で可換、零関数以外との畳み込みは単射。
    $\delta$ 関数が単位元。
    定数関数の $1$ との畳み込みは全範囲での積分を表す。
    $(f*g)' = f'*g = f*g'$

円形膜の振動。太鼓の中心を叩く。
$\frac{\partial^2 u}{\partial t^2} = c^{2}\Delta u$
$u = 0$ on $x^2 + y^2 = a^2$
- 円筒座標系で波動方程式を解く。
  $x = r \cos \theta$
  $y = r \sin \theta$
  $\frac{\partial^2 u}{\partial t^2} = c^{2}\left(\frac{\partial ^2}{\partial r^2} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}\right) u$
  $u = 0$ on $r=a$
  $u(r,\theta) = u(r,\theta + 2\pi)$
- 色々略。添付されているyambiさんのtexノート参照。

無限長の弦の振動
- $\ddot{u} = c^2 \Delta u$
  $u(x,0) = F(x)$
  $\dot{u}(x,0) = G(x)$

一般解は
$u(x,t) = f(x-ct) + g(x+ct)$
$u(x,0) = F(x) = f(x) + g(x)$
$\dot{u}(x,0) = G(x) = -cf'(x) + cg'(x)$
$\int_{a}^{x} \frac{G(\xi)}{c}d\xi - f(a) + g(a) = - f(x) + g(x)$
$2f(x) = F(x) - \int_{a}^{x} \frac{G(\xi)}{c}d\xi + f(a) - g(a)$
$2g(x) = F(x) + \int_{a}^{x} \frac{G(\xi)}{c}d\xi - f(a) + g(a)$
$u(x,t) = \frac{1}{2}\left{ F(x-ct) + F(x+ct) \right} + \frac{1}{2}\int_{x-ct}^{x+ct}G(\xi)d\xi$
初期値が影響を与える領域は時間と共に広がってゆくが、有限の広さを持つ。

円板の平衡温度分布
- 半径 $a$ の円板の円周上での温度分布がわかっているとき、
  平衡状態での内部の温度分布を求める。
  $\Delta u = 0$
  境界条件として、 $u = g \mathrm{on} \Gamma$
- 軸対称な問題なので、極座標で考える。極座標に移れば、
  $\Delta u = \left{\frac{\partial^2}{\partial r^2} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}\right}u=0$
  $u(a,\theta) = g(\theta)$
  変数分離して、 $u = R(r)\Theta(\theta)$
  $\Theta \frac{d^2 R}{dr^2} + \frac{\Theta}{r}\frac{dR}{dr} + \frac{R}{r^2}\frac{d^2 \Theta}{d\theta^2}=0$
  $\frac{r^2}{R}\frac{d^2 R}{dr^2} + \frac{r}{R}\frac{dR}{dr} = -\frac{1}{\Theta}\frac{d^2 \Theta}{d\theta^2} = m^2$
  $\Theta$ は周期性を持つ必要があるので、 $m$ は非負整数。
  $\Theta = A \sin (m\theta) + B \cos(m\theta)$
  $\left(r^2\frac{d^2}{dr^2} + r\frac{d}{dr} - m^2 \right)R = 0$
  $\left{\left(r\frac{d}{dr}\right)^2 - m^2 \right}R = 0$
  $\left(r\frac{d}{dr}- m\right)\left(r\frac{d}{dr} + m\right)R = 0$
  Rが原点で有界であることを課せば、
  $R = Cr^m$
  得られた解の群を重ね合わせて一般解を作れば、
  $u(r,\theta) = B_0 + \sum_{m=1}^{\infty}r^{m}(A_m \sin(m\theta) + B_m \cos(m\theta))$
  境界条件を満たすように係数を決めると、
  $u(a,\theta) = B_0 + \sum_{m=1}^{\infty}a^{m}(A_m \sin(m\theta) + B_m \cos(m\theta)) = g(\theta)$

$u(a,\theta) = \frac{1}{2\pi}\int_{0}^{2\pi}g(\xi)d\xi \frac{a^2-r^2}{a^2 -2ar\cos(\theta - \xi) + r^2}$

ラプラス方程式 $\Delta u = 0$ の解を、調和関数と呼ぶ。
- 「調和関数の平均値の定理」が成り立つ。
  $u(p) = \frac{1}{2\pi a}\oint_{C}u |dr|$ 、 $C$ はpを中心とする半径aの円。
  三次元では、 $u(p) = \frac{1}{4\pi a^2}\int_{\partial \Omega}u |dS|$ 、 $\Omega$ はpを中心とする半径aの球。
- 最大・最小の定理
  定数でない調和関数は領域内部では極大・極小をとらない。
- ラプラス方程式、ディリクレ境界問題の解の一意性。境界での関数値を指定する。
  ２つの解の存在を仮定すると、その差もラプラス方程式を満たし、差は境界で0となる。
  ラプラス方程式の解は境界でのみ最大値・最小値を取れるので、最大値も最小値も0。
  よって２解の差は常に0となり、一致する。

１階常微分方程式の初期値問題の数値解法
$y'(x)=f(x,y(x))\ \wedge\ y(0)=y_0$
- オイラー法
  非常に単純。精度はとても悪い。
  $y' = \frac{dy}{dx} = \lim_{\Delta x\to 0}\frac{y(x+\Delta x)- y(x)}{\Delta x} \sim \frac{y(x+h)- y(x)}{h}$
  $\frac{y_{n+1} - y_{n}}{h} = y'_n = f(x_{n},y_{n})$
  $y_{n+1} = y_{n} + h f(x_{n},y_{n})$
  - 精度の向上
    $y(x+\Delta x) = y(x) + y'(x)\Delta x + o(\Delta x)$
    オイラー法では、解のテイラー展開のうち、
    $\Delta x$ の一次までとって他を無視している。
    もっと高次の項までとれば解の精度がよくなる。
    なるべく手間をかけずに高階微分の項を取り入れる方法を探す。
- 4次のRunge-Kutta法
  xをずらした点の値を使って、高階微分の項を取り入れる。
  - $k_1 = f(x_n,y_n)$
    $k_2 = f(x_n + \frac{h}{2},y_n + k_{1}\frac{h}{2})$
    $k_3 = f(x_n + \frac{h}{2},y_n + k_{2}\frac{h}{2})$
    $k_4 = f(x_n + h,y_n + k_{3}h)$
    $y_{n+1} = y_n + (k_1 + 2k_2 + 2k_3 + k_4)\frac{h}{6}$
数値積分の利用
- $y'(x)=f(x,y(x))$ の両辺を $[x_n,x_{n+1}]$ の区間で積分すると、
  $y_{n+1} - y_{n} = \int_{x_n}^{x_n+1}f(x,y(x))dx$
  右辺には未知関数の積分が入っているのですぐには値を求められない。
  右辺を長方形近似したものがオイラー法。
  $\int_{x_n}^{x_n+1}f(x,y(x))dx \sim f(x_{n},y_{n})\Delta x$
  まずオイラー法で予測子として仮の $y_{n+1}$ の値を求め、
  台形公式で近似するのがホイン法。
  現在の値と１ステップ前の値から直線近似で次の値を求めるのがアダムス・バッシュフォース法

連立微分方程式と高階微分方程式の初期値問題
- 連立微分方程式
  既知の初期値から始まり、微分が既知の関数として与えられている。
  解法は１変数のときと同じ
  $x(t) \in \mathbb{R}^{n}$
  $\dot{x} = f(t,x)$
  $x(0) = x_0$
  - オイラー法
    $x_{n+1} = x_{n} + f(t_{n},x_{n})\Delta t$
  - ホイン法
    $k_1 = f(t_{n},x_{n})$
    $k_2 = f(t_{n}+\Delta t,x_{n}+k_1 \Delta t)$
    $x_{n+1} = x_{n}+(k_1+k_2)\frac{\Delta t}{2}$
  - ４次Runge-Kutta法
    $k_1 = f(t_{n},x_{n})$
    $k_2 = f(t_{n}+\frac{\Delta t}{2},x_{n}+k_1 \frac{\Delta t}{2})$
    $k_3 = f(t_{n}+\frac{\Delta t}{2},x_{n}+k_2 \frac{\Delta t}{2})$
    $k_4 = f(t_{n}+\Delta t,x_{n}+k_3 \Delta t)$
    $x_{n+1} = x_{n}+(k_1+2k_2+2k_3+k_4)\frac{\Delta t}{6}$
- 高階微分方程式
  高階の微分方程式は変数を増やすことで連立微分方程式に帰着可能 $\ddot{x} = f(t,x,\dot{x})$ は $y = \dot{x}$ と置けば
  $\dot{y} = f(t,x,y)$
  $\dot{x} = y$

数値微分
離散的なサンプル値だけが既知の関数の微係数の近似値を計算する。
- 求めたい点の周りの関数値の線形結合で微係数を近似計算する。
  周り３点を使うと、
  $ax(t-k) + bx(t) + cx(t+h) \sim (a+b+c)x + (ch-ak)\dot{x} + (ak^2+ch^2)\ddot{x}$
  $a+b+c = 0$
  $ch-ak = 1$ と選べば近似値が得られ、
  さらに $ak^2 + ch^2$ とすれば精度が上がる。
  $a=0\ hb=-1\ hc=1$ と選んだものを前進差分。
  $ka=1\ kb=-1\ c=0$ と選んだものを後退差分。
  $k=h$ のときに $ha=2^{-1}\ hb=-1\ hc=2^{-1}$ と選んだものを中心差分と呼ぶ。

境界値問題
$\ddot{y}+y+x=0\ (0<x<1)$
$y(0)=y(1)=0$ と、
領域の境界での制約が与えられた問題を境界値問題と呼ぶ。
今回の微分方程式の一般解は $y=A\cos(x)+B\sin(x)-x$
境界条件を課すと $y=\frac{\sin(x)}{\sin{1}}-x$
- 差分法による近似解法
  空間を離散化し、格子点上での近似解を求める。
  各格子点に番号を振り、微分を差分で置き換えて、
  微分方程式を差分方程式に近似する。
  $\ddot{y}_{n} = \frac{y_{n+1}+y_{n-1}-2y_{n}}{h^2}$
  $y_{n+1}+y_{n-1}-2y_{n} + h^{2}(y_{n}+x_{n}) = 0$
  $y_{n+1}-(2-h^{2})y_{n}+y_{n-1}=-h^{3}n\ (1<n<N)$
  $y_{1}=y_{N}=0$
  線形微分方程式を連立一次方程式に変換できた。
  この方程式は対角要素とその上下一つだけが非ゼロの帯状行列。
  対角より下の要素を消すように掃き出し、後退代入で下から値を定めていけば
  $O(n)$ で解ける。

偏微分方程式の差分解法
- １次元拡散方程式の陽解法
  $\frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}\ 0<x<1,\ t>0$
  境界条件 $u(0,t)=u(1,t)=0$ と
  初期条件 $u(x,0)=f(x)$ の
  両方が与えられている問題は初期値境界値問題。
  問題を解く領域を固定して、時間・空間を等間隔に離散化し、
  偏微分を差分に置き換える。
  $u_{j}^{n} = u(x_j,t_n)$
  $\frac{\partial u}{\partial t} = \frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}$
  $\frac{\partial^2 u}{\partial x^2} = \frac{u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{(\Delta x)^2}$
  $\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t} = \frac{u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{(\Delta x)^2}$
  $u_{j}^{n+1} = u_{j}^{n} + \frac{\Delta t}{(\Delta x)^2}(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n})$
  これはnについての漸化式。初期値と境界値が決まれば次々に決めて行くことができる。

差分近似で刻み幅を大きくとりすぎてしまうと解が振動してしまい、
誤差がどんどん大きくなってしまう。
この熱伝導方程式の場合、その境目は0.5。

$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}$ に
$u(t,x) = g(t)e^{-i\xi x}$ の形の解を仮定すると、
$g=e^{-\xi^2 t}$ が得られる。
差分方程式の方にも特解として $u_{j}^{n}=G^{n}e^{\xi j\Delta x}$ の形の解を仮定する。
$\frac{u_{j}^{n+1}}{u_{j}^{n}} = G$ は増幅率と呼ばれる。
差分方程式が有界な解を持つためには、 $|G|\leq 1$ が必要。
$u_{j}^{n+1} = u_{j}^{n} + \frac{\Delta t}{(\Delta x)^2}(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n})$
$G = 1 + \frac{\Delta t}{(\Delta x)^2}(Ge^{\xi i\Delta x}-2+e^{-\xi i\Delta x})$
$G = 1-4\frac{\Delta t}{\Delta x^2}(\sin\frac{\theta}{2})^2$
$\frac{\Delta t}{\Delta x^2}\leq \frac{1}{2(\sin\frac{\theta}{2})^2}\leq \frac{1}{2}$ ならば解は収束する。

陰解法
時間微分を差分に直す際、前進差分ではなく後退差分を取ってみる。
$\frac{\partial u}{\partial t} = \frac{u_{j}^{n}-u_{j}^{n-1}}{\Delta t}$
$\frac{u_{j}^{n}-u_{j}^{n-1}}{\Delta t} = \frac{u_{j-1}^n-2u_{j}^{n}+u_{j+1}^{n}}{(\Delta x)^2}$
これは、漸化式型ではなく、連立一次方程式になっている。
陰解法では $G=\frac{1}{1+4\frac{\Delta t}{(\Delta x)^2}(\sin\frac{\theta}{2})^2}$ 、刻み幅の取り方に依らずに安定。

陰解法と陽解法を組み合わせて、 $\frac{u_{j}^{n}-u_{j}^{n-1}}{\Delta t} = (1-\alpha)\frac{u_{j-1}^n-2u_{j}^{n}+u_{j+1}^{n}}{(\Delta x)^2} + \alpha \frac{u_{j-1}^{n+1}-2u_{j}^{n+1}+u_{j+1}^{n+1}}{(\Delta x)^2}$
と取ると、 $\alpha\geq\frac{1}{2}$ のとき $|G|\leq 1$ になる。
$\alpha=\frac{1}{2}$ のとき時間精度が一つ向上する。

ADI法(Alternating Direction Implicit Method)
二次元の熱伝導方程式の陰解法では、巨大な疎行列がかかった連立一次方程式が現れる。
計算時間を削減するために、特定の方向だけを陰的に、他の方向を陽的にし、
各ステップで陰的にする方向を交互に切り替えて計算すると
計算量を削減しつつ安定性を確保できる。
２次元ならうまくいくが、３次元ではこの方法では安定性は確保できない。

移流方程式の差分解法
- 1次元の波動方程式
  $\ddot{u} = c^{2}\Delta u$ を変形すると、
  $\left(\frac{\partial}{\partial t}-c\frac{\partial}{\partial x}\right)\left(\frac{\partial}{\partial t}+c\frac{\partial}{\partial x}\right)u = 0$
  この方程式は特解として、
  $\left(\frac{\partial}{\partial t}+c\frac{\partial}{\partial x}\right)u = 0$
  という一次元移流方程式を含む。
  移流方程式は厳密に解けて、
  初期条件 $u(x,0)=f(x)$ に対して
  解は $u(x,t)=f(x-ct)$
- 差分解法
  微分を差分で置き換えると $\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+c\frac{u_{j+1}^{n}-u_{j-1}^{n}}{2\Delta x}=0$
  $u_{j}^{n+1} = u_{j}^{n} - \frac{\mu}{2}(u_{j+1}^{n}-u_{j-1}^{n})$ 、 $\mu = \frac{c\Delta t}{\Delta x}$
  さて、元の方程式の厳密な特解を使って安定性を調べる、
  $u_{j}^{n}=G^{n}e^{i\xi j\Delta x}$ と置くと、
  $G=1-i\mu \sin \theta$
  $|G|^2 = 1+\mu^2 \sin^2 \theta > 1$
  この解法では解は常に不安定になってしまう。
- 差分解法その2
  今度は中心差分ではなく後退差分で近似する。
  $\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+c\frac{u_{j}^{n}-u_{j-1}^{n}}{\Delta x}=0$
  この差分は1次精度上流差分と呼ぶ。
  $u_{j}^{n+1} = (1-\mu)u_{j}^{n} - \mu u_{j-1}^{n}$
  再び厳密解を使って安定性を判別する。
  $G=1-\mu e^{-i\theta}$ 、 $\theta = \xi \Delta x$
  これは複素平面上で $(1-\mu,0)$ を中心とした半径 $\mu$ の円を描く。
  $|G|\leq 1$ であり、近似解が安定になる条件は、
  $\frac{c\Delta t}{\Delta x} = \mu\leq 1$
  特性曲線が影響領域に含まれる必要がある。(CFL条件)
- 精度の向上(Lax,Wedroff法)
  時間微分を差分に置き換える際、2次の項まで残し、
  空間微分を中心差分に置き換えると、
  $u_{j}^{n+1}=u_{j}^{n}-\frac{\mu}{2}(u_{j+1}^{n}-u_{j-1}^{n})+\frac{\mu^2}{2}(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n})$
  時間空間ともに2次精度の公式が得られる。
  安定条件は $\mu \leq 1$
- 安定条件の発見的な見つけ方
  1次精度の上流差分近似をしたものをテイラー展開し、2次の項まで残すと、
  $\frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}+ c\frac{u(x,t)-u(x-\Delta x,t)}{\Delta x} = 0$
  $\frac{\partial u}{\partial t}+ c\frac{\partial u}{\partial x} = -\frac{\Delta t}{2}\frac{\partial^2 u}{\partial t^2} + \frac{c\Delta x}{2}\frac{\partial ^2 u}{\partial x^2} = \frac{c\Delta x}{2}(1-\mu)\frac{\partial^2 u}{\partial x^2}$
  となり、近似によって移流法定式が移流拡散方程式になっている。
  もし拡散定数 $1-\mu$ 負なら、波は拡散せず
  徐々に中心に集まって尖ってゆきいずれ発散してしまう。
  厳密解と増幅率を使った安定判別は線形の偏微分方程式にしか使えないが、
  差分をとったものを再びテイラー展開する方法は少々非線形でも使える場合がある。
- クーラン数の意味
  $\mu = \frac{c\Delta t}{\Delta x} = \frac{c}{\frac{\Delta x}{\Delta t}}$
  分母は数値的な情報の伝わる速さ。分子は物理的な情報の伝わる速さ。
  数値的な情報の伝播は、物理的な情報の伝播に追いつかれてはいけない。

差分スキーム
差分方程式で、 $\Delta x$ と $\Delta t$ の間に関数関係を仮定したもの。
差分近似をしたときに、差分演算子 $S$ を使って、
$u(x,t+\Delta t)-S(h(\Delta t),\Delta t)u(x,t) = o(\Delta t)$
にの形になるなら、差分スキームは偏微分方程式に適合している、という。
$|S^{n}|<c$ ならばその差分スキームは安定であるという。
- Laxの同等定理。
  線形部分方程式を近似する差分スキームが、安定で偏微分方程式に適合しているならば、
  $\Delta t\to 0$ のとき、差分スキームの解は微分方程式の解に収束する。

ポアソン方程式の数値解法
$\Delta u = -\rho$

ラプラシアンの意味
ラプラシアンを離散化すると、ラプラシアンの作用は
各点について、周囲の点での値の平均からのずれを表している。
- cf)調和関数の平均値の定理,最大最小の定理
- 拡散方程式は、各点での関数の値が平均化する方向へ変化して行くことを表現している。

優調和関数，劣調和関数
- $u_p \geq \frac{1}{2\pi a}\oint_c u ds$
- $u_p \leq \frac{1}{2\pi a}\oint_c u ds$

差分法で得られる連立1次方程式系の特徴
- 大次元で非常に規則正しい疎行列
  疎行列の逆行列は疎行列になるとは限らない。⇒反復解法
- $Ax = b$ を $x = Mx+c$ と変形し、
  $x_{n+1} = Mx_{n}+c$ という漸化式を反復する。
  $x_{n}$ が収束すれば、収束先の微分法定式の解。
  $M$ のスペクトル半径が1より小さければ収束する。
- ヤコビの反復法
  $Ax = b$
  $A = L+U+D$ と分解。Lは下三角行列、Uは上三角行列、Dは対角行列。
  Dの対角要素には0が現れないように適当に入れ替えておく。
  $Ax = (L+U+D)x = b$
  $x = -D^{-1}(L+U)x + D^{-1}b$
- ガウス・ザイデル法
  $(L+D+U)x=b$
  $x_{n+1} = -(L+D)^{-1}Ux_{n} + (L+D)^{-1}b$
  $x_{n+1} = -D^{-1}Lx_{n+1} - D^{-1}Ux_{n} + D^{-1}b$
  $D^{-1}L$ は対角要素の無い下三角行列。上から順に、陽に計算できる。
  スペクトル半径はヤコビ法の半分ほどになり、半分ほどの反復で求められる。
- SOR法(Successive Over Relaxation Method)
  $x'_{n} = -D^{-1}Lx_{n+1}-D^{-1}Ux_{n}+D^{-1}b$
  $x_{n+1} = (1-\alpha)x_{n}+\alpha x'_{n}$
  $\alpha$ は加速係数。収束するためには $0<\alpha<2$ が必要。
  $\alpha$ を最適値に選ぶと、運が良ければヤコビ法の10〜20倍速くなる。
  $\alpha$ の値は $b$ には依らない。
  理論的に決める方法は無く、試行により決める。

差分法の応用
- 複雑な領域の取扱い
  物理的な意味のある領域形状を、計算しやすい領域に座標変換。
  物理面から計算面への変換。
- 物理面の格子点での差分を、格子点の座標の差分で割れば物理量の差分近似になっている。
  複雑な形状の物理面をパラメタ付けし、パラメタが整数の点での物理量の値と座標値の差分を使う。

数値流体力学(CFD,Computational Fluid Dynamics)
基礎方程式は古くから知られていたが、非線形方程式であるため解析的には解きにくい。
- 非圧縮性流体の基礎方程式
  流速が音速に比べて0.3以下くらいなら非圧縮と近似してよい。
  $\mathrm{div}v = 0$
  $\frac{\partial v}{\partial t} + (v\cdot \nabla)v = - \mathrm{grad}\phi + \nu \Delta v$
  $v$ ：流速
  $\phi=\frac{P}{\rho}$ ：(圧力)/(密度)
  運動量の時間変化が、流れによる運動量の輸送と、
  圧力差による力、粘性力によって決まることを記述している。
- 離散化
  もし圧力項が無ければ、
  $\frac{\partial v}{\partial t} = \frac{v^{n+1}-v^{n}}{\Delta t}$ として、
  $v^{n+1} = v^{n} + \Delta t f(v^{n})$ の形にできる。
  まずは圧力項を無視して1ステップ進め、圧力項の分を補正する。
  $\frac{v^*-v^n}{\Delta t} + (v^{n}\cdot \nabla)v^n = \nu \Delta v^n$
  $\frac{v^{n+1}-v^n}{\Delta t} + (v^{n}\cdot \nabla)v^n = -\mathrm{grad}\phi^n + \nu \Delta v^n$
  $\frac{v^*-v^{n+1}}{\Delta t} = \mathrm{grad}\phi^n$
  $\frac{\mathrm{div} v^*- \mathrm{div} v^{n+1}}{\Delta t} = \Delta \phi^n$
  $\Delta \phi^n = \frac{\mathrm{div} v^*}{\Delta t}$ ポアソン方程式として解くと $\phi^n$ を得られる。
  $v^{n+1} = v^* - \Delta t \mathrm{grad}\phi^n$
  この方法はフラクショナルステップ法と呼ばれる。
  連続の式が近似的なものになることと、ポアソン方程式の収束が遅いことが難点。
  静止壁についての境界条件は、壁と流体の相対速度を0に設定する、圧力に関しては運動方程式から決める。