Introduction

For some f[u], nonlinear wave equations of the form
2
∂
u
∂
2
t
=
2
∂
u
∂
2
x
+f[u]
have solutions that look like staircases, e.g. feature several distinct levels separated by steep inclines. The formation of these "staircases" depends on initial conditions and the form of f[u]. The border between two levels propagates with the speed approximately equal to 1, and gets more and more sharp with time. The region of the staircase solution that is inside the borders the deviation of u from the nearest equilibrium is small, so there the equation can be approximated as linear.
General information about nonlinear wave equations
This chapter discusses some general properties of nonlinear wave equations.

Brief overview

Nonlinear wave equation are equations of the form
2
∂
u
∂
2
t
=
2
∂
u
∂
2
x
+f[u]
(
1
)
Where f[u] is a real-valued function and x goes from -∞ to +∞
The equation (1) is a nonlinear partial differential equation. It is described by the Lagrangian and conserves the energy function, which puts restrictions on its' solutions, and in particular means that u,
u
t
and
u
x
have to remain finite. For any given value of energy for a certain type of f[u], the value of u is bounded to lie in a certain region. The equation can exhibit stationary, oscillatory, unbounded, or singular behavior, depending on initial conditions and the form of f[u]. The equation allows waves to travel along the x-axis, and for the vast majority of f[u] the speed of the waves depends on the waves themselves.

Numerical solutions

