I have been trying to numerically solve the 1d homogeneous wave equation

\begin{align} u_{tt}(x,t) &=& u_{xx}(x,t) \\ u(x,0) &=& f(x) \\ u_t(x, 0) &=& 0 \\ u(0, t) &=& u(1, t) = 0 \end{align}

by casting it into two first-order linear equations

\begin{align} u_t &=& v_x \\ v_t &=& u_x \end{align}

This is what I've tried

  1. Apply forward finite differences for the $t$ derivatives and centered differences for the $x$ derivatives

\begin{align} u^{n+1}_j &=& u^n_j + \frac{\alpha}{2}(v^n_{j + 1} - v^n_{j - 1}) \\ v^{n+1}_j &=& v^n_j + \frac{\alpha}{2}(u^n_{j + 1} - u^n_{j - 1}) \end{align}

where $u^n_j = u(t_n, x_j)$ and $\alpha = \Delta t/\Delta x < 1$. As initial condition $f(x)$ I used a gaussian function of variance 0.1. This is what I get when the pulse approaches the boundaries:

enter image description here Colors label different times, until the wave hits the boundaries, were it gets really bad!

  1. I also followed this document

Lecture notes from Prof. L. Rezzolla

Turns out he introduces yet another variable I don't see the need for. See his Chapter 5. It sort of works here, but still the wave gets really distorted when it reaches the edges of the $x$-domain

  1. Upgraded the order of the finite differences for the $x$-variable, it does not improve.

I am out of ideas, any suggestion is much appreciated

  • $\begingroup$ Please describe, how you deal with $v_0^{n+1}$ and $v_J^{n+1}$ (boundary values), that's crucial $\endgroup$ – uranix Nov 1 '16 at 21:57

I would say the easiest way to discretize the wave equation is to use centered differences for each second derivative, to get

$\frac{u^{n+1}_j - 2u^n_j + u^{n-1}_j}{\Delta t^2} = \frac{u^n_{j+1} - 2u^n_j + u^n_{j-1}}{\Delta x^2}.$

If you do a von Neumann analysis, the stability criteria is $\Delta t \leq \Delta x$. This should give you no issues near the boundary.

I think the point Rezzolla is making in his notes is that you have to be very careful how you discretize the first order system to ensure stability. Generally, schemes for first order equations have to be upwind (or monotone) for stability and convergence. There is no notion (to my knowledge) of upwind scheme for systems of first order PDE, hence the difficulty.

If you want to see why your scheme is unstable, do a von Neumann analysis. The way to do it here is to plug $v^n_{j+1}$ and $v^n_{j-1}$ into the scheme for $u^{n+1}_j$. After a lot of simplification you will get a scheme in terms of only $u$. Then do the usual von Neumann analysis: Look for a separable solution $u^n_j =\lambda_k^n e^{ikj\Delta x}$ and find $\lambda_k$. It should turn out that no mater the grid resolution, at least one eiegnvalue $\lambda_k$ has modulus larger than 1, so the scheme is unstable.

  • $\begingroup$ Hi Jeff, thanks for your prompt reply. I have solved it with centered differences for both time and space coordinates and it definitely works. The problem is that I am implementing open boundary conditions and it seems I require to solve the two first-oder equations I just mentioned in the post $\endgroup$ – caverac Nov 2 '16 at 21:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.