Message Boards

GROUPS:

Staff Picks Mathematics Calculus Equation Solving Wolfram Language Computer-Based Maths Wolfram High School Summer Camp

Rauan Kaldybaev

[WSC20] Exploring staircase solutions to nonlinear wave equations

Rauan Kaldybaev

Posted 6 months ago

1438 Views

2 Replies

5 Total Likes

Follow this post

enter image description here

Introduction

For some f[u], nonlinear wave equations of the form

∂

+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

∂

+f[u]

(

)

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,

and

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.

initConds =  u[0, x]  Exp[-x^2], Derivative[1, 0][u][0, x]  0;

In[]:=

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

∂

(t,x)

. u[0, x] and Derivative[1, 0][u][0, x] represent the initial values of u and

, respectively.

Here and in the rest of the article

is a short notation for

∂u

∂t

. Similarly,

≡

∂u

∂t

≡

∂

, and

≡

∂

∂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

T = 20;R = 32;

In[]:=

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

(t,-R)=

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

boundaryConds= Derivative[0, 1][u][t, -R]0, Derivative[0, 1][u][t, R]0;

In[]:=

Lastly, we have our equation. In general, it is written as

eqnTemplate =( Derivative[2, 0][u][t, x]Derivative[0, 2][u][t, x] + f[u[t, x]]);

In[]:=

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

uInterpolatingFunction

Domain: {{0.,20.},{-32.,32.}}

Output: scalar



Out[]=

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.

solveThePDE[replacementRuleForF_Rule] := NDSolve[ Join[initConds, boundaryConds, {eqnTemplate /. replacementRuleForF}], u, {t, 0, T}, {x, -R, R} ];

In[]:=

Here is how it works:

solveThePDE[fCos]

Also, here is the plotter function:

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 ] ];

In[]:=

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

with respect to u is

, and f = -V' =

-3

Lagrangian

As Stephen Wolfram has shown, equation (1) can be described by the Lagrangian with the density

ℒ=

-V[u]

(

)

where

V[u]=-

∫

f[u']du'

(

)

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,

∞

∫

-∞

ρdx,

=0,ρ=

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

∞

∫

gdt,

∞



,g=

+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

We can write (1) in the form when all operators are taken to the left hand side. In this case, it becomes



∂

)u=f[u]

. From here, it is easy to see that the equation is symmetric relative to hyperbolic rotation in the t-x plane. If we make the change of variables

t'=tcosh[ε]+xsinh[ε]x'=xcosh[ε]+tsinh[ε]

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.

∂ρ

∂t

+∇·J=0

(

)

Since the problem is one-dimensional, del operator is identical to ∂/∂x. Then, using the definition of ρ , we obtain

∂

∂t



+V[u]+

∂J

∂x

Which can be simplified to

+V'[u]

∂J

∂x

Using the equations (1) and (3), we arrive at the following.

∂J

∂x

∂

∂x

(

)J=C-

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.

J=-

(

)

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.

If u[t, x + Δx] = u[t, x], then, using the definition of derivative

lim

h0

u[t,x+h]-u[t,h]

it's easy to show that

for all t and x. Also since u no longer depends on x it can be written as u = u[t]. So, the equation (1) then can be written as

=f[u]

(

)

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

This type of behavior is probably the most boring of all. If at some time

[

]=0f[u[

]]=0

[

]=0

then the value of u will always be equal to

The plot of such a solution looks just like a straight line.

Plot[1,{t,0,1}]

Out[]=

Oscillatory

The equation (6) can be re-written as

d

=-dV

Integrating this, we obtain

=2(C-V[u])

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.

More formally, if at some point in time there exist such

and

that

