# Finite difference scheme for “wave equation”, method of characteristics

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.

• I am a bit confused by your subscripts, but this looks to me to be some sort of finite-difference time-domain formulation . . . perhaps with a staggered mesh formulation (half indices?). – meawoppl Sep 7 '12 at 19:05
• @meawoppl: He just calls his variables $u,v$ instead of $x,t$ as commonly done. (In the usual $u,v$ formulation, they are also rotated by $45^\circ$ in the space-time plane against $x,t$, but that's a separate matter.) – Wolfgang Bangerth Sep 7 '12 at 23:21
• I have edited to clarify (Wolfgang Bangerth´s explanation is what I had in mind). – Willie Wong Sep 10 '12 at 8:48

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).

• "but Google will help you more" Actually that is one of the big problems. I am not exactly sure what to Google for (I suspect the numerical literature may use some different terms from the pure literature). If you can suggest some keywords that I should search for, I'd be grateful. ("Method of lines", for example, is pointing me to a veritable wealth of information [perhaps even a bit much for me to be able to filter through :-) ].) – Willie Wong Sep 10 '12 at 11:05
• @WillieWong - One reference for hyperbolic equations that we commonly cite is LeVeque's Finite Volume Methods for Hyperbolic Problems. I'm not sure if this is the right reference for you to start with, but it will at least provide you an introduction to the terms and techniques in the field. – Aron Ahmadia Sep 10 '12 at 11:35
• Okay, I added some keywords and references. I hope they help. – David Ketcheson Sep 10 '12 at 19:44
• Much thanks for the references! That got me a good start. – Willie Wong Sep 11 '12 at 11:13

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.

• Wow, I can't believe I didn't realize that this was in the CFL paper... – David Ketcheson Sep 17 '12 at 6:43