I've been trying to find an explicit method for the 1D heat-diffusion equation with variable thermal conductivity coefficient $\kappa$:

$$\frac{\partial T}{\partial t}=\frac{\partial}{\partial x} (\kappa \frac{\partial T}{\partial x})$$ However, it seems to me like there is no simple approach as in the case where $\kappa$=const. The explicit methods I found were quite complicated and to be honest, I haven't fully understood them while implicit methods look very similar to the linear ($\kappa$=const.) case.

Why is it so hard to use an explicit method in the nonlinear case? What is the problem?


  • $\begingroup$ You said that $\kappa$ is variable. But variable as a function of what ? Function of $x$? of $t$? of both ? $\endgroup$ – JJacquelin Oct 2 '16 at 5:46
  • $\begingroup$ As a function of x $\endgroup$ – OD IUM Oct 2 '16 at 11:26

When $k=k(x)$ is not constant, the PDE is still linear, just no longer constant coefficient. There are generally 2 approaches to explicit finite difference schemes in this case.

1) Expand the right hand side and write the PDE as

$\frac{\partial T}{\partial t} = k\frac{\partial^2 T}{\partial x^2} + k'\frac{\partial T}{\partial x}.$

Use the standard discretization for the first term, and use centered differences for $\frac{\partial T}{\partial x}$. It sounds like you already know how to solve the constant coefficient heat equation, so you should be good from here. If you need more details I can expand on this.

2) Discretize the equation in divergence form. I'll give more details here. Choose a time step $\Delta t$ and a spatial step $\Delta x$, and write $T^n_j = T(j\Delta x,n\Delta t)$ and $k_j = k(j\Delta x)$. Basically we will use a centered difference for the outer $\frac{\partial }{\partial x}$ using the mid-grid points $j+1/2$ and $j-1/2$. The scheme is

$T^{n+1}_j = T^n_j + \frac{\Delta t}{\Delta x^2} \left( k_{j+1/2}(T^n_{j+1} - T^n_j) - k_{j-1/2}(T^n_j - T^n_{j-1})\right).$

Since $k$ is a given function, there is no difficulty in evaluating $k$ off the grid.

EDIT: Let me add a few words about stability. Normally one does a Von Neumann stability analysis. For parabolic equations, there is a short cut based on monotonicity. For example, for method (2) if we write $s=\Delta t/\Delta x^2$ then

$T^{n+1}_j = (1-s(k_{j+1/2} + k_{j-1/2}))T^n_j + s(k_{j+1/2}T^n_{j+1} + k_{j-1/2}T^n_{j-1}).$

The scheme is monotone if all the coefficients on the right hand side (the numbers in front of $T^n_j$, $T^n_{j-1}$ and $T^n_{j+1}$) are all positive. Then the scheme satisfies a maximum principle and you can prove stability. Since $k$ is the thermal conductivity, $k>0$, so the only term in question is the first one. Hence the scheme is stable when

$1-s(k_{j+1/2} + k_{j-1/2}) \geq 0.$

Rearranging, we require

$\frac{\Delta t}{\Delta x^2} = s \leq \frac{1}{k_{j+1/2} + k_{j-1/2}}.$

If the PDE is uniformly elliptic, then there exists $\theta>0$ so that $k(x)\geq \theta$ for all $x$. Hence, the PDE is stable when $s \leq 1/(2\theta)$, or

$\Delta t \leq \frac{1}{2\theta} \Delta x^2.$

You can perform the same analysis for scheme (1), and you will get the same time step constraint as above, with the additional constraint that

$\Delta x \leq \min_x \frac{k(x)}{|k'(x)|}.$

If we let $M = \max_x |k'(x)|$, then we can also phrase this stability constraint as

$$\Delta x \leq \frac{\theta}{M}.$$

If $k$ is constant, we have $M=0$ and there is no constraint on $\Delta x$. In this respect, perhaps scheme (2) is slightly better, though the constraint on $\Delta x$ above is not burdensome, unless the PDE is degenerate $\theta =0$ or singular $k'=\infty$.

  • $\begingroup$ Thx for your reply. Is there a name for the 1st approach so that I can Google for it? Never seen that before. And which of these methods is more stable? $\endgroup$ – OD IUM Oct 4 '16 at 8:14
  • 1
    $\begingroup$ I don't think there is a specific name for the first scheme. Any book on numerical schemes for PDE will describe both. I've added a short summary of stability for each scheme. There is not much difference in terms of stability, except that (2) requires an additional constraint on $\Delta x$, but this is not a very strict constraint. $\endgroup$ – Jeff Oct 4 '16 at 16:21

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.