Numerical Solution to the Initial Value Problem for the Heat Equation Given Dirichlet Boundary Conditions 📂Numerical Analysis

Numerical Solution to the Initial Value Problem for the Heat Equation Given Dirichlet Boundary Conditions

Example ¹

$$ \begin{cases} u_{t} = \gamma u_{xx} \\ u(t,0) = u(t,l) = 0 \\ u(0,x) = f(x) \end{cases} $$

The given problem has an algebraic solution simple enough to solve, yet it serves as a clear example of why we learn numerical methods to solve differential equations. It shows how solving a simple differential equation of the form $y ' = f(x,y)$ leads to the solution of partial differential equations.

Solution

Unlike in algebraic solutions, the value of $\gamma$ is not important, so let’s say it is $\gamma = 1$ for convenience.

$$ x_{j} = \delta j $$ If we set the step size of the Line Space $x$ to $\delta$, then for $j= 0, 1, \cdots , m$, it can be represented as shown above.

For the second order derivative approximation of $x$ $$ U_{xx} \approx {{ U(t, x_{j+1}) - 2 U(t, x_{j}) + U(t, x_{j-1}) } \over { \delta^2 }} $$ then only for the equation regarding $t$ $$ u_{t} (t, x_{j}) = {{ u(t, x_{j+1}) - 2 u(t, x_{j}) + u(t, x_{j-1}) } \over { \delta^2 }} $$ we can obtain $m-1$ equations. Representing this as a matrix gives us $$ \begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix} = {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix} \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix} $$ and to save space, $$ \mathbf{u}’(t) :=\begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix} $$

$$ \Lambda := {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix} $$

$$ \mathbf{u} (t) := \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix} $$ let’s say it is $\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)$. Now, setting the time step size $t$ to $h$ and calling it $t_{k} = hk$, we can apply various methods to find the numerical solution.

■

Lesson: Stiff Differential Equations

However, all these are theoretical discussions, and actually finding a numerical solution poses significant challenges.

To know the eigenvalues of $\Lambda$, let’s calculate the determinant of matrix $\Lambda - \lambda I_{m-1}$.

Determinant of Tridiagonal Matrices: For the tridiagonal matrix $X_{n} := \begin{bmatrix} x & 1 & 0 & \cdots & 0 & 0 \\ 1 & x & 1 & \cdots & 0 & 0 \\ 0 & 1 & x & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & x & 1 \\ 0 & 0 & 0 & \cdots & 1 & x \end{bmatrix}$, $\displaystyle | x_{n}| = U_{n} \left( {{x} \over {2}} \right)$

$$ \Lambda - \lambda I_{m-1} = {{1} \over { \delta^2 }} \begin{bmatrix} -2 - \delta^2 \lambda & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 - \delta^2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 - \delta^2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 - \delta^2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 - \delta^2 \end{bmatrix} $$ since it is a tridiagonal matrix $(m - 1) \times (m-1)$, the determinant is calculated as $$ \det \left( \Lambda - \lambda I_{m-1} \right) = {{1} \over { \delta^2 }} U_{m-1} \left( {{-2 - \delta^2 \lambda } \over {2}} \right) $$. Here, $U_{m-1}$ means the second kind Chebyshev polynomials of order $m-1$.

Roots and Symmetry of the Second Kind Chebyshev Polynomials:
[1]: Roots of $\displaystyle U_{n} (X)$ are $\displaystyle x_{k} = \cos \left( {{k} \over {n+1}} \pi \right)$, $k=1, \cdots , n$
[2]: $U_{n} (-x) = (-1)^{n} U_{n} (x)$

Extracting $-1$ from the bracket gives us $$ {{(-1)^{m-1}} \over { \delta^2 }} U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right) $$ and to find the eigenvalues, only $$ U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right) $$ needs attention. The $j$th root of $\displaystyle U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right)$ for $j=1, \cdots , (m-1)$ is $$ {{2 + \delta^2 \lambda_{j}} \over {2}} = \cos \left( {{j} \over {m}} \pi \right) $$ and multiplying both sides by $2$ gives $$ \begin{align*} \displaystyle 2+\delta^2 \lambda_{j} =& 2 \cos \left( {{j} \over {m}} \pi \right) \\ =& 2 \cos \left( 2 \cdot {{j} \over {2m}} \pi \right) \\ =& 2 \left[ 1 - 2 \sin^2 \left( {{ j \pi } \over { 2m }} \right) \right] \end{align*} $$ Therefore, the eigenvalues are obtained as $$ \lambda_{j} = - {{4} \over { \delta^2 }} \sin^2 \left( {{j \pi } \over {2m}} \right) $$ However, $$ \begin{align*} \left| \lambda_{m-1} \right| =& \left| {{4} \over { \delta^2 }} \sin^2 \left( {{(m-1) \pi } \over {2m}} \right) \right| \\ \approx & {{4} \over {\delta^2}} \\ =& 4 m^2 \end{align*} $$

thus, as $\displaystyle m = {{1} \over {\delta}}$ gets larger, the magnitudes of the eigenvalues also hugely increase. If you chose a small $\delta$ to split $[0,1]$ into many parts to accurately know the numerical solution, it becomes exceedingly challenging to compute. Even if you set it to $m=100$, $\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)$ becomes a Stiff problem, making it quite difficult to solve numerically.