<u[

>u[

]V[

]>V[u[

]]+

[

]

,V[

]>V[u[

]]+

[

]

then u will forever remain bounded within that interval.

<u[t]<

The plot of u[t] then looks something like this.

oscillatorySolnTo6 = NDSolve[{u''[t]  -Sin[u[t] - 1], u[0]  0, u'[0]  1.7}, u, {t, 0, 16}];

In[]:=

Plot[u[t]/.oscillatorySolnTo6,{t,0,16}]

Out[]=

Because u and

are connected by the law of conservation of energy, the function u[t] must be periodic.

Unbounded

If either of

mentioned earlier do not exist, the mass can go off to either negative infinity or infinity, respectively. If

and

are both undefined, the mass can escape in any direction. If the mass goes off to infinity and

a≤

lim

u∞

V[u]≤b

where a and b are some numbers, then

will be bounded. In other words, the value of u increases indefinitely but does so at a restricted rate. The plot of an unbounded solution looks somewhat like this.

unboundedSolnTo6 = NDSolve[{u''[t]  -Sin[u[t] - 1], u[0]  0, u'[0]  1.78}, u, {t, 0, 16}];

In[]:=

Plot[u[t]/.unboundedSolnTo6,{t,0,16}]

Out[]=

Singular

If the point is not bounded by potential barriers and goes off, let’s say, to infinity, and

lim

u∞

V[u]=-∞

then the point will be moving infinitely fast. An example of a singular solution is plotted below.

singularSolnTo6 = NDSolve[{u''[t]  u[t], u[0]  0, u'[0]  1.78}, u, {t, 0, 16}];

In[]:=

Plot[u[t]/.singularSolnTo6,{t,0,16}]

Out[]=

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.

=2(C-V[u])

If u goes off to infinity, the value of

will always remain positive. Then,

=√(2(C-V[u]))

From here,

√(2(C-V[u]))

And

t[u]=

√2

∫

u[0]

du'

√(C-V[u'])

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.

∞

√2

∞

∫

u[0]

du'

√(C-V[u'])

As we can see, a solution has a singularity if the integral

∞

∫

u[0]

du'

√(C-V[u'])

converges.

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.

Rasterize[Plot3D[1,{t,0,1},{x,-1,1}]]

Out[]=

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.

V1 = Function[u, u^2];

In[]:=

Plot[V1[u],{u,-1,1}]

Out[]=

In the case of V[u] =

, for example, the solution looks like this.

myPlot[solveThePDE[fVtoF[VV1]]]

Out[]=

Unbounded

Just like (6), the equation (1) can allow u to escape to infinity with a finite speed.

V2=Function[u, -u^2/(1+u^2)];

In[]:=

Plot[V2[u],{u,-4,4}]

Out[]=

If the plot of V[u] looks somewhat like this, then the solution increases with a finite speed.

myPlot[solveThePDE[fVtoF[VFunction[u,-u^2/(1+u^2)]]]]

Out[]=

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.

V3 = Function[u, -u^2];

In[]:=

Plot[V3[u],{u,-3,3}]

Out[]=

singularSoln=solveThePDE[fVtoF[VV3]];

In[]:=

myPlot[singularSoln]

Out[]=

An important thing about singular behavior is that if the u-function develops a singularity in one particular point, u also explodes in the neighboring points. To see why that's the case, let us plot J =

at time 20.

JforSingularSoln[x_] := -First[Derivative[1, 0][u][20, x] * Derivative[0, 1][u][20, x] /. singularSoln];

In[]:=

Plot[JforSingularSoln[x],{x,-R,R},PlotLegends{"J"},PlotRangeAll,PlotStyleOrange]

Out[]=

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.

VectorPlot[{JforSingularSoln[x],0},{x,-20,20},{y,-1,1},VectorPointsTable[{x,0},{x,-7,7,1}],AspectRatio1/6,ImageSizeLarge]

Out[]=

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.

∞

∫

-∞

+V[u]dx

One obvious way to make u large is to make it large at one point while making it zero everywhere else. However, as it's easy to show, that would make the integral

∞

∫

-∞

infinite. Spreading u over a large interval wouldn't be a particularly good solution too, as that would make

∞

∫

-∞

V[u]dx

large. Then, we conclude that the optimal arrangement lies somewhere in-between those two. The figure below shows the qualitative plots of the three solutions: too sharp, too smooth, and optimal.

Plot[{If[Abs[x]<0.1,0.5,0],0.75Exp[-x^2/8],Exp[-Abs[x]]},{x,-3,3},PlotStyleThick,AxesOrigin{-2,0},PlotLegends{"Too sharp","Too wide","Optimal"},ExclusionsNone]

	Too sharp
	Too wide
	Optimal

Out[]=

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.

From the expression for energy, it is obvious that at t = 0,

must be zero everywhere: since

only enters the integral under the square sign, making

anything other than zero would increase the energy while not increasing the maximal value of u. Later in this section we will practically forget about the time dependency of u and, for convenience, call u[x] = u[0, x]. With u and

things are more complicated: we need to find a balance between making the solution too sharp and giving too much energy to the V[u] term. For that, we can use the Euler-Lagrange equation. Since

∞

∫

-∞

ρdx

, the equation will be written in terms of ρ.



∂ρ

∂

=

∂ρ

∂u

Since

ρ=

+V[u]

∂ρ

∂

and

∂ρ

∂u

=-f[u]

. Thus,

[

]=-f[u]

So, finally, we get an equation describing the optimal solution.

=-f[u]

It can be simplified if we introduce an auxiliary variable

y≡

. THe latter can then be re-written as follows.

=-f[u]

ydy=-f[u]du

But, since y dy =

d



and -f[u] du = dV,

d

=dV

=2(V+C)

Where C is the constant of integration. Substituting the definition of y, we obtain

=2(V[u]+C)

We can understand the meaning of C by setting u to 0. Since V[0] =

∫

f[u]du

= 0, C =

u=0

. Here

u=0

means the square of

at the point when u = 0. So, C basically measures how steep the plot of u is when it crosses the u = 0 surface.

By the initial conditions, the value of u and

at infinity equals zero. One way to satisfy these conditions is to set C to zero. This, however, is not always viable, as the shape of the solution would then be governed by the equation

=2V[u]

. This equation only allows the value of u to change with x if V[u] > 0. What this means is that the value of u cannot surpass points where V[u] is zero. But since we also want to compute an optimal solution for arbitrary peak value of u, we need to find a different approach.

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.

Plot[UnitStep[Pi/2-Abs[x]](Pi/8+Cos[x]-Abs[x]/4),{x,-3,3},PlotStyleThick]

Out[]=

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

=±√(2(V[u]+C))

The sign of

depends on which side you are at: on the left it's positive and on the right it's negative. It is clear that the solution will have a mirror symmetry around the origin, so let us focus on the right half of the solution where

is negative.

=-√(2(V[u]+C))

Since

=1

, we can find x as a function of u by flipping the latter expression and integrating.

x[u]=

∫

max

du'

-√(2(V[u']+C))

√2

max

∫

du'

√(V[u']+C)

After that, we can get u[x] as an inverse function of the latter integral. After that, energy as a function of maximal u can be found by integrating the energy density. Finally,

max

can be found as an inverse function of

max

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

f[u]=-cu

And (1) takes the form

∂



∂

-cu

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.

As it's known, nonlinear dispersion relations cause different types of waves to travel with different speeds.Below you can see two plots, both showing numerical solutions to the same PDE

∂



∂

-24u

on the same interval of values of t and x. The horizontal axis is x (from left to right), and the vertical axis is t (from bottom to top).The only difference between the two plots are the initial