Numerical solution of non-constant coefficient diffusion equation via finite-difference method

I would like to solve numerically a diffusion equation: $$\frac{\partial c}{\partial t} = \frac{\partial}{\partial x}\left(k\frac{\partial c}{\partial x}\right),$$ given some boundary conditions and initial value.

My idea is to solve it via discretization of the spatial variable using finite-difference method and then solve the resulting system of ODEs with Matlab, that is, use the Method of Lines approach. However, as the term $k=k(x)$ then I do not know how to discretize the right hand side of the equation. Could you suggest any idea?

Thanks!

If $k$ is continuously differentiable, then there is a simple approach:

$$\frac{\partial}{\partial x} \left ( k \frac{\partial c}{\partial x} \right ) = k'(x) \frac{\partial c}{\partial x} + k(x) \frac{\partial^2 c}{\partial x^2}.$$

Now you can discretize these two spatial derivatives as you see fit.

More generally, it is better to think about fluxes. (Indeed it is essential to think about fluxes in hyperbolic problems like the transport and wave equations, even though it is not so crucial in parabolic problems like the diffusion equation.) Consider a point $x$ on the line with neighbors $x-\Delta x$ and $x+\Delta x$. The flux balance at this point over a time step of length $\Delta t$ should look like this:

\begin{align} x-\Delta x \to x & : k(x-\Delta x/2) \Delta t c(x-\Delta x) \\ x \to x-\Delta x & : k(x-\Delta x/2) \Delta t c(x) \\ x \to x+\Delta x & : k(x+\Delta x/2) \Delta t c(x) \\ x+\Delta x \to x & : k(x+\Delta x/2) \Delta t c(x+\Delta x)\end{align}

(plus corrections of higher order, of course).

This would give the ODE form

$$\frac{\partial c}{\partial t}(t,x)=k(x-\Delta x/2)c(t,x-\Delta x)+k(x+\Delta x/2)c(t,x+\Delta x) \\ -c(t,x)(k(x-\Delta x/2)+k(x+\Delta x/2)).$$

If $k$ is continuous then this will work, though it could be slow if $k$ is especially irregular. If $k$ is not even continuous, then you must be even more careful than this, and will need to write down a weak formulation of the problem in order to write down a sensible method.

• The scheme you obtained using the flux balance method does not seem to match the scheme obtained from doing the chain rule and discretising. However, the discretising after the chain rule does match with Julian Aguirre's answer. Could you please explain why you have difference between the scheme from the flux balance method and Julian Aguirre's scheme? – adamG Jun 11 '17 at 19:12
• @adamG I see looking now that I should have been measuring the conductivity at the interface not the start node, i.e. the relevant values are really $k(x-\Delta x/2)$ and $k(x+\Delta x/2)$. I will fix it later. I think after the correction the methods will agree up to a higher order correction. – Ian Jun 11 '17 at 19:52
• That makes sense. Thanks for correcting! – adamG Jun 12 '17 at 19:59
• @adamG See e.g. this post for complements. – EditPiAf Apr 15 '19 at 21:34
• How would it look like if the PDE was instead $$\frac{\partial c}{\partial t} = k(t,x)\left(\frac{\partial^2 c}{\partial x^2}\right) ?$$ – Parseval Feb 12 at 21:06

The straight forward way would be to approximate the right hand side by $$\frac1h\Bigl(k(x+h)\frac{c(x+2h)-c(x+h)}{h}-k(x)\frac{c(x+h)-c(x)}{h}\Bigr).$$

• Thanks for the answer. However, I do not understand how did you manage to get there. My idea was to do the chain rule and then discretize so I get an expression divided by $h^2$. – DOMiguel Nov 2 '15 at 12:53
• I forgot to divide by $h$. Approximate first the outer derivative and then the inner one. – Julián Aguirre Nov 2 '15 at 15:34
• Ok, eskerrik asko. – DOMiguel Nov 2 '15 at 20:22
• @DOMiguel Alternatively, one could write the r.h.s. as $\frac{1}{h}\left(f(x+h/2) - f(x-h/2)\right)$, where $$f(x+h/2) = k(x+h/2) \frac{c(x+h)-c(x)}{h} \, ,$$ which is spatially symmetric. – EditPiAf Dec 6 '17 at 14:37