Need Help Discretizing this Heat Equation

Question

This might be a fairly easy question but since I haven't done much numerical PDE's before I don't really know what to do. I know that for a normal heat equation i.e. $\large \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}$ (for $D$ arbitrary) we are able to discretize this equation to obtain an implicit numerical form:

$\Large \frac{u^{k+1}_{j} - u^{k}_{j}}{\Delta t} = D \frac{u^{k+1}_{j+1} - 2u^{k+1}_{j} + u^{k+1}_{j-1}}{\Delta x^2}$

which can then be rewritten as

$\large u^k_j = (1 + 2S)u^{k+1}_{j} - S (u^{k+1}_{j+1} + u^{k+1}_{j-1})$ where $S = \frac{D\Delta t}{\Delta x^{2}}$

and we can then use this to solve the numerical PDE through the use of a tridiagonal matrix and by solving a matrix system.

However, say for instance I now have a modified heat equation such as:

${\Large \frac{\partial u}{\partial t} = D \frac{1}{ \frac{\partial^2 u}{\partial x^2}} - D}$

If I try and discretize this equation, I now get

$\Large \frac{u^{k+1}_{j} - u^{k}_{j}}{\Delta t} =D \frac{1}{\frac{u^{k+1}_{j+1} - 2u^{k+1}_{j} + u^{k+1}_{j-1}}{\Delta x^2}}$

Following the same method as I wrote above, I now get:

$\Large u^k_j = u^{k+1}_j - D \frac{1}{\frac{u^{k+1}_{j+1} - 2u^{k+1}_{j} + u^{k+1}_{j-1}}{\Delta x^2}}$

The problem I'm having now however is that I'm not sure how to write this into a program since I'm not able to split up the $u$ terms in a neat way to allow the computer to work with the matrix system and there doesn't seem like there's a way to create a tridiagonal matrix. Did I write this the correct way?

Thanks in advance.

Answer 1

That's quite a strange PDE. I'd rather not call it a heat equation. Do you really need an implicit discretization for it? Anyway, you've obtained a set of nonlinear algebraic equations for $u^{k+1}_j$ and may utilize some nonlinear solver, like Newton's method.

Note that the system of equations is still tridiagonal, that is $$ \left\{ \begin{gathered} &f_1(u_1, u_2) = u_1^{k}\\ &f_2(u_1, u_2, u_3) = u_2^{k}\\ &f_3(u_2, u_3, u_4) = u_3^{k}\\ &\vdots\\ &f_j(u_{j-1}, u_{j}, u_{j+1}) = u_j^k\\ &\vdots\\ &f_{n-1}(u_{n-2}, u_{n-1}, u_{n}) = u_{n-1}^{k}\\ &f_{n}(u_{n-1}, u_{n}) = u_{n}^{k} \end{gathered} \right. $$ I've removed $k+1$ index for simplicity.

The Newton's method linearizes this system around some approximate solution $u_j^{(s)}$. To start with, one can take $u^{k}_j$ as starting $s = 0$ approximation $u_j^{(0)}$.

Let $u^{(s+1)}_j = u^{(s)}_j + \delta u_j$. Then linearized equation $$ f_j(u_{j-1}, u_j, u_{j+1}) = u^k_j $$ becomes $$ f_j(u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) + \frac{\partial f_j}{\partial u_{j-1}} (u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) \delta u_{j-1} + \frac{\partial f_j}{\partial u_{j}} (u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) \delta u_{j} + \frac{\partial f_j}{\partial u_{j+1}} (u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) \delta u_{j+1} = u_j^k. $$ Assembling that altogether we get a linear system for $\boldsymbol \delta \mathbf u$ $$ \mathbf f^{(s)} + J^{(s)} \boldsymbol \delta \mathbf u = \mathbf u^{k} $$ with $J$ being a tridiagonal Jacobi matrix $$ J = \frac{\partial \mathbf f}{\partial \mathbf u}. $$

So the time stepping from $k$ to $k+1$ algorithm becomes:

1. Use $\mathbf u^k$ as initial $s = 0$ approximation for $\mathbf u^{k+1}$
2. Compute $\mathbf f^{(s)}$ Jacobi matrix $J^{(s)}$ using last know approximation
3. Solve linear system for $\boldsymbol \delta \mathbf u$
4. Use $\mathbf u^{(s+1)} = \mathbf u^{(s)} + \delta \mathbf u$ as new approximation and increase $s$. If $\boldsymbol \delta \mathbf u$ is sufficiently small, use $\mathbf u^{(s+1)}$ as $\mathbf u^{k+1}$, else continue from step 2.