My question is: how come that certain formulations of the wave equation can be time integrated more efficiently then others?

Le me expand a bit on that. Consider the wave equation:

$$ \frac{d^2 p(t,x)}{dt^2} -c^2 \frac{d p(t,x)}{dx^2} = 0 $$

with boundary conditions $ p(t,0) = \sin(\omega t) $ and initial conditions $ p(0,x) = 0 $ and $ p'(0,x) = 0$. One thing that is a bit of a concern is that these ics and bcs are not consistent.

This equation can be written as a system of first order equations. In fact, there are multiple versions the of systems of first order equations the original equation can be re-written in. Here is a system of first order equations that is equivalent to the above equation, let's call it A. We introduce an auxiliary variable $v$:

$$ \frac{d p(t,x)}{dt} - v(t,x) = 0 $$ $$ \frac{d v(t,x)}{dt} -c^2 \frac{d p(t,x)}{dx^2} = 0 $$

with boundary conditions $ p(t,0) = \sin(\omega t) $ and $ v(t,0) = \omega \cos(\omega t) $ and initial conditions $ p(0,x) = 0 $ and $ v(0,x) = 0$.

Next, is a second version the original equation can be re-written as, let's call that B:

$$ \frac{d p(t,x)}{dt} + c^2 \frac{d v(t,x)}{dx} = 0$$ $$ \frac{d v(t,x)}{dt} + \frac{d p(t,x)}{dx} = 0$$

with boundary conditions $ p(t,0) = \sin(\omega t) $ and $ v(t,0) = 1/c \sin(\omega t) $ and initial conditions $ p(0,x) = 0 $ and $ v(0,x) = 0$.

I am observing that B can be time integrated much more efficiently than A. Let me show that in Mathematica code.

We generate a single FEM mesh that is used in all cases. We set up some constants:

$HistoryLength = 0;
c = 343;
omega = 5025;
tend = 0.002;
mesh = ToElementMesh[Line[{{0}, {1}}], "MaxCellMeasure" -> 0.025]

This generates a 1D mesh with about 40 second order line elements and it used in all examples below.

Set up the original equations, bcs, and ics:

SecondOrderModel = D[p[t, x], {t, 2}] - c^2*D[p[t, x], {x, 2}] == 0;
SecondOrderBC = p[t, 0] == Sin[omega t];
SecondOrderIC = {p[0, x] == 0, Derivative[1, 0][p][0, x] == 0};
pde1 = {SecondOrderModel, SecondOrderBC, SecondOrderIC};

Set up A:

NDSolveFirstOrderModel = {
   D[p[t, x], {t, 1}] - v[t, x] == 0, 
   D[v[t, x], {t, 1}] - c^2*D[p[t, x], {x, 2}] == 0};
NDSolveFirstOrderIC = {p[0, x] == 0, v[0, x] == 0};
NDSolveFirstOrderBC = {p[t, 0] == Sin[omega t], 
   v[t, 0] == D[Sin[omega t], t]};
pde2 = {NDSolveFirstOrderModel, NDSolveFirstOrderIC, 

Set up B:

FirstOrderModel = {
   D[p[t, x], {t, 1}] + c^2*D[v[t, x], {x, 1}] == 0, 
   D[v[t, x], {t, 1}] + D[p[t, x], {x, 1}] == 0};
FirstOrderIC = {p[0, x] == 0, v[0, x] == 0};
FirstOrderBC = {p[t, 0] == Sin[omega t], 
   v[t, 0] == 1/c Sin[omega t]};
pde3 = {FirstOrderModel, FirstOrderIC, FirstOrderBC};

In order to do a fair timing comparison, we solve all equations once in order to warm up the jit compiler, load stuff etc.

NDSolveValue[pde1, p, {t, 0, tend}, Element[{x}, mesh]];
NDSolveValue[pde2, p, {t, 0, tend}, Element[{x}, mesh]];
NDSolveValue[pde3, p, {t, 0, tend}, Element[{x}, mesh]];

Solving the original second order equation:

 sol1 = NDSolveValue[pde1, p, {t, 0, tend}, Element[{x}, mesh]];]

