I want to integrate this equation in time:
$$\frac{\partial u(x,t)}{\partial t} = \frac{\partial}{\partial x} \kappa(x) \frac{\partial u(x,t)}{\partial x}$$
with initial condition
$$u(x,0) = \frac{1}{2} \left [\text{erf}\left ( \frac{x-L/4}{s} \right ) - \text{erf}\left ( \frac{x+L/4}{s} \right ) \right ]$$
over the domain $x \in [-L/2,L/2]$. The diffusivity is given by
$$\kappa(x) = \frac{1}{\sqrt{2 \pi} \sigma} \left [ \text{exp} \left [ - \frac{(x+L/4)^2}{2 \sigma^2}\right ] +\text{exp} \left [ - \frac{(x-L/4)^2}{2 \sigma^2}\right ] \right ].$$
I use $L=2$ and $s=1/20$, so the initial condition looks like this:
For the diffusivity I use $\sigma = 1/20$, so the profile looks like this:
Since the problem is periodic over the length $L$, I compute derivatives using the discrete Fourier transform (DFT). I found a description for how to do this at math.mit.edu/~stevenj/fft-deriv.pdf, and I use Algorithm 3 suggested in this document. I do not understand the full background, but it is written that in order to have a self-adjoint operator the algorithm consists of:
Compute $\frac{\partial u}{\partial x}$ the conventional spectral way. This means that the wave numbers are: $$k = \frac{2\pi}{L}(0, 1, 2, \dots , N/2-2, N/2-1, 0, 1-N/2, 2-N/2 , \dots , -2, -1),$$ where $N$ is the number of grid points (I use $N=10^4$), and the derivative is: $$\frac{\partial u}{\partial x} = \mathcal{F}^{-1} i k \mathcal{F} u$$ at each grid point. $\mathcal{F}$ and $\mathcal{F}^{-1}$ is the forward and inverse DFT. Before applying the inverse transform I save the coefficient of the Nyquist frequency $\hat{u}_{N/2}$.
Compute $\kappa \frac{\partial u}{\partial x}$ at each grid point.
Compute $\frac{\partial }{\partial x} \kappa \frac{\partial u}{\partial x}$ the conventional spectral way with this modification: Before applying the inverse transform I set the Nyquist frequency to $$-\langle \kappa \rangle \left (\frac{N \pi}{L} \right )^2 \hat{u}_{N/2},$$ where $\langle \kappa \rangle$ is the mean value of $\kappa$.
Now that I can compute $\frac{\partial u}{\partial t}$, I integrate forward in time. However, here I am doing domething wrong. For the first three time steps using $\Delta t = 10^{-6}$, nothing happens. But at the fourth time step, huge oscillations appear, as seen in this in-zoom (compare with the intial condition above, that looks identical to $u(x,3\Delta t)$):
For the time integtration, I have tried both one-step forward Euler and 4:th order Runge-Kutta. Both methods give these oscillations.
What am I doing wrong? How can one integrate the heat equation with a varying diffusivity numerically?