Let's have a look at several numerical solutions to (1) to get the feel of it. I am going to use the NDSolve function and plot solutions using Plot3D.
Firstly, let us specify the initial conditions. The choice is arbitrary, however I'm going to start with a Gaussian. This function is defined for all x, symmetric with respect to x = 0, and decays quickly while being infinitely differentiable for all values of x.
In[]:=
initConds = ​​ u[0, x]  Exp[-x^2],​​ Derivative[1, 0][u][0, x]  0​​;
Generally, u[t, x] denotes the value of u at time t and position x, and Derivative[n, m][u][t, x] is the same as
2
∂
u
n
∂
t
m
∂
x
(t,x)
. u[0, x] and Derivative[1, 0][u][0, x] represent the initial values of u and
u
t
, respectively.
Here and in the rest of the article
u
t
is a short notation for
∂u
∂t
. Similarly,
u
x
≡
∂u
∂t
,
u
tt
≡
2
∂
u
∂
2
t
, and
u
tx
≡
2
∂
u
∂t∂x
.
Below is the plot of initial values of u. As you can see, it decays very quickly, and even for x = 3 it’s already almost zero.
Plot[Exp[-x^2],{x,-3,3}]
Now, let us define the boundaries. Ideally, (1) should be solved within an infinite interval of x. However, for practical reasons, let us instead solve the equation in the finite interval x ∈ [ -R, R ] and t ∈ [ 0, T ], where R and T some positive numbers. I will take
In[]:=
T = 20;​​R = 32;
Finally, we have the boundary conditions. In the original problem, the domain of x is infinite, and we don't have to worry about those. In the numerical approximation, however, boundary conditions are compulsory. I will choose them to be
u
x
(t,-R)=
u
x
(t,R)=0
, for a number of reasons. Firstly, they are simple and familiar. They will not introduce any unexpected motion or disturbance. Secondly, they serve as a decent approximation to the missing semi-infinite string. Thirdly, as we will discover soon, these boundary conditions will prevent energy flux in or our of the studied interval.
In[]:=
boundaryConds=​​ Derivative[0, 1][u][t, -R]0,​​ Derivative[0, 1][u][t, R]0​​;
Lastly, we have our equation. In general, it is written as
In[]:=
eqnTemplate =( Derivative[2, 0][u][t, x]Derivative[0, 2][u][t, x] + f[u[t, x]]);
When we will need to study the case of a specific f[u], we can just use replacement like that:
eqnTemplate/.fSin
Or like that:
eqnTemplate/.fFunction[u,u-u^3]
The Function[u, func] thing means that we are defining a function called func with the argument u. Note: Function[u, u - u^3] gives the same result as Function[x, x - x^3].
This is how we use NDSolve to obtain a numerical solution for (1):
soln1 = NDSolve[​​ Join[initConds, boundaryConds, {eqnTemplate /. fFunction[u, -u]}],​​ u, {t, 0, T}, {x, -R, R}​​ ] (* the curly brackets around f[u[t, x... are here because of the rules of Wolfram language; ​​ they don't change the meaning *)
Out[]=
uInterpolatingFunction
Domain: {{0.,20.},{-32.,32.}}
Output: scalar

Now, we can plot our solution, just like that:
Plot3D[​​ u[t,x]/.soln1,​​ {t,0,T},{x,-R,R},​​ PlotRangeAll,​​ PlotPoints50​​ ]
For future use, I've defined a function that automatically uses NDSolve to give numerical solutions, removing the need to specify a ton of arguments every time we want to solve our beloved equation.
In[]:=
solveThePDE[replacementRuleForF_Rule] := NDSolve[​​ Join[initConds, boundaryConds, {eqnTemplate /. replacementRuleForF}],​​ u, {t, 0, T}, {x, -R, R}​​ ];
Here is how it works:
solveThePDE[fCos]
Also, here is the plotter function:
In[]:=
Options[myPlot] = {Rasterize  True};​​myPlot[solutionToPlot_, OptionsPattern[]] := Module[{plt},​​ plt = Plot3D[​​ u[t, x] /. solutionToPlot,​​ {t, 0, T}, {x, -R, R},​​ PlotRange  All,​​ PlotPoints  75​​ ];​​ ​​ If[OptionValue["Rasterize"],​​ Rasterize[plt],​​ plt​​ ]​​ ];
This "Rasterize" option controls whether the plot will be converted to an image or not. Converting dynamic plots to static images allows to reduce the load on the computer, however sometimes we might want to be able to spin and play around with the original plot.
Lastly, in the rest of this article I'll be working mainly not with f[u] but with V[u], as it gives more insight into the problem. NDSolve, however, only accepts actual differential equations and thus needs f[u]. Thus, I've created a function to calculate f[u] from V[u].
VtoF[Vrule_] := With[{expr = -D[(V /. Vrule)[u], u]},ToExpression["Function[u, " <> ToString[InputForm[expr]] <> "]"]]
It looks complicated, but basically what it does is take the derivative and multiplies the result by -1. Here is an example.
VtoF[VFunction[u,u^3]]
The derivative of
3
u
with respect to u is
3
2
u
, and f = -V' =
-3
2
u
.

Lagrangian

As Stephen Wolfram has shown, equation (1) can be described by the Lagrangian with the density
ℒ=
2
u
t
-
2
u
x
2
-V[u]
(
2
)
where
V[u]=-
u
∫
0
f[u']du'
(
3
)

Symmetry, Noether’s theorem, and conservation laws

Time translation symmetry

It can be seen that (2) is translationally symmetric with respect to t and x. The invariant corresponding to time symmetry, energy, is, as Stephen Wolfram has shown,
E=
∞
∫
-∞
ρdx,
dE
dt
=0,​​ρ=
2
u
t
+
2
u
x
2
+V[u]

Space translation symmetry

However, because (2) is also translationally symmetric with respect to x, there is one more energy-like invariant. It can be found using the same formula for Hamiltonian used to determine the energy function, except that instead of t there should be x.
G=
∞
∫
0
gdt,
dG
dx
=-
u
t
u
x
∞

0
,​​g=
-
2
u
t
-
2
u
x
2
+V[u]
The physical meaning of G is not entirely clear to me. In addition, since it is not constant throughout the x-axis, it's hard to find a practical application for it.

Hyperbolic rotation symmetry

It can be shown that the equation (1), as well as , are written in exactly the same way!

Energy flow

Let us derive an expression for the flow of energy, that I will denote J. In order to do that, let us use the law of conservation of energy.
Since the problem is one-dimensional, del operator is identical to ∂/∂x. Then, using the definition of ρ , we obtain
Which can be simplified to
Using the equations (1) and (3), we arrive at the following.
Where C is an arbitrary constant. From physical reasoning, the flow of energy must be zero if u and all its’ derivatives are zero. Then, the expression for the flow of energy is the following.

Stability of space-independent form of (1)

In this section, I'm studying the space-independent form of (1) to better understand the original equation. I've identified four types of behavior, one stable and two unstable.
Understanding the stability of (1) is an important stage in understanding the equation as a whole. But first, let us study the stability of space-independent form of (1), when u[t, x + Δx] = u[t, x] for any Δx. As we will later see, the general form of (1) exhibits similar types of behavior to its' space-independent form.
This is the equation of motion of a unit mass experiencing a force f. It can exhibit four types of behavior: stationary, oscillatory, unbounded, and singular. Let us look at each of those in more detail.

Stationary

Oscillatory

The equation (6) can be re-written as
Integrating this, we obtain
Where C is an arbitrary constant.
​
From this equation, it can be seen that the mass can only be moving if V[u] < C. If the mass is surrounded from both sides by potential barriers higher than C, it will never overcome any of them and instead will be oscillating somewhere in-between.
then u will forever remain bounded within that interval.
The plot of u[t] then looks something like this.

Unbounded

Singular

If the point is not bounded by potential barriers and goes off, let’s say, to infinity, and
then the point will be moving infinitely fast. An example of a singular solution is plotted below.
Actually, there are two types of a singular solution: those that actually have a singularity and those that just rise quickly. To see why is that, let us use the following equation derived earlier.
From here,
And
Let us find the time it takes for the solution to go to infinity. If the solution actually has a singularity, u will reach infinity in finite time. Else, the solution is simply a quickly increasing function.
However, for our purposes quickly rising solutions and singular solutions are not too different from one another, so they both fell under one classification slot.

Stability of general form of (1)

When the derivative of u with respect to x is not zero, the equation (1) can no longer be reduced to (6). However, the stability classification system developed for (6) is still applicable. As earlier, there are four types of behavior, listed below.

Stationary

Here, everything is clear without words.

Oscillatory

In the case of (1), the criterion for oscillatory behavior is not as strict as that for (6). However, if the plot of V[u] looks like a pit, the behavior of (1) is likely to be oscillatory.

Unbounded

Just like (6), the equation (1) can allow u to escape to infinity with a finite speed.
If the plot of V[u] looks somewhat like this, then the solution increases with a finite speed.
You can see that the solution growth with a steady rate bounded from above by the law of conservation of energy.

Singular

Some functions V[u] just yield singularities. If the value of u manages to escape to infinity and the limit of V[u] is negative infinity, u starts increasing infinitely fast.
As you can see, the flux of energy is positive to the right from the center and negative to the left. What this means is that to the right of the center the energy is flowing to the right, and to the left of the center it's flowing to the left, which basically implies that energy is being drained from the center.
This figure shows the direction of flux of energy around the center.

Upper and lower bounds for u

The law of conservation of energy determines the maximal and minimal possible value u can have for a given energy E. In this subsection, I will prove the latter statement and find the bounds for u. For the sake of shortness, I will only find the upper bound and not the lower bound, as the process of their derivation is practically identical.

Why are there any bounds for u?

Due to the law of conservation of energy, the value of u is restricted to some finite interval. To see why that is, let us examine the expression for the energy of the system.

Finding the bounds

Now, the question is how to find an optimal function u[t, x]. The word an here is important, as due to translation symmetry we can get as many optimal functions u[t, x] from one of them by shifting it either by either the t or x axis. To avoid ambiguity, let us restrict ourselves to such u[t, x] that have their global maximum at t = 0, x = 0.
Energy is the limiting factor in the process of increasing u. Thus, the optimal solution will be such that uses the minimal amount of energy.
So, finally, we get an equation describing the optimal solution.
Or
Where C is the constant of integration. Substituting the definition of y, we obtain
One possible solution to the problem is to make C nonzero and forcefully equate u to zero at some two points and further away from the origin. This would satisfy the initial conditions and allow us to reach any values of u. Below you can see an illustration explaining how this would work. The solution is smooth everywhere except the points x = 0, x = π / 2 and x = -π / 2.
In this case, the value of C should be chosen so that u would reach a certain value.
Anyway, the shape of this ideal solution is described by the equation

Speed of signal propagation

For the vast majority of f[u], the speed at which waves travel depends on the shape of the waves themselves. In this section the latter is demonstrated on the example of linear f[u]. To prove the point, I perform mathematical manipulations and create some numerical solutions.
For the vast majority of f[u], the speed at which waves travel depends on the shape of the waves themselves. To see this, we can study the case of linear f[u], when
And (1) takes the form
Since this is a linear equation, its' solution will be a sum of exponential functions of the form u[t, x] = C Exp[i(ω t - k x)], where i is the imaginary unit. Solving the characteristic equation, we arrive at the following nonlinear relation between ω and k.
The explanation for the above results is the following. From the previous dispersion relation, we can find the group speed, which equals
Below is the code I've used to generate the two plots. The function solveThePDEwithCustomInitConds does the same thing as solveThePDE, except that it also allows to specify the initial conditions.
myDensityPlot is similar to myPlot in a sense that it, too, plots u[t, x]. The difference between them, however, is that while myPlot plots u by drawing a surface in a 3D-space, myDensityPlot does it by drawing a 2D-canvas, where the color of a pixel is determined by the value of u at that point. Blue corresponds to low values of u, while orange denotes high values.
The code below basically applies the aforementioned functions to generate two density plots, which are then merged together.
Staircase solutions
This chapter defines staircase solutions and explores their formation and properties.

What is it?

Staircase solutions are solutions to the equation (1) that feature several distinct "levels" separated by steep inclines, e.g. look like staircases.
An important feature of staircase solutions is that while the levels get wider with time, their height remains unchanged.

Numerical solutions

In this short section, I will introduce some tools for numerically solving (1) that will be used later in the chapter.
The process of computing a staircase solution is almost the same as before, except that more accuracy is required.
To achieve it, I will use the function solveThePDEwithCustomStep , which is the same as solveThePDE, except that it allows to specify MaxStepSize. Here is how it works.
And this is what the plot looks like.
You can see that the value of u suddenly jumps from one level to another. Below is presented the slice of the above plot at t = 20.
Sharp transitions are what makes staircase solutions hard to compute numerically. To achieve sufficient accuracy it's necessary to use small time steps and fine grid.

Why this happens

The process of formation of staircase solutions

There are are two parallel processes that need to happen in order for a staircase to form. These are:
◼
  • Formation of borders
    • ◼
  • Formation of levels
    • Let us investigate each of them in more detail.
    As an example to support my claims, I will be using a numerical solutions to (1) computed with f[u] = Cos[4 u] / 4. Below, you can see the code I used to compute it.
    Here MaxStepSize was set to 0.01 to get sufficient accuracy.

    1. Formation of borders

    One of the key features of staircase solutions is the presence of sharp borders between adjacent levels. Let us see in detail how they form.
    The plot is the same as it would have been in the case of an ordinary wave equation, except for the f[u] term. If f[u] were zero, then the initial bump would simply split into parts that would propagate away from the origin without changing their shape. However, the force term drives them to deform and become increasingly steep as time goes by.
    However, since there are no floors in this solution, it is not considered a staircase solution.

    2. Formation of levels

    A process that goes in parallel with the formation of sharp borders is the formation of levels. Levels are regions of the solution where the value of u oscillates near some equilibrium value that doesn't change with time.

    Analogy with gas combustion

    The abovementioned process is in some sense similar to the process of subsonic gas combustion. Let us explore this analogy in more detail.

    Chain reaction

    Relative peace in the middle

    Behind the combustion wave, things are relatively calm - combustion products simply wander around at an almost fixed temperature, and nothing of much interest happens. The temperature and pressure are inhomogeneous, however not nearly as much as at the combustion edge.
    ​
    When it comes to PDE, things are very similar. At the boundary the values of derivatives of u are huge and the behavior of u[t, x] is highly complex. In the middle, however, the values of the derivatives are much smaller, and the shape of u is relatively simple.

    Diffusion counteracts contrast

    Parts of the gas that are currently burning have higher temperature than the cold regions, and higher temperature makes the gas burn, further increasing the temperature. If there was no heat transfer, the temperature inside the gas would abruptly jump from cold to very hot. However, since heat transfer counteracts extreme gradients of temperature, the transition from cold to hot remains relatively smooth.
    ​
    The same thing happens with the PDE. If not the second space derivative term, the value of the u-function would simply stick near one of the local minimums of potential energy. But, because of the spatial derivative term, the border is more or less smooth.

    What functions f[u] can lead to this

    There are two main criteria f[u] must satisfy to give rise to a staircase solution.

    Criterion 1: Roughness

    To give rise to a staircase solution if the initial conditions are Gaussian, a function f[u] must have a sharp transition between 0 and 1. In other words, the perfect plot of f[u] looks somewhat like this:
    And here is what the solution it yields look like:
    Let me explain why this is.
    If the function f[u] transitions too smoothly, a staircase does not form.

    Criterion 2: Energy

    Let's see why that's true. As we saw earlier, (1) obeys the law of conservation of energy.
    So,
    Just look at this plot! It looks like the first level is desperately trying to form, but keeps failing.
    It satisfies the roughness criterion. Just look at this plot of f[u]. f goes from negative to positive quite quickly.
    The secret ingredient is the energy criterion!
    This time, the plot of f is practically the same:
    However, the staircase finally forms!

    Studying the border in more detail

    In this section, I will be focusing on processes happening directly at the border between two levels. First, as a preparation, I will compute the position of one of the borders as a function of time. Then, I will study its' propagation and evolution with time.
    I will demonstrate my claims on the example of a particular numerical solution to (1), however the main points are true for all staircase solutions. The numerical solution I will be using as an example is soln3; it corresponds to the case when f[u] = Sin[4 u] / 4. Below you can see the 3D-plot and the density plot of this solution.

    Quantifying the position of the border

    In this subsection, I will quantify the position of one of the borders in soln3.
    Here is what the first four elements of the output look like:
    Here is the graph of the function.
    As you can see, the plot is practically linear throughout most of its' length. However, close to t = 0 it's visibly concave.

    Speed of propagation of the border

    The speed at which the border between two level propagates is almost constant and equal to 1. To demonstrate this, I will quantify the motion of one particular border and then use linear approximation to describe it. After that, I will provide mathematical explanation to the observed regularities.
    You can see that the two lines match almost exactly, and for the majority of the length of the plot the orange line matches the blue one so closely that it overshadows the latter. However, near the origin they differ quite significantly. In the next subsubsection we will see why things happen as they do and what consequences this brings.

    Mathematical explanation for the observations

    In this subsubsection, I will give explanation to what was shown empirically in the previous two, and also analyze the implications.
    From here, we can deduce that the speed at which the border between two levels propagates is approximately equal to 1. And as the border gets sharper, it's speed will approach 1 even closer. Conversely, if the border is not steep enough, the latter approximation becomes inaccurate and the speed of the border starts to differ significantly from 1.

    Evolution of the border with time

    The border between two levels typically gets sharper with time. I will make a change of variables in (1) and show that the resulting equation most often leads to sharpening of the borders. Then I will make some computations to support my claims.

    Changing variables

    In this subsubsection, I will be introducing delayed position, that makes it easier to study the change in the shape of the border as it travels.
    Now that we know that the border moves with a constant speed, let us study its’ motion with the help of delayed position, defined as
    Where c is the speed with which the border travels. As we saw in the previous section, the value of c is not precisely defined, however as the border gets steeper it approaches 1. We will assume that c changes with time slowly and can be approximated as constant.
    ​
    By the definition of z, the border has a certain value of z associated with it that remains almost unchanged as the border propagates. On the density plot below, you can see the numerical solution to (1) with f[u] = Sin[4 u] / 4, with a mesh over it. Horizontal lines show the lines of constant t, while the diagonal lines are the lines of constant z. (taking c to be 1.035)
    So, making the change of variables from t-x to t-x, we obtain
    Let us define a new function w, that is the same as u, except that it’s defined in terms of t and z instead of t and x.
    In terms of the new variables the equation (1) can be re-written as follows.
    Or
    Since, for steep inclines, c ≈ 1,
    Going back to t and x,
    This is a pretty neat relation, however we can go further than that. By the fundamental theorem of calculus,
    So we can write

    Assessing accuracy

    In this subsubsection, I will be evaluating the accuracy of the approximations derived above by plotting the right hand side and the left hand side of (8) and comparing them.
    Let us evaluate the accuracy of (8) and (9). For that, I've defined two interpolating functions, one for the left hand side and one for the right hand side of (8). Since to get from (8) to (9) I haven't used any additional approximation, the error of (8) equals that of (9).
    As you can see, the two curves match quite closely for large t. Thus, we conclude that, for large values of t, (8) and (9) are good approximations to (1).

    Implications

    The equation (9) has significant implications, as it tells the rate of change of the steepness of the border. In this subsubsection, I will be analyzing what this means for the rest of the article.
    If after infinite time the shape of the border approaches a certain limit, or equilibrium, then the equation (9) can tell us what conditions it must satisfy. The equation is written as follows.

    Studying the inside in mode detail

    The equation (1) inside the region bounded by the borders can be approximated as being linear.
    I will demonstrate my claims on the example of a particular numerical solution to (1), however the main points are true for all staircase solutions. The numerical solution I will be using as an example is soln3; it corresponds to the case when f[u] = Sin[4 u] / 4. Below you can see the 3D-plot of this solution.

    Linearity

    On each of the levels, the value of u does not differ very much from the nearest level. Then, we can expand f[u] into a Taylor series.
    One of the solutions to this equation is an exponential function
    Where ω and k are related as follows.
    You can see that the inner part of the solution to the equation (1) with f[u] = Sin[4 u] / 4 is very similar to that with f[u] = -u, as at the first level u = Pi / 4, and f'[Pi / 4] = -1.
    Even though that the solution on the left is visibly more rugged at the edges, it is clear that the look of the waves inside is almost identical.
    Bibliography
    ◼
  • Stephen Wolfram, A New Kind of Science
    • ◼
  • Stephen Wolfram, notes from A New Kind of Science