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
$$D=\frac{1}{x}+\frac{1}{1-x}$$
( 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.