Finite difference scheme for "wave equation", method of characteristics

Question

Consider the following problem $$ W_{uv} = F $$ where the forcing term can depend on $u,v$ (see Edit 1 below for the formulation), and $W$ and its first derivatives. This is a 1+1 dimensional wave equation. We have initial data prescribed at $\{u+v = 0\}$.

I am interested in the solution inside the domain of dependence of an interval $$\{ u+v = 0, u \in [- u_M,u_M]\}$$ and am considering the following finite difference scheme.

• The goal is to evolve $W_u$ by $W_u(u,v+1) - W_u(u,v) = F(u,v)$ and similarly $W_v(u+1,v) - W_v(u,v) = F(u,v)$. This scheme is integrable in the sense that $$ W(u,v) + W_u(u,v) + W_v(u+1,v) = W(u+1,v+1) = W(u,v) + W_v(u,v) + W_u(u,v+1)$$ so I can consistently compute $W$ from the initial data by integrating upwards; hence I only really need to look at the evolution equations for $W_v$ and $W_u$.
• For the initial data, we need the compatibility condition $W_u(u,v) - W_v(u+1,v-1) = W(u+1,v-1) - W(u,v)$. Which suggests that I can compute the initial data by using the forward (in $u$) finite difference of $W$ on the initial time with the values of given $W_t$ at half-integer points $(u+0.5,v-0.5)$.

Question:

1. Is this a well known scheme? In particular, where can I find analysis of this scheme?
2. Any thing obvious I should look out for?

Background: Pretend I know next to nothing (which is probably true, as I am a pure mathematician trying to learn a little bit of computation machinery).

---

Edit 1: Just to clarify (to address some comments): the equation in $x$ $t$ coordinates would be $$ W_{tt} - W_{xx} = F $$ and $u$ and $v$ are ¨null coordinates¨ given by (up to some renormalising factors of 2) $u = t+x$ and $v = t-x$. So the initial data at $\{u+v = 0\}$ is in fact at $\{t = 0\}$.

So instead of a mesh adapted to $(t,x)$ I consider a mesh adapted to $(u,v)$ which is ¨rotated 45 degrees¨. Compared to the $(t,x)$ where $t,x$ take integer values, one can think of the $u,v$ mesh as having additional points where both (but not just one of) $t$ and $x$ take half-integer values.

Answer 1

There is definitely literature on schemes like this. Two keywords are

• Modified method of characteristics
• Semi-Lagrangian schemes

After 20 minutes of googling: some possibly important papers are http://dx.doi.org/10.1137/0719063 and http://dx.doi.org/10.1137/0728024 (search forward from there). Those probably aren't the best references out there, but they should be a starting point to get you into the right literature.

I think of this as a rotated method of lines with dimensional splitting. Presumably you're very well aware of the equivalence of your equation and the usual form of the wave equation $$W_{tt}-W_{xx} = F$$ under the transformation $$u = t+x, \ \ \ \ v = t-x.$$ For me it is useful to think of your scheme in terms of this traditional form of the wave equation. What the scheme does is integrate first along one set of characteristics, then along the other. The integration is done using dimensional splitting and Euler's method, both of which are first order accurate.

Of course, since you're integrating along characteristics, your scheme would be exact in the case $F=0$. That is, the numerical errors in your scheme will be due only to numerical integration of $F$ (this may be obvious, but is perhaps useful to point out to those who are accustomed to more traditional numerical methods). Furthermore, your scheme is unconditionally stable for the case $F=0$. Nothing more can be said about its stability without knowing some properties of $F$. In general, the scheme will be stable only under some finite step size restriction (since Euler's method is explicit). If the Jacobian of $F$ has any purely imaginary eigenvalues, the scheme will be unstable.

The general discretization approach of reducing a PDE to a system of ODEs (as in your method) is known as the method of lines. As with any method of lines discretization, you could increase the order of accuracy by using a higher-order ODE solver and you could improve the stability by using an appropriate implicit ODE solver (with the attendant increase in computational cost per step).

Answer 2

Starting from where David Ketcheson left me in his answer, a little bit more search revealed some historical notes.

The scheme I outlined above was considered already back in 1900 by J. Massau, in Mémoire sur l'intégration graphique des équations aux dérivées partielles. The work is republished in 1952 by G. Delporte, Mons.

The first (albeit brief) modern analysis of its convergence and such was given by Courant, Friedrichs, and Lewy's in their classic 1928 paper in Math. Ann.