Solve a fourth order differential equation

I want to solve $$\frac{\partial^2}{\partial t^2}u(z,t) + a\frac{\partial^2}{\partial z^2}u(z,t) + k\frac{\partial^4}{\partial z^4}u(z,t) = 0$$

with $$u(z,0) = 1+0.1e^{-\frac{z^2}{2}}$$. I'd like to study the dynamic of the wave.

What is the simplest way to do it?

I tried with the finite difference method but my solution diverges very quickly!

I use this iteration

$$f_j^{n+1} = 2f_j^n - f_j^{n-1} + \frac{\Delta t^2}{\Delta z^2}a(f_{j+1}^n - 2f_j^{n} + f_{j-1}^n) +k \frac{\Delta t^2}{\Delta z^4}(f^n_{j-2}-4f^n_{j-1}+6f^n_{j}-4f^n_{j+1} + f^n_{j+2})$$ where $$n$$ it's the time step and $$j$$ the space step.

Improving my question after comments:

I used a very small time step and this is what I obtain: How can avoid the reflection with this iteration method?

• What are your boundary conditions? – Jesse Chan Sep 8 '14 at 19:55
• OK - they probably aren't the issue on second glance. You may just need to use an incredibly small timestep? A rough guess is that fourth order equations imply a CFL condition of $\Delta t < C\frac{\Delta z^2}{\sqrt{k}}$ for some constant $C$. – Jesse Chan Sep 8 '14 at 22:46
• I believe for a 4th order equation it's even worse, $\Delta t < \mathcal{O}(\Delta z^4)$. Either use a REALLY small time step or an implicit method. – Doug Lipinski Sep 9 '14 at 0:50
• Reflections are likely a result of your boundary conditions - can you specify them? – Jesse Chan Sep 9 '14 at 13:51
• @JLC the script O is just \mathcal{O} – Doug Lipinski Sep 10 '14 at 2:05

You haven't specified boundary conditions, but it sounds like you want to solve the Cauchy problem (i.e., the initial value problem on the whole real line $\infty < z < \infty$).

In that case, there is no need to use finite differences. You can simply use Fourier analysis (separation of variables) to write down the exact answer. Once you have that, you can evaluate it numerically by truncating the domain and the set of Fourier modes.

If you insert the ansatz

$$u(z,t) = \exp(i(\xi x - \omega t))$$

into your wave equation, you find there are two frequencies of propagation for each wavenumber $\xi$:

$$\omega(\xi) = \pm \sqrt{a \xi^2 - k \xi^4}.$$

These correspond to left- and right-going waves. To write down the full solution you need to determine which part of the initial data goes each way. That depends on $u_t(z,0)$, which you haven't specified. Without that information, the problem is not well-posed.