This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

Collect various “interesting” 1+1 and 1+2 PDEs and PDE systems from the literature (or self-constructed ones) and visualize the qualitative shape of their solutions.

We start with a simple system of two coupled equations depending from the two independent variables t and g [✶]. Derivative terms appear only for the variable t. This means, this system is basically a parametrized coupled system of ordinary differential equations. The solution shows a spiral-like behavior.

Module[{d = 0.4,  ppL = 5000, ndsol, x, y, t, g},
ndsol = (* pseudo PDE; no derivatives with respect to g *)
NDSolve[{D[x[t, g], t] == -1/2 Sinh[2 Pi y[t, g]]/
                               (Cosh[2 Pi y[t, g]] - Cos[2 Pi x[t, g]] + d^2),
         D[y[t, g], t] == +1/2 Sin [2 Pi x[t, g]]/
                               (Cosh[2 Pi y[t, g]] - Cos[2 Pi x[t, g]] + d^2),
         x[0, g] == -1/2 + g, y[0, g] == 0}, 
        {x, y}, {t, -1, 1}, {g, 0, 1},
         Method -> {"MethodOfLines", "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> 4, 
                    "MaxPoints" -> {ppL}, "MinPoints" -> {ppL}}}];
(* 3D plot of {y[t, g], y[t, g], t} *)
ParametricPlot3D[Evaluate[{x[t, g], y[t, g], t} /. ndsol[[1]]],
                 {g, 0, 1}, {t, -1/2, 1}, PlotRange -> All,
                 PlotPoints -> {240, 60}, Axes -> False]]

The following 1+2 PDE system generates interestingly structured labyrinth-like solutions [✶]. It is used to model brain growth. The model has four parameters γ, β, ν, and δ.

eqs1[{γ_, β_, ν_, δ_}] := 
 {D[ψ[t, x, y], t] == γ (η[t, x, y] - ψ[t, x, y]^3/3 + ψ[t, x, y]) + 
                      D[ψ[t, x, y], x, x] + D[ψ[t, x, y], y, y],
  D[η[t, x, y], t] == -1/γ (ψ[t, x, y] + ν + β η[t, x, y]) + 
                      δ (D[η[t, x, y], x, x] + D[η[t, x, y], y, y])}

The function ric implements random periodic initial conditions.

ric[n_, m_, seed_][x_?NumberQ, y_?NumberQ] := (SeedRandom[seed]; 
Sum[Random[Real, {-1, 1}] Cos[Random[Integer, {1, m}] 2Pi x + 2 Pi Random[]]*
                          Cos[Random[Integer, {1, m}] 2Pi y + 2 Pi Random[]],
    {j, 0, n}])

The function solveAndMakeGraphics solves the 2 PDEs and generates 3D plots of the two components.

solveAndMakeGraphics1[{γ_, β_, ν_, δ_}, {{nψ_, mψ_, seedψ_}, {nη_, mη_, seedη_}},
                     L_, T_, ppL_, do_, ndsolveOpts___, graphicsArrayOpts___] := 
