1

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!

DOMiguel
  • 514

2 Answers2

2

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). $$

Julián Aguirre
  • 77,704
  • 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 02 '15 at 12:53
  • I forgot to divide by $h$. Approximate first the outer derivative and then the inner one. – Julián Aguirre Nov 02 '15 at 15:34
  • 1
    Ok, eskerrik asko. – DOMiguel Nov 02 '15 at 20:22
  • 1
    @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 06 '17 at 14:37
2

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.

Ian
  • 104,572
  • 1
    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? – Neurobro 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! – Neurobro 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 '21 at 21:06
  • @Parseval Then all the $k$'s are just evaluated at $x$ (at least in the leading order approximation). – Ian Feb 12 '21 at 21:20
  • So there would not be any $\Delta t$ and $Delta x$ there? Or what do you mean evaluated at $x$? – Parseval Feb 12 '21 at 22:18
  • @Parseval The appearances of $k$ would all be $k(x)$ instead of $k(x \pm \Delta x/2)$. – Ian Feb 12 '21 at 22:19
  • So what happens to the $t$? K is a function of both $t$ and $x$. – Parseval Feb 12 '21 at 22:20
  • @Parseval In this question it wasn't, actually, but if it were then it would be getting evaluated at $t$ all the time either way. This type of time resolution is called the "method of lines", where you perform a spatial discretization to obtain an ODE system that you then consider to be a "black box" rather than tailoring the space and time discretizations to one another. – Ian Feb 12 '21 at 22:21
  • @Ian - So If I understand you correctly, there is no difference if $k$ is a constant or if $k$ is a function of $t,x$? For example in this scheme: https://flothesof.github.io/heat-equation-cook-my-meat.html Can I just replace $D$ by $k(t,x)?$ – Parseval Feb 12 '21 at 22:24
  • @Ian Yes, I know. I kind of posted my own question here but it has no answers yet: https://math.stackexchange.com/questions/4023602/discretization-scheme-of-heat-equation-with-non-constant-diffusion-coefficient – Parseval Feb 12 '21 at 22:25
  • @Parseval You can basically do the same thing but of course your $k$ will need to change as you change which point you are computing the time derivative at. You just don't need the values $k(x \pm \Delta x/2)$. That being said, the form in the OP (divergence form) is what typically appears in physics. – Ian Feb 12 '21 at 22:27