# 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? Jun 11, 2017 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, 2017 at 19:52
• That makes sense. Thanks for correcting! Jun 12, 2017 at 19:59
• @adamG See e.g. this post for complements. Apr 15, 2019 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) ?$$ Feb 12, 2021 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$. Nov 2, 2015 at 12:53
• I forgot to divide by $h$. Approximate first the outer derivative and then the inner one. Nov 2, 2015 at 15:34
• Ok, eskerrik asko. Nov 2, 2015 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. Dec 6, 2017 at 14:37