{1.4, Null}

Solving A:

 sol2 = NDSolveValue[pde2, p, {t, 0, tend}, Element[{x}, mesh]];]
{1.4, Null}

In fact A is what Mathematica automatically generates from the original equation. So no surprise here.

Solving B:

 sol3 = NDSolveValue[pde3, p, {t, 0, tend}, Element[{x}, mesh]]]
{0.271091, Null}

So this is much faster to time integrate and solutions are all comparable. The reason that time integration is faster is that the automatic time step selection uses much lager steps. My question is why can B integrated with larger time steps than A. I am looking for references where I could read up on that.

  • 1
    $\begingroup$ Are you absolutely sure you're getting the same answer each time? Ignoring Mathematica and just solving the equations I get $p_1 = \frac{1}{2}\left(\sin\left(\omega\left(t+\frac{x}{c}\right)\right)+\sin\left(\omega\left(t-\frac{x}{c}\right)\right)+A\cos\left(\omega\left(t+\frac{x}{c}\right)-A\cos\left(\omega\left(t-\frac{x}{c}\right)\right)\right)\right)$ (you appear to be missing a BC for this one) $p_2 = \frac{1}{2} (\sin\left(\omega\left(t+\frac{x}{c}\right)\right)+\sin\left(\omega\left(t-\frac{x}{c}\right)\right))$ and $p_3 = \sin\left(\omega\left(t-\frac{x}{c}\right)\right) $ $\endgroup$
    – origimbo
    Apr 29, 2019 at 16:19
  • $\begingroup$ @origimbo, thank you for your comment. I'll check the solutions. You mention that p1 or p2 is missing a boundary condition, why do you think one is missing. One issue is that the initial conditions and the boundary conditions are not consistent for the first two examples and I need to find a way to sensibly fix that. $\endgroup$
    – user21
    Apr 30, 2019 at 8:19

1 Answer 1


It's not really meaningful to talk about integrating the equation in form A or B, since one way to integrate A is to first transform to B and then discretize. You can only really compare the actual discretizations, and I don't know what Mathematica is using.

That said, we can still make an educated guess. Since A is first order in time but second order in space, Mathematica probably discretizes it in a way that requires $\Delta t \propto \Delta x^2$. Meanwhile B is first order in both, so Mathematica uses a discretization that requires only $\Delta t \propto \Delta x$.

In short, B is a much more useful way to write the wave equation than A, both for purposes of analysis and computation. Mathematica is not really a great tool for solving PDEs.

  • $\begingroup$ Thank you for your time to answer. I appreciate that. I have updated the post to indicate that in all cases the same second order FEM mesh is used. Since the same mesh is used I do not think that your guess about the spatial complexity is correct. What interests me the most is some reference as to why B is more useful than A in terms of computation. Part of the issue I see is that the initial condition and boundary condition are not consistent. $\endgroup$
    – user21
    Apr 30, 2019 at 8:16
  • $\begingroup$ I disagree with your statement about Mathematica. As does the fact that the differential-equations tag is the third highest over at Methematica Stackexchange. Perhaps you should have a go at it once more, since the last time you used it... $\endgroup$
    – user21
    Apr 30, 2019 at 8:17
  • $\begingroup$ @user21 I did not say that they use different spatial meshes. Actually I imagined that they would both use the same mesh. If I have time later, I will give a more complete explanation, but in case I don't you can get the ideas by reading about discretizations of the heat equation and the wave equation in an introductory numerical analysis text. $\endgroup$ Apr 30, 2019 at 10:29
  • 1
    $\begingroup$ I don't think a debate about Mathematica is worthwhile here. $\endgroup$ Apr 30, 2019 at 10:30
  • 1
    $\begingroup$ A qualitative property like stability can only depend on some dimensionless number. The dimensionless number that can be extracted from A involves $(\Delta t)/(\Delta x^2)$ The dimensionless number extracted from B involves $(\Delta t)/(\Delta x)$ This justifies David Ketcheson's analogy. $\endgroup$
    – Philip Roe
    May 2, 2019 at 20:15

Your Answer

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

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