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

講義日程-2007年度冬学期

数理手法Ⅱ †

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

↑

内容 †

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

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