I'm currenly working on my own PDE solver for non-linear simulations in python. I've done succesfully simulations for KdV and Fisher's equation, but now I'm playing with second order derivatives in time.

The equation I'm using is 1D wave equation, with $\Delta x = 0.5$ and $\Delta t = 0.01$. For simplicity I'm using $c = 1$ and solving by defining an auxiliar field $v(x,t)$, so the equation looks like

$$\partial_t u = v(x,t); \qquad \partial_t v = c^2 \partial_{xx} u$$

Currently I'm using my own RK4 integrator for time and finite differences up to $O(h^2)$ for space. My problem is that as time passes by, with $u(x,0) = \sin(\alpha x)$ and $\partial_t u(x,0) = 0$ the solution seems to grow until it numerically diverges. I tried with gaussian initial conditions, such as $u(x,0) = \exp(-\beta x^2)$, $\partial_t u(x,0) = 0$ and $u(x,0) = 0$, $\partial_t u(x,0) = \exp(-\gamma x^2)$and the solutions seem to irradiate waves which are not part of the analytical solution for these cases.

The boundary conditions are not fixed for the gaussian cases, I just tried these initial condions in order to find what's going on.

Does someone faced the same issues? Could it be for the order on finite differences? I'm pretty new at simulations, so anything would be helpful :)

Here are the plots for the situations mentionatres above, respectively.

enter image description here enter image description here enter image description here

  • $\begingroup$ If you're new to this, the CFL condition is something to be aware of ... en.wikipedia.org/wiki/…. It would be helpful also to know what exactly you are doing at the boundaries to close off the system. $\endgroup$ Commented Mar 11, 2021 at 17:13
  • $\begingroup$ A practical way to quickly verify that your solution is not facing time step related instability is to simply try with a time step twice or ten times smaller and see if the evolution is the same. Alternatively, you can construct a "quasi-exact" solution by using integrators with adaptive time stepping and tight error tolerances (see solve_ivp for instance). If the refined solution looks better, then you most certainly have a stability issue. Though, judging from your graphs, the solution does not explode, so the time step is most likely not the problem. $\endgroup$
    – Laurent90
    Commented Mar 11, 2021 at 20:51
  • $\begingroup$ Good advice from @Laurent90, very direct and practical. A concern I have is about how you’ve implemented boundary conditions since, done incorrectly, can cause issues like energy being reflected back into the region when you expect it to radiate. But using a narrow Gaussian for the initial conditions is probably good enough. $\endgroup$ Commented Mar 12, 2021 at 0:46
  • 1
    $\begingroup$ Thank you for the suggestions, finally it was a problem with de boundary conditions and now I have done simulations for my nonlinear case (Parametrically Driven Nonlinear Schrödinger Equation) with the CFL condition in mind and works just fine. $\endgroup$ Commented Mar 15, 2021 at 16:40
  • $\begingroup$ Good luck @RafaelRiverosÁvila! $\endgroup$ Commented Mar 15, 2021 at 16:41


Your Answer

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