Boundary value technique for heat equation

Question

My heat equation is $$ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}, \quad x \in [0,1], \quad t \in (0,0.1] $$ with initial condition $u(x,0)=\sin(\pi x)$ and homogeneous Dirichlet boundary conditions $u(0,t) = u(1,t) = 0$. The exact solution is $u(x,t) = e^{-\pi^2t}\sin(\pi x)$. I approximate the equation with forward in time numerical scheme $$ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{u_{i+1}^n - 2 u_i^n + u_{i-1}^n}{(\Delta x)^2} $$ I would like to use boundary value technique to numerically solve this equation, i.e. i would like to write the following system of $MN$ unknowns in matrix form and solve it in Matlab.

I partition time as $0 = t_0 < t_1 < \dots < t_{N+1} = 0.1$. Step size is of length $h=\Delta t =\frac{0.1}{N+1}$, i.e. $t_n = t_0 + nh$, $n=1,2,\dots,N+1$.

I partition space as $0 = x_0 < x_1 < \dots < x_{M+1} = 1$. Step size is of length $h'=\Delta x =\frac{1}{M+1}$, i.e. $x_i = x_0 + ih'$, $i=1,2,\dots,M+1$.

The system to solve is then \begin{align*} \frac{u_i^{n+1} - u_i^n}{\Delta t} &= \frac{u_{i+1}^n - 2 u_i^n + u_{i-1}^n}{(\Delta x)^2}, \quad i=1,2,\dots,M, \quad n=0,1,2,\dots,N \\ u_0^n &= u_M^n = 0, \quad n=0,1,\dots,N+1 \\ u_i^0 &= \sin((i\Delta x) x), \quad i=0,1,\dots,M+1 \end{align*} I have big problems how to write this in matrix form. Any suggestions?

Answer 1

Let $$A = \left(\begin{array}{ccccc} -2 & 1 \\ 1 & -2 & 1 \\ &&\ddots\\ &&1 & -2 & 1\\&&&1&-2\end{array}\right)$$ be the standard Laplacian matrix. Your system can be written in traditional time-stepping form as $$u^{n+1} = \left(I + \frac{h}{(h')^2}A\right)u^n \doteq Mu^n,$$ where $$u^n = \left(\begin{array}{c} u^n_1 \\ \vdots \\u^n_M\end{array}\right).$$ Now define $$u = \left(\begin{array}{c} u^1 \\ \vdots \\ u^{N}\end{array}\right),$$ and $$B = \left(\begin{array}{ccccc} I \\ -M & I \\ &-M &I\\&&&\ddots\\&&&-M &I\end{array}\right)$$ and $$b = \left(\begin{array}{c} Mu^0 \\ 0 \\\vdots\\0\end{array}\right).$$ The solution for all time $u$ is then determined by the solution to the linear system $$Bu = b.$$

---

It's worth noting that this is a highly inefficient way to solve the heat equation. A naive approach to solving this system would require $\mathcal{O}(N^3M^3)$ computations. In fact, any method short of multigrid would require significantly more than $\mathcal{O}(NM)$.

Conversely, formulating the system with time-stepping results in $\mathcal{O}(NM)$ computations using the relatively naive approach of applying the stencil pointwise at each time step.