拡散型偏微分方程式
ある関数 $u(x,t)$ の、空間（$0 \leq x \leq 1$）・時間（$t \geq 0$）発展が、次のような拡散型偏微分方程式にしたがっているとします。
$$\frac{\partial u}{\partial t} = \frac{\partial ^2 u}{\partial x^2} \tag{1}$$
ここで、初期条件（$t=0$）および境界条件（$x=0,1$）として
$$u(x,0) = f(x),\ \ \ \ \ u(0,t) = f(0),\ \ \ \ \ u(1,t) = f(1) \tag{2}$$
をとります。解 $u(x,t)$ のふるまいを数値的に調べてみましょう。
空間メッシュを $I$ 等分して、それぞれのメッシュに番号をつけます（$x_i;\ x_0 \equiv 0,\ x_I \equiv 1$）。すなわち、隣り合うメッシュの間隔を $\Delta x = 1/I$ とすると
$$x_i \equiv i \Delta x\ \ \ \ \ (i=0, \cdots, I)$$
となります。同様に時間軸も差分化して番号をつけます（$t_j;\ t_0 = 0$）。時間間隔を一定にとり、その値を $\Delta t$ と書くと
$$t_j \equiv j \Delta t\ \ \ \ \ (j=0, 1, \cdots)$$
となります。以下、時刻 $t_j$ における $x_i$ での $u$ の値を $u_{i,j}$ と表すことにします。
さて、時刻 $t_j$ における $u$ の値 $u_{i,j}$ が全てわかったものとして、次の時間ステップにおける $u$ の値 $u_{i, j+1}$ を求める場合について考えましょう。このとき、基本微分方程式 (1) の差分の仕方には、**陽解法（explicit method）および陰解法（implicit method）**の2つの方法があります。以下、それぞれの方法について簡単に説明します。
陽解法（explicit method）
式 (1) の微分方程式の右辺を、古い時刻（$t_j$）で評価すると
$$\frac{u_{i j+1} - u_{i,j}}{\Delta t} = \frac{u_{i+1, j} - 2u_{i,j} + u_{i-1, j}}{(\Delta x)^2},\ \ \ \ \ (i = 1, \cdots, I-1) \tag{3}$$
となります。また境界条件 (2) より
$$u_{0,j+1} = u_{0,j},\ \ \ \ \ u_{I, j+1} = u_{I, j}$$
となります。さて式 (3) の中で未知数は、各 $i$ について $u_{i, j+1}$ だけであることに注意してください。すなわち式 (3) を変形して
$$u_{i,j+1} = u_{i,j} + \frac{\Delta t}{(\Delta x)^2} (u_{i+1,j} - 2u_{i,j} + u_{i-1,j}),\ \ \ \ \ (i = 1, \cdots, I-1)$$
から $u_{i, j+1}$ を求めることができます（右辺は全て既知数）。これが陽解法（explicit method）です。プログラミングは実に簡単で、単に上式を書き下すだけです。
ところが、この方法には重大な欠点があります。あまり大きな時間ステップを取ると、数値不安定を起こして解けなくなってしまうのです。この数値不安定を避けるためには、時間ステップ $\Delta t$ を以下のとおり小さく取らないといけません。
$$\Delta t \leq \frac{(\Delta x)^2}{2}$$
陰解法（implicit method）
式 (1) の微分方程式の右辺を、新しい時刻（$t_{j+1}$）で評価すると
$$\frac{u_{i j+1} - u_{i,j}}{\Delta t} = \frac{u_{i+1, j+1} - 2u_{i,j+1} + u_{i-1, j+1}}{(\Delta x)^2},\ \ \ \ \ (i = 1, \cdots, I-1) \tag{4}$$
となります。境界条件は陽解法と同じです。さて式 (4) において未知数は、各 $i$ につき $u_{i-1, j+1}$, $u_{i, j+1}$, $u_{i+1, j+1}$ と3つあります。すなわち式 (4) は空間的にグローバルに結合していて、そのままでは解けません。そこで、プログラミングに特別な工夫が必要となります。
今 $b \equiv (\Delta x)^2 / \Delta t,\ a \equiv 2+b$ とすると、
$$u_{i-1, j+1} -a u_{i,j+1} + u_{i+1, j+1} = -b u_{i,j}$$
となり、式 (4) は次のような行列方程式に帰着されることになります。
\left(
\begin{array}{ccccccc}
b & 0 & 0 & \ldots & 0 & 0 & 0 \\
-1 & a & -1 & \ldots & 0 & 0 & 0 \\
0 & -1 & a & \ldots & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & \ldots & -1 & a & -1 \\
0 & 0 & 0 & \ldots & 0 & 0 & b
\end{array}
\right)
\left(
\begin{array}{c}
u_{0, j+1} \\
u_{1, j+1} \\
u_{2, j+1} \\
\vdots \\
u_{I-1, j+1} \\
u_{I, j+1}
\end{array}
\right) = b
\left(
\begin{array}{c}
u_{0, j} \\
u_{1, j} \\
u_{2, j} \\
\vdots \\
u_{I-1, j} \\
u_{I, j}
\end{array}
\right)
これを $I$ 行 $I$ 列の行列 $A$ と、$1$ 行 $I$ 列の行列 $U_j$, $U_{j+1}$ を用いてシンボリックに
$$A \cdot U_{j+1} = bU_{j}$$
と書くと、解は
$$U_{j+1} = b A^{-1} \cdot U_j$$
となります。このように、陰解法（implicit method）を用いるとプログラミングは格段に複雑になりますが、数値不安定の条件を考慮しなくても、方程式を安定に解くことができるようになります。
課題
以下の初期条件のもとで、拡散方程式 (1) を陽解法および陰解法で解きなさい。
$$f(x) = \exp \left[-\frac{(x-0.5)^2}{0.25}\right]$$
また以下の初期条件のもとで、拡散方程式 (1) を陽解法および陰解法で解きなさい。
f(x) = \left\{
\begin{array}{l}
1.0\ \ \ \ \ \mathrm{for}\ \ 0 \leq x \leq 0.5 \\
0.0\ \ \ \ \ \mathrm{for}\ \ 0.5 < x \leq 1.0
\end{array}
\right.
それぞれ、陽解法の場合に時間ステップを大きく取ると数値不安定を起こすことを確認しなさい。