I am a physicist and for my work I have to solve a non-linear, one dimensional, diffusion problem numerically.

So the problem can be formulated simply as

$$ \frac{\partial x}{\partial t} = \frac{\partial }{\partial z} D(x(z,t),t) \frac{\partial x}{\partial z} = - \frac{\partial}{\partial z} J(x(z,t),t) $$

$x$ is for my purpose some conserved and bounded quantity. Let's say $x \in (0,1)$. Furthermore, let's assume that our system is a box with non-permeable walls. That is $z \in [0,d]$ with $d>0$. We will have some, not necessarily everywhere continuous, initial $x_0(z) = x(z,0)$ function. The non-permeabilty of the walls can be expressed by the no-flow boundary conditions. That is $J|_{z=0}=0=J|_{z=d}$

In what I have to do, the function $D$ may change a lot, so I would like to keep this part as general as possible. In the simplest case, that is in my problem the interdiffusion of two species, the function may look like


( that is why $x \neq 0$ and $x\neq 1$) But as already noted, I would like to discuss this problem for somewhat more general functions $D(x,t)$. Perhaps be $D \in C^{\infty}$ for $x\in(0,1)$ and $t\in R$.

Since I am not an expert for numerical problems, I first just tried several approaches to solve this problem. Among them Crank-Nicholson, finite volume methods and Lax-Friedrichs method. The first two led to the development of strong oscillations that formed on the boundaries and propagated through the whole system. The Lax-Friedrich method led to the problem that for different time steps I chose for the simulation, the resulting curves differed. That is if I chose a time step $\Delta t$ and a time step $\frac{\Delta t}{2}$, the result was not the same.

So, my question is, what is the right numerical approach to solve such a problem? Is there a general approach or does it depend explicitly on the form of the function $D$? Perhaps someone can elaborate on how to find a fitting approach and how to test its applicability?

Thanks in advance :)


For diffusion type problems you want to use an implicit time stepping method, such as Crank - Nicholson. An explicit method will not have suitable convergence properties for your equation, in particular the Lax-Friedricks method requires $\Delta t < c \Delta x^2$ for some constant $c$ to be convergent, which is why you get different results for different temporal resolutions, the time steps are too big in both cases.

As to why the Crank-Nicolson approach breaks, you say that oscillations are generated near the boundary, which indicates that your boundary conditions are not being applied properly. I would recommend sticking with this approach and trying to figure out what you've done wrong, or posting another question detailing your approach and seeing if anyone has any feedback.


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.