Module[{ndsol},
(* solve PDE system *)
ndsol = NDSolve[Flatten[{eqs1[{γ, β, ν, δ}], 
         ψ[0, x, y] == ric[nψ, mψ, seedψ][x/L, y/L], 
         η[0, x, y] == ric[nη, mη, seedη][x/L, y/L],
         ψ[t, -L, y] == ψ[t, +L, y], ψ[t, x, -L] == ψ[t, x, +L],
         η[t, -L, y] == η[t, +L, y], η[t, x, -L] == η[t, x, +L]}],
         {ψ, η}, {t, 0, T}, {x, -L, L}, {y, -L, L}, 
         Sequence @@ {ndsolveOpts}, MaxSteps -> 10^5,
         Method -> {"MethodOfLines", "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                    "MaxPoints" -> {ppL, ppL}, "MinPoints" -> {ppL, ppL}}}];
(* make graphic of ψ and η *)
Show[GraphicsArray[
Plot3D[Evaluate[ndsol[[1, #, 2]][T, x, y]], {x, -L, L}, {y, -L, L}, 
       PlotRange -> All, PlotPoints -> 2 ppL, Mesh -> False,
       DisplayFunction -> Identity, Axes -> False,
       Axes -> {False, False, True}]& /@ {1, 2}], 
       Sequence @@ {graphicsArrayOpts}, PlotLabel -> do]]

Here are two examples of solutions. Because of the periodic boundary conditions, we use the Pseudospectral option setting.

Table[solveAndMakeGraphics1[{4.9331, 2.9971, 0.014014, 1.62681}, 
                           {{20, 4, k}, {20, 4, k + 1}},
                            100, 6, 120, Pseudospectral,
                           {PrecisionGoal -> 8, AccuracyGoal -> 8}], {k, 0, 1}]

The second example is a 1+1 PDE that is used for describing surface evolution due to step bunching [✶], [✶], [✶].

eq2[γ_] := D[h[t, z], t] == -D[h[t, z], {z, 2}] + γ D[h[t, z], {z, 3}] -
           D[h[t, z], {z, 4}] + D[D[h[t, z], z]^2, {z, 2}]

The function solveAndMakeGraphics solves the differential equation with the parameter value γ and the initial conditions ic. We use periodic boundary conditions.

solveAndMakeGraphics2[γ_, ic_, T_, pp_, do_] := 
Module[{ndsol, T1},
(* solve PDE *)
ndsol = NDSolve[{eq2[γ], h[0, z] == ic[z], h[t, -Pi] == h[t, Pi]},
                 h[t, z], {t, 0, T}, {z, -Pi, Pi}, 
                 PrecisionGoal -> 3, AccuracyGoal -> 3, 
                 MaxSteps -> 10^5, Method -> StiffnessSwitching,
 Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid",
            "DifferenceOrder"-> do, "MinPoints" -> pp, "MaxPoints" -> pp}}];
(* visualize solution *)            
Show[Graphics3D[      
T1 = ndsol[[1, 1, 2, 0, 1, 1, 2]];  
Plot3D[Evaluate[{Abs[h[t, z]] /. ndsol[[1]],
                 Hue[(Arg[h[t, z]] + Pi)/(2Pi) /. ndsol[[1]]]}], 
                 {z, -Pi, Pi}, {t, 0, 0.9 T1}, 
                PlotPoints -> 2 pp, Mesh -> False, PlotRange -> All, 
                DisplayFunction -> Identity]] /.
        {h_Hue :> SurfaceColor[h, h, 2.4], p_Polygon :> {EdgeForm[], p}},
        DisplayFunction -> $DisplayFunction, Lighting -> True]]

For real values of γ, the solutions smooth the original oscillations. Here are two examples.

SeedRandom[111];
Show[GraphicsArray[
Block[{$DisplayFunction = Identity}, 
Table[solveAndMakeGraphics2[0.3, Function[z, Evaluate[
 (* random periodic real initial condition *)
 Sum[Random[Real, {-1, 1}] Cos[k z + 2 Pi Random[]], {k, 6}]]], 
     0.01, 100, "Pseudospectral"], {2}]]]]

For complex values of γ, the solutions frequently also smooth the original oscillations. But sometimes complicated patterns and singularities form. Here are two examples.

(* suppress error messages *) Off[NDSolve::eerr];
Show[GraphicsArray[
Block[{$DisplayFunction = Identity}, 
Function[seed, SeedRandom[seed]; 
solveAndMakeGraphics2[Exp[I 2Pi Random[]], Function[z, Evaluate[
 Sum[Random[Complex, {-1 - I, 1 + I}] Cos[k z + 2 Pi Random[]], 
     {k, 0, Random[Integer, {3, 6}]}]]], 0.025, 200, 
     "Pseudospectral"]] /@ {2165915158297496362, 95157975033230510106}]]]

For exactly integrable solutions with a variety of different solutions, see [✶], [✶].

We continue with a variant of the Swift–Hohenberg equation [✶]. It is a one-parameter 1+1 PDE.

eq3[η_] :=  D[Ψ[x, t], t] + η Ψ[x, t] + D[Ψ[x, t], {x, 2}] + 
            D[Ψ[x, t], {x, 4}] +  2 Ψ[x, t] D[Ψ[x, t], x] == 0

Again, the function solveAndMakeGraphics3 solves the PDE and visualizes the solution as a contour plot. The last argument specifies the function to be displayed and the contour values.

solveAndMakeGraphics3[η_, L_, T_, ppL_, do_, {F_, cls_}] := 
Module[{ndsol},
(* solve PDE *)
ndsol = NDSolve[{eq3[η], Ψ[x, 0] == Sum[Random[Real, {-1, 1}] Sin[k x/L Pi], {k, 16}],
                     Ψ[0, t] == Ψ[L, t]},
         {Ψ}, {x, 0, L}, {t, 0, T}, MaxSteps -> 10^5,
         PrecisionGoal -> 3, AccuracyGoal -> 3, 
         Method -> {"MethodOfLines", Method -> StiffnessSwitching,
                    "SpatialDiscretization" -> 
                     {"TensorProductGrid", "DifferenceOrder" -> 6, 
                      "MaxPoints" -> {ppL}, "MinPoints" -> {ppL}}}];
(* visualize solution *)            
ContourPlot[Evaluate[F[Ψ[x, t] /. ndsol[[1]]]], {x, 0, L}, {t, 0, T},
            PlotRange -> All, PlotPoints -> 300, ContourLines -> False,
            Contours -> cls, ContourStyle -> {Thickness[0.0001]},
            FrameTicks -> {True, True, False, False}, 
            ColorFunction -> (Hue[0.8 #]&)]];

Here are the solutions for two sets of parameters η and the domain size L. The left graphic shows the regions Ψ(x,t)≶0 and the right graphic shows |Ψ(x,t)|. The qualitative behaviors of the solutions are quite different.

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
 {SeedRandom[888];
  solveAndMakeGraphics3[1.400, 180,  40, 2400, 6, {Identity, {0}}],
  SeedRandom[666];
  solveAndMakeGraphics3[0.087, 110, 100, 2400, 6, {Abs, 50}]}]]]

Next, we use a 2D version of the Swift–Hohenberg equation [✶], [✶].

eq4[ϵ_, α_] =  D[ψ[x, y, t], t] == ϵ ψ[x, y, t] - (* biLaplacian *)
                 Nest[(# + D[#, x, x] + D[#, y, y])&, ψ[x, y, t], 2]
                     α ψ[x, y, t]^2 - ψ[x, y, t]^3;

Starting with random initial condition in the inner part of the x,y-domain, yields in most cases characteristic labyrinth-like structures.

ricIUS[k_, {x_?NumberQ, y_?NumberQ}] :=
If[Max[Abs[{x, y}]] > 1/2, 0, 
   SeedRandom[ToString[N[k (x + Sqrt[2] y)]]]; Random[Real, {-1, 1}]]

solveAndMakeGraphics4[{ϵ_, α_, kSeed_}, L_, T_, ppL_, do_] := 
Module[{ndsol},
(* solve PDE *)
ndsol = NDSolve[{eq4[ϵ, α], ψ[x, y, 0] == ricIUS[kSeed, {x/L, y/L}],
                 ψ[-L, y, t] == 0, ψ[L, y, t] == 0, 
                 ψ[x, -L, t] == 0, ψ[x, L, t] == 0,
                 Derivative[1, 0, 0][ψ][-L, y, t] == 0,
                 Derivative[1, 0, 0][ψ][+L, y, t] == 0,
                 Derivative[0, 1, 0][ψ][x, -L, t] == 0,
                 Derivative[0, 1, 0][ψ][x, +L, t] == 0},
         {ψ}, {x, -L, L}, {y, -L, L}, {t, T, T}, MaxSteps -> 10^5,
         PrecisionGoal -> 2, AccuracyGoal -> 2, 
         Method -> {"MethodOfLines", Method -> Automatic,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL, ppL}, "MinPoints" -> {ppL, ppL}}}];
(* visualize solution *)   
ContourPlot[Evaluate[ψ[x, y, T] /. ndsol[[1]]], {x, -L, L}, {y, -L, L}, 
       PlotRange -> All, PlotPoints -> 160, Contours -> 20, 
       ContourLines -> False, FrameTicks -> False, 
       ColorFunction -> (Hue[0.8 #]&)]]

Here are two examples, each starting from two different random initial conditions. With a higher difference order or/and more points, we could, at the expense of a longer calculation time, obtain a more reliable solution.

Show[GraphicsArray[#]]& /@ 
Block[{$DisplayFunction = Identity},
{Table[solveAndMakeGraphics4[{1.154, 0.482, k}, 50, 15, 80, 2], {k, 2}],
 Table[solveAndMakeGraphics4[{1.135, -1.4, 3k}, 50, 20, 80, 2], {k, 2}]}];

The next PDE is of type 2+1 and has fourth order spatial derivatives [✶].

eq5[ϵ_] =  ϵ^2 D[A[x, y, t], t] ==  ϵ^2 A[x, y, t] - 
               (A[x, y, t] + 2 (D[A[x, y, t], {x, 2}] + D[A[x, y, t], {y, 2}]) +
                D[A[x, y, t], {x, 4}] + D[A[x, y, t], {x, 2}, {y, 2}] + 
                D[A[x, y, t], {y, 4}]) +
           1/3 (D[#[[1]], x] + D[#[[2]], y])&[
                 {D[A[x, y, t], x], D[A[x, y, t], y]}*
                  (D[A[x, y, t], x]^2 + D[A[x, y, t], y]^2)];

We will start with random initial conditions. The function randomPeriodicICs generates random initial conditions with boundary values that are periodic in [-1,-1]×[-1,-1].

randomPeriodicICs[k_, {x_?NumberQ, y_?NumberQ}] :=
If[x == -1. || x == 1. || y == -1. || y == 1., 1., 
   SeedRandom[ToString[N[k (E x + Pi y)]]]; Random[Real, {-1, 1}]]

Next we solve the PDE for the value ε=-0.2833. While a difference order 2 does not suffice to deal with the spatial fourth order derivative, NDSolve is smart enough to use the smallest order that will work. While carrying out the following input, we get some NDSolve::eerri messages. Taking the noncontinuous nature of the initial conditions into account, these messages are expected.

ndsol = 
Module[{ppL = 80, L = 40, T = 4, ϵ = -0.2833},
NDSolve[{eq5[ϵ], A[x, y, 0] == randomPeriodicICs[1, {x/L, y/L}],
                 A[-L, y, t] == A[L, y, t], A[x, -L, t] == A[x, L, t]},
         {A}, {x, -L, L}, {y, -L, L}, {t, 0, T}, MaxSteps -> 10^5,
         PrecisionGoal -> 2, AccuracyGoal -> 2, 
         Method -> {"MethodOfLines", Method -> StiffnessSwitching,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> 2, 
                     "MaxPoints" -> {ppL, ppL}, "MinPoints" -> {ppL, ppL}}}]];

Here are a sequence of graphics that show the evolution of the system from the random initial values to an ordered square pattern [✶].

Module[{pp = 160, L = 40, T = 4},
Show[GraphicsArray[#]]& /@ Partition[
ContourPlot[Evaluate[A[x, y, #] /. ndsol[[1]]], {x, -L, L}, {y, -L, L}, 
       PlotRange -> All, PlotPoints -> pp, Contours -> 20, 
       ContourLines -> False, DisplayFunction -> Identity,
       FrameTicks -> False, ColorFunction -> (Hue[0.8 #]&)]& /@
       (* list of times for frames *)
       Prepend[T Table[2^-k, {k, 7, 0, -1}], 0], 3]]

Our next example is the Meinhardt equation set, a coupled system of three nonlinear reaction diffusion equations [✶], [✶] that are used to model animal coats [✶]. The system has twelve numerical parameters.

eqs6[{{Da_, ρa_, κa_, Ds_, ρs_, μs_, σs_, ρY_, μY_, σY_}, {κs_, κY_}}] = 
{D[a[x, y, t], t] == Da (D[a[x, y, t], x, x] + D[a[x, y, t], y, y]) +
                     ρa (a[x, y, t]^2 s[x, y, t]/(1 + κa a[x, y, t]^2) -
                     a[x, y, t]),
 D[s[x, y, t], t] == Ds (D[s[x, y, t], x, x] + D[s[x, y, t], y, y]) +
                     σs/(1 + κs Y[x, y, t]) -  μs s[x, y, t] - 
                     ρs a[x, y, t]^2 s[x, y, t]/(1 + κa a[x, y, t]^2),
 D[Y[x, y, t], t] == ρY Y[x, y, t]^2/(1 + κY Y[x, y, t]^2) - 
                     μY Y[x, y, t] + σY a[x, y, t]};

This time we use initial conditions that are constant for two of the dependent variables and mostly constant for the third. The third variable has a nonvanishing value only on a few random points.

(* some nonvanishing point values; boundary value is 0 *)
a0IC[{A_, o_, ppL_}, {x_?NumberQ, y_?NumberQ}] := 
If[Max[Abs[{x, y}]] == 1., 0, If[Random[] < o 1/(ppL - 2)^2, A, 0]]

The function solveAndMakeGraphics6 solves the Meinhardt system and generates a contour plot of the last variable.

solveAndMakeGraphics6[{T_, L_}, pdePars_, {A_, o_}, ppL_, do_] :=
Module[{ndsol},
(* solve pde system *)
ndsol = NDSolve[{eqs6[pdePars], a[x, y, 0] == a0IC[{A, o, ppL}, {x, y}/L], 
                 s[x, y, 0] == 0, Y[x, y, 0] == 0,
                 a[-L, y, t] == a[L,  y, t], a[x, -L, t] == a[x, L, t], 
                 s[-L, y, t] == s[L,  y, t], s[x, -L, t] == s[x, L, t], 
                 Y[-L, y, t] == Y[L,  y, t], Y[x, -L, t] == Y[x, L, t]},
         {a, s, Y}, {x, -L, L}, {y, -L, L}, {t, T, T}, MaxSteps -> 10^5,
         PrecisionGoal -> 2, AccuracyGoal -> 2, 
         Method -> {"MethodOfLines", Method -> BDF, "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> ppL {1, 1}, "MinPoints" -> ppL {1, 1}}}]; 
(* make contour plot of Y[x, y]^2 *)
ContourPlot[Evaluate[Abs[Y[x, y, T]]^2 /. ndsol[[1]]], {x, -L, L}, {y, -L, L}, 
       PlotRange -> All, ContourLines -> False, FrameTicks -> False, 
       PlotPoints -> 120, Contours -> 30, ColorFunction -> (Hue[0.8 #]&)]]

Here are three examples. The first two show a Voronoi diagram-like structure and the third shows a mixture between spots and stripes. We use 80 points in each spatial direction. While this suffices to obtain a qualitatively correct solution, it does not suffice to obtain a quantitatively correct one within the specified accuracy and precision goals. For the price of larger calculation time, one could calculate a better solution.

Show[GraphicsArray[
Block[{$DisplayFunction = Identity}, 
     (SeedRandom[111]; Print[#];
      solveAndMakeGraphics6[#1, #3, {#2, 20}, 80, 1])& @@@ 
 (* parameter data found by searching *)
 {{{648, 66}, 5.58964, 
   {{0.015803, 0.070456, 0.003254, 0.826182, 0.001030, 0.006167, 
     0.088085, 0.351054, 0.003161, 0.010544}, {2.07226, 3.111023}}},
  {{641, 57}, 9.0855, 
   {{0.065498, 0.247513, 0.002189, 0.391285, 0.070575, 0.005826, 
     0.458718, 0.189077, 0.001017, 0.008231}, {2.38141, 0.235835}}},
  {{437, 20}, 9.85503, 
   {{0.001890, 0.039474, 0.002056, 0.887802, 0.003273, 0.006249, 
     0.537041, 0.006767, 0.115995, 0.140463}, {3.04301, 1.058715}}}}]]]

A rich source of complicated patterns are the cubic, quintic, … complex Ginzburg–Landau equations. Here is a generalized relative of it.

eq7[d0_, cList_] := 
D[A[x, t], t] == d0 D[A[x, t], x, x] +
                 A[x, t] (Abs[A[x, t]]^Range[0, Length[cList] - 1]).cList

We start with a localized smooth initial condition. The function solveAndMakeGraphics7 solves the PDE and generates a contour plot of the absolute value of its solution.

solveAndMakeGraphics7[{d0_, cList_}, ppL_, do_] := 
Module[{T = 80, L = 60, ndsol, arg = If[# == 0., 0, Arg[#]]&,
        cOpts = Sequence[PlotPoints -> 200, Contours -> 30, 
        ContourLines -> False, FrameTicks -> False]},
 (* solve PDE *) ndsol = 
 NDSolve[{eq7[d0, cList], (* localized initial condition *)
          A[x, 0] == If[Abs[x/L] < 1/2, Cos[Pi x/L]^2, 0],  
          A[-L, t] == A[L,  t]}, {A}, {x, -L, L}, {t, 0, T}, 
          MaxSteps -> 10^5, PrecisionGoal -> 3, AccuracyGoal -> 3, 
         Method -> {"MethodOfLines", Method -> StiffnessSwitching,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL}, "MinPoints" -> {ppL}}}];
(* visualize solution *)
Block[{$DisplayFunction = Identity},
{(* absolute value *)
 ContourPlot[Evaluate[Abs[A[x, t] /. ndsol[[1]]]], {x, -L, L}, {t, 0, T}, 
            Evaluate[cOpts], ColorFunction -> (Hue[0.8 #]&),
            PlotRange -> All],
 (* argument *)
 ContourPlot[Evaluate[arg[A[x, t] /. ndsol[[1]]]], {x, -L, L}, {t, 0, T}, 
            Evaluate[cOpts], ColorFunction -> Hue,
            PlotRange -> {-Pi, Pi}]}]]

Here are some of the resulting solutions. We show the absolute value and the argument for six sets of parameters. Once more we use a moderate number of discretization points and the lowest order to calculate a qualitatively correct solution. Using more points (and/or other methods) will give quantitatively more reliable solutions.

Off[NDSolve::eerr]; 

Show[GraphicsArray[Transpose[#]]]& /@ 
Partition[solveAndMakeGraphics7[#, 200, 1]& /@ 
(* parameters found by random searching *) {
{ 0.562585, { 0.472726 + 0.331368 I,  0.801779 - 2.480545 I,  0.474796 - 
              0.776168 I, -2.519074 - 0.952473 I, -1.014533 + 0.210578 I}},
{ 1.112343, { 0.333686 - 0.048153 I,  0.051441 + 0.346931 I,  0.262625 + 
              0.254911 I,  0.059902 + 0.067229 I, -0.125889 + 0.072792 I}},
{-0.457228, { 0.064266 + 0.219516 I,  0.103752 + 0.492474 I,  0.042513 - 
              0.062272 I, -0.190156 - 0.219106 I, -0.182677 - 0.551021 I}},    
{ 5.078614, {-0.043256 + 0.262005 I,  0.269991 + 0.278785 I,  0.371664 - 
              0.118355 I, -0.573401 + 0.766421 I,  0.132533 + 0.425168 I}},  
{ 0.896841, {-0.004116 - 0.029930 I,  0.289541 - 0.120667 I, -0.230067 - 
              0.511332 I,  2.387745 - 0.110776 I, -1.770957 - 1.860567 I}},
{ 0.613237, { 0.188025 + 1.353056 I, -0.093262 - 1.959692 I,  0.977863 - 
              1.689762 I,  1.183884 + 0.821903 I, -2.001599 + 0.123309 I}}}, 3]
                               
On[NDSolve::eerr];

Our next example is a system of two nonlinear 1+1D PDEs [✶].

eq8[A_, B_, δ_] = 
{D[u[x, t], t] == D[u[x, t], x, x] - u[x, t] v[x, t]^2 + A (1 - u[x, t]),
 D[v[x, t], t] == δ^2 D[v[x, t], x, x] + u[x, t] v[x, t]^2 - B v[x, t]};

This time we start with a smooth localized bump in the center of the integration domain.

startBump[x_, α_, L_] := If[Abs[x/L] < 1/α, Cos[Pi/2 (x/L)/(1/α)]^3, 0]

Like always, the function solveAndMakeGraphics8 solves the PDEs and generates a graphic of the solution, concretely, a contour plot of the second component.

solveAndMakeGraphics8[{{A_, B_, δ_}, {L_, T_}, α_, do_}, ppL_] := 
Module[{ndsol}, 
(* solve PDEs *)
ndsol = NDSolve[{eq8[A, B, δ], 
                 v[x, 0] == startBump[x, α, L], u[x, 0] == 1 - startBump[x, α, L], 
                 Derivative[1, 0][u][-L, t] == 0, Derivative[1, 0][u][+L, t] == 0,
                 Derivative[1, 0][v][-L, t] == 0, Derivative[1, 0][v][+L, t] == 0},
         {u, v}, {x, -L, L}, {t, 0, T}, 
         MaxSteps -> 10^5, PrecisionGoal -> 2, AccuracyGoal -> 2, 
         Method -> {"MethodOfLines", Method -> BDF,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL}, "MinPoints" -> {ppL}}}];
(* visualize solution *)
ContourPlot[Evaluate[v[x, t] /. ndsol[[1]]], {x, -L, L}, {t, 0, T}, 
       PlotRange -> All, PlotPoints -> 400, Contours -> 20, 
       ContourLines -> False, FrameTicks -> False, 
       ColorFunction -> (Hue[0.8 #]&)]]

Here are three examples. The special property of these PDEs is that the localized initial structure repeatedly copies itself [✶], [✶], [✶], [✶]. (For a variety of examples of such structure in expanding domains L=L(t), see [✶], [✶], [✶].)

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
      solveAndMakeGraphics8[#, 2400]& /@ 
       (* three sets of parameters *)
       {{{0.011051, 0.021635, 0.026161}, { 8, 4000}, 16, 1},
        {{0.373018, 0.095230, 0.003531}, { 5, 2000}, 29, 1},
        {{0.591115, 0.273925, 0.061422}, {10,  300}, 27, 1}}]]]

We continue with a system of two coupled nonlinear PDEs in 2+1 dimensions [✶], [✶].

eqs9[ϵ_, α_, A_] = 
{D[θ[x, y, t], t] == (ϵ^2 (D[θ[x, y, t], x, x] + D[θ[x, y, t], y, y]) -
                      θ[x, y, t]^3 + θ[x, y, t] + η[x, y, t])/α,
 D[η[x, y, t], t] == (D[η[x, y, t], x, x] + D[η[x, y, t], y, y]) -
                      θ[x, y, t] - η[x, y, t] + A};

This time, we start with a ring-shaped initial condition.

θ0[x_, y_, θ00_, L_] := (* ring-shaped elongation *)
With[{γ = 1/8}, If[L/2 - γ L <= Sqrt[x^2 + y^2] <= L/2 + γ L, 
                   Sin[3 Pi x/L] Cos[Sqrt[x^2 + y^2]/L/γ Pi/2]^2, 0]]

solveAndMakeGraphics9[{{ϵ_, α_, A_}, {L_, T_}, θ00_}, do_, ppL_] := 
Module[{ndsol},
 (* solve PDEs *)
 ndsol = NDSolve[{eqs9[ϵ, α, A], 
                 θ[x, y, 0] == θ0[x, y, θ00, L], η[x, y, 0] == 0,
                 θ[-L, y, t] == 0, θ[L, y, t] == 0, 
                 θ[x, -L, t] == 0, θ[x, L, t] == 0, 
                 η[-L, y, t] == 0, η[L, y, t] == 0, 
                 η[x, -L, t] == 0, η[x, L, t] == 0},
         {θ, η}, {x, -L, L}, {y, -L, L}, {t, T, T}, MaxSteps -> 10^5,
         PrecisionGoal -> 2, AccuracyGoal -> 2, 
         Method -> {"MethodOfLines", Method -> Automatic,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL, ppL}, "MinPoints" -> {ppL, ppL}}}];
 (* visualize η part of the solution *)
 ContourPlot[Evaluate[Abs[η[x, y, T] /. ndsol[[1]]]], {x, -L, L}, {y, -L, L}, 
       PlotRange -> All, PlotPoints -> 80, Contours -> 20, 
       ContourLines -> False, FrameTicks -> False, 
       ColorFunction -> (Hue[0.8 #]&)]]

For many parameter values the initial ring disintegrates into closed, deformed pieces. Here are four examples.

Show[GraphicsArray[#]]& /@ 
Block[{$DisplayFunction = Identity}, 
 Partition[solveAndMakeGraphics9[#, 8, 120]& /@ 
   (* four sets of parameter values *)
   {{{ 0.500521, 0.798381, -0.0086344}, { 22, 68}, -1.666285},
    {{-0.120092, 0.084926,  0.2385211}, { 13, 77}, -0.878614},
    {{ 0.503609, 1.048002, -0.0223839}, { 25, 42}, -1.632547},
    {{-1.138091, 0.023999, -0.0019694}, {173, 16} ,  1.704992}}, 2]]

We continue with a slightly nonlinear wave equation that models two different media [✶], [✶].

(* position-dependent c and ω *)
eq10[{c1_, c2_}, {ω1_, ω2_}] =
D[u[x, t], t, t] - If[x > 0, c1^2, c2^2] D[u[x, t], x, x] +
                   If[x > 0, ω1^2 Tan[u[x, t]], ω2^2 Sin[u[x, t]]] == 0;

(* two localized initial packets *)
startPackets[{𝒶_, 𝒷_}, x_] := 
 Which[Abs[x - 1/3] <= 𝒷, 𝒶 Cos[Pi/2 (x - 1/3)/𝒷]^2,
       Abs[x + 1/3] <= 𝒷, 𝒶 Cos[Pi/2 (x + 1/3)/𝒷]^2,
       True, 0];

solveAndMakeGraphicsTwoMediaWaveEquation[{{{c1_, c2_}, {ω1_, ω2_}}, 
                                 {𝒶_, 𝒷_}, {L_, T_}, g_, do_}, ppL_] := 
Module[{ndsol}, 
(* solve PDEs *)
ndsol = NDSolve[{eq10[{c1, c2}, {ω1, ω2}], 
                 u[x, 0] == startPackets[{𝒶, 𝒷}, x/L], 
                 Derivative[0, 1][u][x, 0] == 0, 
                 u[ L, t] == 0, u[-L, t] == 0},
         {u}, {x, -L, L}, {t, 0, T}, 
         MaxSteps -> 10^5, PrecisionGoal -> g, AccuracyGoal -> g, 
         Method -> {"MethodOfLines", Method -> Automatic,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL}, "MinPoints" -> {ppL}}}];
(* visualize solution *)
ContourPlot[Evaluate[Abs[u[x, t]] /. ndsol[[1]]], {x, -L, L}, {t, 0, T}, 
       PlotRange -> All, PlotPoints -> 240, Contours -> 20, 
       ContourLines -> False, FrameTicks -> False, 
       ColorFunction -> (Hue[0.8 #]&)]]

We clearly see reflections of the two spreading packets at the interface z=0.

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
solveAndMakeGraphicsTwoMediaWaveEquation[
                    {#1, #2, {10, 12}, 3, Pseudospectral}, 81]& @@@
  (* two parameter sets *)                  
  {{{{0.51314, 0.42573}, {1.5459, 0.0775}}, { 1.9945, 0.11295}},
   {{{0.09219, 0.38770}, {1.3690, 1.8004}}, {-0.6500, 0.082157}}}]]]

Our next example is the classical scalar wave equation with a piecewise constant speed of light. But now we consider the 2+1 case. The two materials are located above and below the line z≷0.

(* 2D wave equation with position-dependent c *)
eq11[c_] = If[z > 0, 1, 1/c^2] D[u[x, z, t], t, t] == 
               D[u[x, z, t], x, x] + D[u[x, z, t], z, z];

We start with a simple localized wave packet [✶]. (For other packets, see [✶], [✶], [✶], [✶].) To avoid the use of complex arithmetic, instead of a factor exp(i k.r), we will use two packets of the form cos(k.r) and sin(k.r) and add their contributions at the end. We then display the square of the absolute value of the field.

(* field and derivative at t == 0 *)
startPacket2D[{σl_, σr_, α_, ℓ_, k0_, β_}, {x_, z_}] := 
Module[{x0 = ℓ Cos[Pi - α], z0 = ℓ Sin[Pi - α], 
       𝓀0 = k0 {Cos[-α], Sin[-α]}, 
       ξD = {Cos[-α], -Sin[-α]}, ηD = {Sin[-α], Cos[-α]}, ξd, ηd},
       ξd = ({x, z} - {x0, z0}).ξD; ηd = ({x, z} - {x0, z0}).ηD;
      (* elongated moving Gaussian packet *)
      Exp[-ξd^2/σr^2 - ηd^2/σl^2] Cos[𝓀0.{x, z} + β]]   
      
startPacket2DDerivative[{σl_, σr_, α_, ℓ_, k0_, β_}, {x_, z_}] := 
Module[{t}, D[startPacket2D[{σl, σr, α, ℓ, k0, β}, 
                     {x - Cos[α] t, z + Sin[α] t}], t] /. t -> 0]

This time, we repeatedly solve the wave equation over a certain time interval and create a visual snapshot of the state and the end of this time interval. This method allows to follow the wave packet for a longer time without saving the potentially quite large interpolating function generated by NDSolve.

solveAndMakeGraphicsListTwoMediaWaveEquation2D[
     {c_, {σl_, σr_, α_, ℓ_, k0_}, g_, do_}, {L_, T_}, ppL_, frames_] := 
Module[{ΔT = T/frames, state, pps, xy0 = L {Cos[Pi - α], Sin[Pi - α]},
        stateC, stateS}, 
(* set up PDE;
   use two initial conditions with shifted modulation to model
   a complex field  *)
{stateC, stateS} = NDSolve`ProcessEquations[{eq11[c], 
                 u[x, z, 0] == startPacket2D[{σl, σr, α, ℓ, k0, #}, {x, z}], 
                 Derivative[0, 0, 1][u][x, z, 0] == 
                 startPacket2DDerivative[{σl, σr, α, ℓ, k0, #}, {x, z}], 
                 u[ L, z, t] == 0, u[-L, z, t] == 0,
                 u[x,  L, t] == 0, u[x, -L, t] == 0},
         {u}, {x, -L, L}, {z, -L, L}, t, 
         MaxSteps -> 10^5, PrecisionGoal -> g, AccuracyGoal -> g, 
         Method -> {"MethodOfLines", Method -> Automatic,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL, ppL}, 
                     "MinPoints" -> {ppL, ppL}}}][[1]]& /@ {0, Pi/2};
(* propagate solution in time and make visualizations *)
Table[(* solve wave equation in time steps ΔT *)
   NDSolve`Iterate[stateC, j ΔT]; NDSolve`Iterate[stateS, j ΔT]; 
   (* form current solution *)
   {ppsC, ppsS} = NDSolve`ProcessSolutions[#, "Forward"]& /@ {stateC, stateS};   
   (* visualize solution *)
   Show[Graphics[ContourPlot[Evaluate[Abs[ppsC[[1, 2]]]^2] + Abs[ppsS[[1, 2]]]^2, 
                 {x, -L, L}, {z, -L, L}, Contours -> 20, 
                 ContourLines -> False, DisplayFunction -> Identity,
                 FrameTicks -> False, ColorFunction -> (GrayLevel[1 - #]&),
                 PlotRange -> All, PlotPoints -> 160, Epilog -> 
       (* lines along the reflected and refracted rays *)
       {{Hue[0], Thickness[0.008], Line[{{-L, 0}, {L, 0}}]},
        {GrayLevel[0], Thickness[0.004], 
         Line[2 {xy0, {0, 0}, L {Cos[α], Sin[α]}}]}}]],
            PlotRange -> L {{-1, 1}, {-1, 1}}, Frame -> True], {j, frames}]]

The following two rows of graphics show the temporal evolution of two packets impinging on the surface under different angles. The material on the left has a smaller phase velocity than the material on the right hand side. In the left row, one sees that a large fraction of the initial wave packet gets reflected. The two black lines represent the classical reflection law for the center of the wave packet. For larger times, the reflected wave packet is clearly to right of the ray prediction (analog to the Goos–Hänchen effect in optics [✶], [✶], [✶], [✶]).

Show[GraphicsArray[#]]& /@ Transpose[
 solveAndMakeGraphicsListTwoMediaWaveEquation2D @@@
     {{{1/Sqrt[4], {1/2, 1/2, Pi/8, 3, 12}, 3, 1}, {4, 7}, 140, 5},
      {{1/Sqrt[2], {1/2, 1/2, Pi/3, 3, 16}, 3, 1}, {4, 7}, 140, 5}}]

Show[GraphicsArray[#]]& /@ Transpose[
 solveAndMakeGraphicsListTwoMediaWaveEquation2D @@@
     {{{1/Sqrt[4], {1/2, 1/2, Pi/8, 3, 12}, 3, 1}, {4, 7}, 140, 24},
      {{1/Sqrt[2], {1/2, 1/2, Pi/3, 3, 16}, 3, 1}, {4, 7}, 140, 24}}]

For refraction on negative n materials, see [✶], [✶]. For a nonlinear reaction-diffusion system with two different media, see [✶].

We continue with the local induction approximation PDE for vortices [✶], [✶], [✶]. Given a curve X(s), where s is the arc length, the time evolution of the curve is assumed to be

((∂X(t,s))/(∂t))=(Subscript[X, s])(t,s)×(Subscript[X, s s])(t,s).

The length of a curve stays constant under this time-evolution.
Using any parametrization instead of the arc length, we obtain the similar PDE

((∂X(t,σ))/(∂t))=(((Subscript[X, σ])(t,σ)×(Subscript[X, σ σ])(t,σ))/(Sqrt[((((Subscript[X, σ])(t,σ).(Subscript[X, σ])(t,σ)))^3)])).

We define the function solveAndMakeGraphicsLIAEquation to solve the solution numerically and to display the time-evolved curve as a surface sept out.

X[t_, σ_] := {x[t, σ], y[t, σ], z[t, σ]};

eq12[α_, {t_, σ_}] = 
Thread @ (D[X[t, σ], t] == α Together @
          Cross[1/Sqrt[D[X[t, σ], σ].D[X[t, σ], σ]] D[X[t, σ], σ],
                1/Sqrt[D[X[t, σ], σ].D[X[t, σ], σ]]*
                 D[1/Sqrt[D[X[t, σ], σ].D[X[t, σ], σ]] D[X[t, σ], σ], σ]]);

solveAndMakeGraphicsLIAEquation[{curve_, σ_}, α_, ppd_, T_, ppg_, opts___] :=
Module[{ndsol, Τ},
   ndsol = (* solve local induction approximation PDE *)
    NDSolve[Flatten[{eq12[1, {t, σ}], Thread[X[0, σ] == curve],
                                      Thread[X[t, 0] == X[t, 2Pi]]}],
         {x, y, z}, {t, 0, T}, {σ, 0, 2Pi},
         MaxSteps -> 10^5, PrecisionGoal -> 4,
         (* use "Pseudospectral" because of periodic boundary conditions *)
         Method -> {"MethodOfLines", "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral", 
                    "MaxPoints" -> {ppd}, "MinPoints" -> {ppd}}}];
   Τ = ndsol[[1, 1, 2, 1, 1, 2]]; 
   (* visualize evolved curve as a parametrized surface *)
   Show[ParametricPlot3D[Evaluate[Append[X[t, σ] /. ndsol[[1]], 
                               {EdgeForm[Thickness[0.001]],
                                SurfaceColor[Hue[0.8 t/Τ], Hue[0.8 t/Τ], 2]}]],
                 {σ, 0, 2Pi}, {t, 0, Τ}, PlotPoints -> ppg,
                 DisplayFunction -> Identity] /. Polygon[l_] :> 
                 Polygon[Plus @@@ Partition[Append[l, First[l]], 2, 1]/2],
     opts, DisplayFunction -> $DisplayFunction, 
     PlotRange -> All, Axes -> False, Boxed -> False]]

Here are the evolution of two knots.

solveAndMakeGraphicsLIAEquation[
 {{Sin[σ] + 2 Sin[2σ], Cos[σ] - 2 Cos[2σ], Sin[3σ]}, σ}, 1, 120, 20, {90, 60}]

solveAndMakeGraphicsLIAEquation[
 {{Sin[σ] + 3/2 Sin[3σ], Cos[σ] - 3/2 Cos[3σ], Sin[4σ]}, σ}, 1, 120, 10, {180, 50}]

And here is the evolution of a Lissajou curve.

solveAndMakeGraphicsLIAEquation[
 {{Sin[2σ], Cos[3σ], Cos[σ]}, σ}, 1, 120, 3, {120, 50}]

By letting the curve change its topology when two curve segments are near to each other, complicated vertex tangles can be formed [✶], [✶], [✶].

We continue with a partial differential equation that contains mixed spatial and temporal derivatives. Solving such differential equations is currently not yet supported. We start with a relatively simple equation, a Born–Infeld-type wave equation [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶].

(* simplest Born–Infeld wave equation for a scalar field *)
eqBornInfeld[] := (1 + D[ϕ[x, t], x]^2) D[ϕ[x, t], t, t] - 
                   2 D[ϕ[x, t], x] D[ϕ[x, t], t] D[ϕ[x, t], x, t] -
                  (1 - D[ϕ[x, t], t]^2) D[ϕ[x, t], x, x] == 0
 
(* two smooth initial wave packets *)                                  
startPackets[β_, x_] := 
 Which[Abs[x - 1/3] <= β, Cos[Pi/2 (x - 1/3)/β]^2,
       Abs[x + 1/3] <= β, Cos[Pi/2 (x + 1/3)/β]^2, True, 0];                                  

solveAndMakeGraphicsBornInfeld[{β_, {L_, T_}, g_, do_}, ppL_] := 
Module[{ndsol}, 
(* solve PDEs *)
ndsol = NDSolve[{eqBornInfeld[], ϕ[x, 0] == startPackets[β, x/L], 
                 Derivative[0, 1][ϕ][x, 0] == 0,
                 ϕ[ L, t] == ϕ[-L, t]},
         {ϕ}, {x, -L, L}, {t, 0, T}, 
         MaxSteps -> 10^5, PrecisionGoal -> g, AccuracyGoal -> g, 
         Method -> {"MethodOfLines", Method -> Automatic,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL}, "MinPoints" -> {ppL}}}];
(* visualize solution *)
ContourPlot[Evaluate[ϕ[x, t] /. ndsol[[1]]], {x, -L, L}, {t, 0, T}, 
       PlotRange -> All, PlotPoints -> 240, Contours -> 20, 
       ContourLines -> False, FrameTicks -> False, 
       ColorFunction -> (Hue[0.8 #]&)]]
       
(* example solution *)       
solveAndMakeGraphicsBornInfeld[{1/6, {10, 25}, 3, Pseudospectral}, 121]

(The relativistic diffusion equation also exhibits mixed spatial and temporal derivatives [✶].)

The next example generates graphics that show how an initially localized smooth perturbation splits into peakon trains [✶], [✶], [✶], [✶].

(* peakon PDEs *)
eqPeakon[u_, {x_, t_}, α_] := 
 Function[m, D[m, t] + u[x, t] D[m, x] + 2 D[u[x, t], x] m == 0][
                                   u[x, t] - α^2 D[u[x, t], x, x]];
(* smooth initial bump *)                                  
startBump[β_, x_] := If[Abs[x] >= β, 0, Cos[Pi/2 x/β]^2];      

solveAndMakeGraphicsPeakon[{{α_, β_}, {L_, T_}, g_, do_}, ppL_] := 
Module[{ndsol, u}, 
(* solve PDEs *)
ndsol = NDSolve[{eqPeakon[u, {x, t}, α], u[x, 0] == startBump[β, x/L], 
                 u[ L, t] == u[-L, t]},
         {u}, {x, -L, L}, {t, 0, T}, 
         MaxSteps -> 10^5, PrecisionGoal -> g, AccuracyGoal -> g, 
         Method -> {"MethodOfLines", Method -> Automatic,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                     "MaxPoints" -> {ppL}, "MinPoints" -> {ppL}}}];
(* visualize solution *)
ContourPlot[Evaluate[u[x, t] /. ndsol[[1]]], {x, -L, L}, {t, 0, T}, 
       PlotRange -> All, PlotPoints -> 240, Contours -> 20, 
       ContourLines -> False, FrameTicks -> False, 
       ColorFunction -> (Hue[0.8 #]&)]]
       
(* solutions for two parameter sets *)       
Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
 {solveAndMakeGraphicsPeakon[{{-1.61667, 0.234245}, {121, 80}, 3, 8}, 80],
  solveAndMakeGraphicsPeakon[{{+1.31838, 0.624857}, { 62, 80}, 3, 8}, 80]}]]]

For further partial differential equations that contain mixed spatial and temporal derivatives, see [✶] (describing a the motion of a vibrating hanging shoe string).