Also known as the Forward Euler Finite Difference Method.
Consider the heat equation
$$\dfrac{\partial{u}}{\partial{t}} = \dfrac{\partial^2{u}}{\partial{x^2}}$$
Now we must establish some boundary conditions in order to obtain a solution just as we did with simple ODEs. Let our sample problem in this case be a bi-metallic bar (or an insulated rod) which is maintained at zero degrees on its outer circumference. Cross-sectionally,
$$0 \leq x\leq 1, \quad t \geq 0$$
is the range of reference points from left to right across the section. Think: $x=0.5$ is the hot core (perhaps from a previous thermal or nuclear process) whereas $x=0$ is the exterior cold periphery on the left (top inner boundary in figure below) and $x=1$ the cold periphery on the right (bottom inner boundary in figure below) both of which are maintained at zero degrees C.
$$
\begin{array}{|rrrrrrr|}
\hline
\qquad ^{0^o C} & & & & & & ^{0^o C} \\
\hline \\
——&——&——&—\updownarrow u(x,t) = u(0.5,t) = f(x) \updownarrow —&——&——&——\\
& & & & & & \\
\hline
& \qquad _{0^o C} & \qquad & \qquad & \qquad & \qquad_{0^o C} \\
\hline
\end{array}
$$
$$u(x,t) \equiv \text{temperature at (x, t)}$$
such that at $$t=0, \quad u(x,0) = f(x)$$ and $$ u(0,t) = u(1,t) = 0$$
Continuing from the discrete approximations obtained for the partial derivatives of the heat equation we can create a simple method for future value of $u(x,t)$. Starting with the heat equation
$$ \dfrac{u(x,t + \Delta t ) \; – \; u(x,t)}{\Delta t} \simeq \dfrac{\partial{u}}{\partial{t}} = \dfrac{\partial^2{u}}{\partial{x^2}} \simeq \dfrac{u(x + \Delta x, t )\; – \; 2u(x,t) \; + \; u(x – \Delta x,t)}{(\Delta x)^2}$$
rearranging,
$$ u(x,t + \Delta t ) \simeq u(x,t)+ \dfrac{\Delta t}{(\Delta x)^2} \Biggl( u(x + \Delta x, t )\; – \; 2u(x,t) \; + \; u(x – \Delta x,t) \Biggr) \tag{1}$$
For a considerably neater notation, write
$$ h \equiv \Delta x \; , \qquad k \equiv \Delta t$$
and let
$$h = \frac{1}{M} \;, \qquad M\in \mathbb{Z_+}$$
be the space discretisation time-step as illustrated below. $M$ is a positive integer, hence $ h\equiv \Delta x $ is small!
$$
\begin{array}{c|rrrrrrr}
t & \qquad & \qquad & \qquad & \qquad & \qquad & \qquad \\
\downarrow \qquad & * & * & * & * & * & * \\ \\
k \qquad & & & & & & \\ (\Delta t)\quad & & & & & & &\quad {} \\
\hline
\uparrow \qquad \mathbf{0} &\leftarrow h\rightarrow * &\leftarrow h\rightarrow * &\leftarrow h\rightarrow * &\leftarrow h\rightarrow * & \leftarrow h\rightarrow * & \leftarrow h\rightarrow \mathbf{1} & \quad{x}\\
& \qquad & \qquad & \qquad & \qquad & \qquad & _{(M=6)}
\end{array}
$$
Now we need to define some spatial and time boundary conditions. At $t=0$, let the peripheral temperature of the insulated bar or rod be $0^oC$.
$$U_j^o = u(jh, 0) \cdots \text{known initial temperature at x=0 and x = 1 of 0 deg C}$$
Illustratively, at time $t = 0$ we have a known discrete temperature distribution $u(x + \Delta x, t=0 )$ across the bar’s cross-section
$$u(0,0)_{\text{bound } x=0}\;\underbrace{ \mathbf{|| \blacktriangleright}\; u(h, 0)_{j=1} \; | \; u(2h, 0)_{j=2} \;|\; u(3h, 0)_{j=3}\;|\; u(4h, 0)_{j=4}\;|\; u(5h, 0)_{j=5}\;|\; u(6h, 0)_{j=6 = \text{bound } x=1} \;\mathbf{\blacktriangleleft||\;}\;}_{\text{t=0 discrete space (bar cross-section) datums x=0, j=1, …,j=6(x=1)}}\;$$
Define $U_j^1$ as the temperature evolution after 1 time step,
$$U_j^1 \simeq u(jh, k) $$
Assume the following spatial boundary conditions after a single time-step,
$$U_o^1 = U_M^1 = 0 \quad (\text{peripheral temperature})$$
We now use the finite approximation of the internal $j=1, \cdots, j=5$ temperatures at the future time $t=1$.
$$U_j^1 = U_j^o + \dfrac{k}{h^2} (U_{j+1}^o – 2U_j^o + U_{j-1}^o) \qquad 1 \leqslant j \leqslant M-1$$
Generalisation
The example above can be generalised for any number of time-steps $i$, and discrete time-space increments $j$, provided we have the same initial conditions
$$U_o^{i+1} = U_M^{i+1} = 0 \quad (\text{peripheral temperature})$$
We can determine the temperature at any time $i+1$
\begin{align}U_j^{i+1} &= U_j^i + \mu (U_{j+1}^i – 2U_j^i + U_{j-1}^i) \qquad 1 \leqslant j \leqslant M-1\\ \\
&= \mu U_{j+1}^i + (1-2\mu)U_j^i + \mu U_{j-1}^i
\end{align}
where $\mu = \dfrac{k}{h^2}$.
The Matrix Form of Explicit Euler
The equation for $U_j^{i+1}$ may be expressed in matrix form as
\begin{align}
\mathbf{U^{i+1}} &= T\mathbf{U^i} \\ \\
T &=
\left(\begin{array}{ccccc} 1-2\mu & \mu & & & \\ \mu & 1-2\mu & \mu & \\ & \ddots & \ddots & \ddots & \\ & & \mu & 1-2\mu & \mu \\ & & & & 1-2\mu\end{array}\right) \quad \in \mathbb{R}^{(M-1) \times (M -1)}\\ \\
\mathbf{U^{i}} &= \left(\begin{array}{c}U^i_1 \\ \vdots \\ U^i_{(M-1)}\end{array}\right) \qquad \in \mathbb{R}^{(M-1)}
\end{align}
where $T$ is a tridiagonal symmetric Toeplitz (TST) matrix. In practice $T$ is sometimes simply written as
$$T(1-2\mu, \mu) \quad \in \mathbb{R}^{(M-1) \times (M -1)}$$
Multiple Time-Steps
It is fairly easy to observe that for time step $U^{i+2}$,
\begin{align}
\mathbf{U^{i+2}} &= T\mathbf{U^{i+1}} \\
&= T^2\mathbf{U^i} \end{align}
obtained by simply substituting the prior expression for $U^{i+1}$. Effectively this is the same as
\begin{align}
\mathbf{U^{i}} &= T\mathbf{U^{i-1}} \\
&= T^2\mathbf{U^{i-2}} \end{align}
Now we can continue going back until the $U^o$ temperature time step;
\begin{align}
\mathbf{U^{i}} &= T\mathbf{U^{i-1}} \\
&= T^2\mathbf{U^{i-2}}\\
\mathbf{U^{i}} &= T^i\mathbf{U^{o}} \end{align}
Theorem
- If $\mu \leqslant \frac{1}{2}, \; T \longrightarrow 0 \text{ and } \mathbf{U^{i}} \longrightarrow 0. \textbf{ STABLE}$
- If $\mu \geqslant \frac{1}{2}, \; \mathbf{U^{i}} \longrightarrow \infty. \textbf{ UNSTABLE}$