explicit finite difference scheme for nonlinear diffusion

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?

regards.

• You said that $\kappa$ is variable. But variable as a function of what ? Function of $x$? of $t$? of both ? – JJacquelin Oct 2 '16 at 5:46
• As a function of x – 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$.

• 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? – OD IUM Oct 4 '16 at 8:14
• 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. – Jeff Oct 4 '16 at 16:21