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.


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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.