# Stability in Numerical Heat Equation

The heat equation $u_t=u_{xx}$ in $\Omega = (0,T] \times(0,1)$, $u(0,x) = u_0(x)$, $u(t,0)=g_1(t)$, $u(t,1)=g_2(t)$ can be discretized by the schemes:

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} = \theta\frac{u_{i-1}^{n+1}-2u_i^{n+1}+u_{i+1}^{n+1}}{\Delta x^2} + (1-\theta)\frac{u_{i-1}^n-2u_i^n+u_{i+1}^n}{\Delta x^2}$$

where $u_i^0=u_0(x_i), u_0^n=g_1(t_ n), u_N^n=g_2(t_n)$, for $\theta \in [0,1]$.

How can we show that the order of consistency is $O(\Delta t + \Delta x^2)$ for $\theta \neq \frac{1}{2}$ and $O(\Delta t^2 + \Delta x^2)$ for $\Theta = \frac{1}{2}$ ? For which $\theta$ is the scheme stable?

• How can I do this Taylor expansion? There are so many terms, I'm confused on how to do that. Commented May 17, 2018 at 19:08

Here are some pointers on what you need to do.

Let's call $u_i^n = u(x,t)$. Then

$$u_{i+i'}^{n+n'} = u(x+i'\Delta x,t + n'\Delta t)$$

where $\Delta x,\Delta t$ is the grid-spacing (for example $u_{i+1}^n$ corresponds to $i'=1$ and $n'=0$). We can Taylor expand $u(x+i'\Delta x,t + n'\Delta t)$ (see this for formulas) $$u(x+i'\Delta x,t + n'\Delta t) = u(x,t) + (i'\Delta x) u_x(x,t) + (n'\Delta t) u_t(x,t) + \frac{1}{2}[u_{xx}(x,t)(i'\Delta x)^2+2u_{xy}(x,t)(i'\Delta x)(n'\Delta t)+u_{tt}(x,t)(n'\Delta t)^2] + \ldots$$

and so on to as high order as you need. Now start with your discretized PDE and insert the Taylor expansion for each of the terms and simplify. You will end up with $u_{t}(x,t) = u_{xx}(x,t) +$error-terms where the error terms will tell you the order of accuracy.

As an example we have the left hand side is given by:

$$\frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{[u + u_t\Delta t + \frac{1}{2}u_{tt}\Delta t^2] - [u]}{\Delta t} = u_t + \frac{u_{tt}}{2}\Delta t$$ which means this expression is accurate to $\mathcal{O}(\Delta t)$.

As for stability see Von Neumann stability analysis. The example given on that page is for the heat equation with the discretization corresponding to $\theta = 0$. Should be straight forward to modify this analysis to your general case.

• I tried to do this but I can't see how the order of consistency can change with the value of $\theta$! It seems to me that is should just be $O(\Delta t + \Delta x^2)$! Commented May 18, 2018 at 22:53
• @Pdqueijo The terms on the left and right hand side can cancel. The leading order error term on the left is $u_{tt} \Delta t$ and we should have $\propto u_{xxt}\Delta t$ error term coming from the first term on the right hand side expanded to third order. Since $u_{xxt} = \partial_t [u_{xx}] = \partial_t [u_t] = u_{tt}$ when $u$ is a solution to the PDE these can cancel! I suspect this is what happens Commented May 18, 2018 at 23:11
• ...and this is indeed what happens. You will get $u_t - u_{xx} \propto \frac{u_{tt}\Delta t}{2} - \theta u_{xxt} \Delta t + \mathcal{O}(\Delta t)^2 + \mathcal{O}(\Delta x)^2$ and the two first terms will conspire when $\theta = 1/2$ since $u_{xxt} = u_{tt}$ Commented May 18, 2018 at 23:27