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:

The initial condition: $u(x,0)$

For the diffusivity I use $\sigma = 1/20$, so the profile looks like this:

The diffusivity $\kappa(x)$

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:

  1. 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}$.

  2. Compute $\kappa \frac{\partial u}{\partial x}$ at each grid point.

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

$u(x,t)$ after four time steps.

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?



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.