Solving  Evolution Equations with Mathematica

Rob Knapp
Wolfram Research

Abstract

The NDSolve command in Mathematica has a framework for solving evolution equations using the numerical method of lines.  In this presentation, the developer of the framework will explain how it works and then show several examples of how to tune the method to get better solutions in less time.

Initialization

Off[General::spell];
Off[General::spell1];

TalkDirectory  =  "C:\\Web\\RD_HTML\\RKnapp\\Talks\\India"

\(\($Path\  = \ 
        Join[$Path, \ 
          Map[ToFileName[
                TalkDirectory, \ #] &, \ {"Data", \ "Packages"}]];\)\)

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

Needs["GUIKit`"]

monitorwidget = Widget["Frame", {
  "title" -> "Solution Monitor",
  WidgetSpace[5],
  { 
    WidgetSpace[5],
      Widget["TextField", {"text" -> "", "editable"->False}, Name->"time"],
     WidgetFill[],
     WidgetSpace[5], Widget["Button", {"text" -> "Stop",
      BindEvent["action",
        Script[
          stopRequested = True;
          ]
       ]}],
       WidgetSpace[5]},

  Widget["ImageLabel", {
    "preferredSize" -> Widget["Dimension", {"width"->288, "height"->288}]
    }, WidgetLayout -> {"Stretching" -> {Maximize, Maximize}}, Name -> "imagePanel"],
    
  Script[
    stopRequested = False;
    ]
    
  }, Name -> "monitorFrame"];

textmonitorwidget := Widget["Frame", {
  "title" -> "Solution Monitor",
  WidgetSpace[5],
  { 
    WidgetSpace[5],
      Widget["TextArea", {"text" -> "", "editable"->True, "rows"->4, "columns"->20}, Name->"time"],
     WidgetFill[],
     WidgetSpace[5], Widget["Button", {"text" -> "Stop",
      BindEvent["action",
        Script[
          stopRequested = True;
          ]
       ]}],
       WidgetSpace[5]},
    
  Script[
    stopRequested = False;
    ]
    
  }, Name -> "monitorFrame"];

The numerical method of lines (MOL) is a technique for solving partial differential equations that works by discretizing in all but one dimension, and then integrating the semi-discrete problem as a system of ODEs or DAEs.

Burgers` equation,

is a simplified model for some aspects of gas dynamics.  When the viscosity parameter, ν,  is small, the nonlinear wave speed of the equation leads to the formation of sharp fronts. It is a good example to consider since even though the formation of fronts can lead to numerical challenges, there is a transformation (Cole-Hopf) that can be made to find analytic solutions such as

These can be used to check the accuracy of numerical approximations.

The exact solution u(x, t) = 
        (1/(1 + (e^(((2 x  -  t)/(4 ν)))))) is effectively a front which keeps a consistent shape and moves to the left with increasing t.

(Subscript[u, exact])  =  
        Function[{t, 
            x},   (1/(1  +  Exp[((2 x  -  t)/(4 ν))]))];

Block[{ν  =  0.01}, 
        Plot3D[(Subscript[u, exact])[t, x], {x, 0, 1}, {t, 0, 2}]];

Here is the equation and boundary conditions as Mathematica input.

burgers  =  
        D[u[t, x], t]  ==   ν D[u[t, x], x, x]  -  
            D[u[t, x]^2, x]/2 ;
binit  =  
        u[1/2, x]  ==  
          (1/(1  +  Exp[((  2  x  -  1/2 )/(4  ν))]));
bbvals  =  {u[t, 0]  ==  
            (1/(1  +  Exp[- (( t )/(4  ν))])),  
          u[t, 1]  ==  
            (1/(1  +  Exp[(( 2 - t )/(4  ν))]))};

The next two commands compute a numerical solution for a particular value of ν and plots the error of the numerical solution as a function of x for a particular t.

(Subscript[u, nds])  =  
      Block[{ν  =  0.01}, 
        First[u  /.  
            NDSolve[{burgers, binit,  bbvals}, 
              u,  {t, 1/2, 3/2},  {x, 0, 1}]]]

Show[
        Block[{ν  =  0.01,  
            t  =  1,  $DisplayFunction  =  Identity}, {Plot[
              (Subscript[u, exact])[t, x] -  (Subscript[u, nds])[t, x], {x, 0, 1},  
              PlotRange -> All],  X  =  (Subscript[u, nds])@"Coordinates"[2]; 
            ListPlot[Transpose[{X, (Subscript[u, exact])[t, X]  -  (Subscript[u, nds])[t, X]}],  
              PlotStyle -> RGBColor[1, 0, 0]]}]];

There are two possibly significant sources of the error you see above, the spatial discretization error and the temporal discretization error.  By varying discretization parameters it is possible to see which one has the largest contribution.

With a particular grid, it is relatively easy to create a discretization by using replacement rules to replace the continuous function u(t, x) with a discretized vector function U(t).

Consider a uniform grid with  n  +  1 points and spacing (Subscript[h, n])  =  1/n for a particular n.

n  =  10; (Subscript[h, n])  =  (1/n);

Let U be the vector function (of t) of the discretized values.

U[t_]  =  Table[(Subscript[u, i])[t], {i, 0, n}];

Ultimately, the goal is to enable the use of difference formulas of any asymptotic order, but at this stage, to keep things simple, consider the standard second order difference formulas.

One easy way to implement such difference formulas in Mathematica is to use ListCorrelate with the following kernels

k[1]  =  {-1, 0, 1}/(2 (Subscript[h, n]));
k[2]  =  {1, -2,  1}/((Subscript[h, n])^2);
DiscreteD[j_][U_]  :=  
        Join[{0}, ListCorrelate[k[j],  U], {0}];

Note that this only discretizes the interior.  The boundary conditions will be handled separately.

The following rules implement the (semi-)discretization

drules  =  {u[t_, x] :>  
            U[t], 
Derivative[i_, 0][u][t, 
              x]  :>  Derivative[i][U][
              t], 
Derivative[0, j_][u][t, 
              x]  :>  DiscreteD[j][
              U[t]], 
x -> (Subscript[h, n]) 
              Range[0, n]};

ultimately, the goal will be to enable difference formulas of any order.

deqns  =  Thread[burgers  /.  drules]

dinit  =  Thread[binit  /.  drules]

To handle the boundary conditions, we can replace the right hand sides of first and last equations with the time derivatives of the boundary conditions

deqns[[{1, -1},  2]]  =  D[bbvals,  t][[All, 
          2]]

These ODE equations are in a form which can be solved with NDSolve

lines  =  
        Block[{ν  =   .1}, 
          NDSolve[{deqns,  dinit},  U[t], {t, 1/2, 3/2}]];

lineplot  =  
        Show[
Block[{$DisplayFunction  =  
                Identity}, 
Table[
              ParametricPlot3D[{i (Subscript[h, n]), t,  First[(Subscript[u, i])[t]  /. lines]}, {t, 
                  1/2, 3/2}], {i, 0, n}]],  PlotRange -> All,  
          AxesLabel -> {"x", "t", "u"}];

To compute with a different grid, you can change the definition of n above and run through the discretization again.  This is not very adaptable in that to change the grid you have to work through the entire discretization again.  Using a symbolic representation for the grid allows going to the next level.

A function which can generate finite difference formulas for any order derivative to any degree of approximation for spatial derivatives makes the process of discretization easier and more general.  Using algorithms for generating weights, such as that of Fornberg, it is possible, though not simple, to write such a function in Mathematica.

That is not necessary, since in Mathematica 5, there is a function, NDSolve`FiniteDifferenceDerivative which does this on a grid.

fdd2  =  
      NDSolve`FiniteDifferenceDerivative[
        Derivative[2],  (Subscript[h, n]) Range[0, n]]

This effectively generates a function which acts on values defined on a grid:

TableForm[fdd2[U[t]]]

At the edges of the grid, one-sided derivatives are used.

The default is to use fourth order differences, however, the order can be changed by using an option

TableForm[NDSolve`FiniteDifferenceDerivative[
          Derivative[2], (Subscript[h, n]) Range[0, n],  DifferenceOrder -> 2][
        U[t]]]

More detail can be found in the section on Finite Differences in the NDSolve  Advanced Documentation in the Mathematica help browser.

A more complete symbolic discretization can be achieved by using a function that symbolically represents the finite difference object and is not evaluated until the grid and difference order are specified.  Let X represent the grid symbolically and U represent the discretized function.

Clear[U,  X];

Then it is simple to define rules which implement the discretization

rhsrules =  {
u[t_, x] :>  
            U, 
Derivative[0, j_][u][t, 
              x]  :> FiniteDifference[j,  X,  m][
              U], 
x ->  
            X
};

drhs  =  burgers[[2]]  /.  rhsrules

The boundary conditions can be handled much as before. Note the use of an intermediate variable to ease the replacement of the appropriate parts.

drhsfun  =  
      Function[{t, U},  
          Block[{U1  =  drhs1}, U1[[{1, -1}]]  =  dbvals;  
            U1]]  /.  {drhs1 -> drhs,  
          dbvals  ->  D[bbvals,  t][[All, 2]]}

This function is not sufficient as is for a proper symbolic discretization because, evaluation with symbolic arguments t and U will lose necessary information.  A wrapper needs to be added which will prevent evaluation for symbolic t and U.  The solution is to use downvalues, for example

drhsd[t_?NumberQ,  U_?VectorQ]  :=  drhsfun[t, U]

However, the problem with this is that the grid dependent values (in particular the finite difference weights) have to be recomputed each time the function is evaluated.  What is needed is to get the body of the Function above into the downvalues.  Doing this is a little bit tricky because of the need to prevent evaluation.  One handy way to do it is to take advantage of the way that patterns substitute in Mathematica

makedrhs[
        fun_]  := (drhsd[t_?NumberQ,  U_?VectorQ]  := 
          fun[t, U])

The desired downvalue is created simply by calling this on the previously produced Function:

makedrhs[drhsfun]

You can check that the definition is correct by using the DownValues command

DownValues[drhsd]

You can see this is not quite adequate since what we really want to do is set up the downvalue once the grid and difference order are defined.  This, however, can be done with one more layer of rules:

makespecificdrhs[fun_,  XX_?VectorQ,  mm_Integer]  :=  
      Module[{orders,  rules},  
orders  =  
          Union[Cases[fun,  FiniteDifference[j_, __] ->  j, 
              Infinity,  
              Heads -> True]]; 
rules  =  
          Map[ (FiniteDifference[#,  X,  m] ->  
                  NDSolve`FiniteDifferenceDerivative[Derivative[#],  XX,  
                    DifferenceOrder -> mm]) &,  
            orders]; 
makedrhs[(fun  /. 
                rules)  /.  {X -> XX,  m -> mm}]]

Then, for example

(Subscript[h, n]) Range[0, n]

makespecificdrhs[drhsfun,  (Subscript[h, n]) Range[0, n],  6]

generates DownValues

DownValues[drhsd]

Producing an efficient evaluation for a vector U.

All that is needed for the complete symbolic discretization is inclusion of the initial conditions, which can readily be included in a function which generates the equations to integrate.

sdisc[XX_,  
        mm_]  :=  (
makespecificdrhs[drhsfun,  XX,  
          mm]; 
{U'[t]  ==  
            drhsd[t, 
              U[t]], (binit  /.   u[t_,  x] :>   U[
                    t])  /.  x -> XX})

Then, for example

sdisc[(Subscript[h, n]) Range[0, n],  2]

can be solved by NDSolve

sol  =  
      Block[{ν  =  0.1}, 
NDSolve[
          sdisc[(Subscript[h, n]) Range[0, n],  2],  U,  {t, 1/2, 3/2}]]

This is a vector-valued InterpolatingFunction object.

One interesting way to visualize it is to plot the data points.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

ifun  =  First[U  /.  sol];
times  =  
        First[InterpolatingFunctionCoordinates[
            ifun]];
values  =  
        InterpolatingFunctionValuesOnGrid[ifun];
X  =  N[(Subscript[h, n]) Range[0, n]];
grid  =  Outer[List,  times,  X];
points  =  
        MapThread[Point[Append[Reverse[#1],  #2]] &,  {grid,  values},  
          2];

Show[Graphics3D[points]];

There was little done in the above section which was special to Burgers` equation.  Basically, the only structure which was needed was an equation with time derivatives on the left, and initial and Dirichlet boundary conditions separated.  Periodic boundary conditions are even easier to handle because NDSolve`FiniteDifferenceDerivative has a PeriodicInterpolation option which can be set to true for this case.

Neumann and mixed boundary conditions can be handled by using the one-sided derivative approximations and Solve with the boundary conditions to find the appropriate differential equation to use at the boundary, but are beyond the scope of this paper.  The general technique is described in the Boundary Conditions section in the NDSolve  Advanced Documentation in the Mathematica help browser.

Thus, using the symbolic ideas developed above, it is possible to write a fairly general discretization function which transforms from a PDE to the system of ODEs resulting from discretization.  For brevity, the function is defined below with only minimal error checking.

NDSolve works using the generalized ideas with a few additions.  Most importantly, NDSolve sets up a symbolic spatial error estimate as a function of the grid and determining the grid so that the spatial error on the initial condition satisfies the AccuracyGoal and PrecisionGoal tolerances.  Description of the technique for obtaining error estimates is described in the Spatial Error Estimates section in the NDSolve  Advanced Documentation in the Mathematica help browser.

The systems of ODEs generated by the method of lines can be quite stiff.  
An example of an equation that leads to stiff ODEs is the KdV equation,

Lets consider it with periodic boundary conditions:

KdV  =  
        D[ψ[t, x],  
            t]  ==  -α ψ[t, x] D[ψ[t, x],  x] - 
            D[ψ[t, x],  {x, 3}];
KdVinit  =  ψ[0, x] ==  Φ[
            x];
KdVbc  =  ψ[t, 0] ==  ψ[t, 40];

An interesting initial condition to consider is a combination of solitons

S[{α_,  c0_,  x0_}, t_,  
          x_ ]  :=  ((3 c0)/α) (((Sech[ ((Sqrt[c0])/2) (x - x0 - c0 t)]))^2);

Since we are using periodic boundary conditions add up all the period contributions so that the initial condition really does appear periodic on the grid.

Φ[x_] = 
        S[{α, 3, 5}, 0,  x]  + 
          S[{α, 3/2, 15}, 0,  
            x] + 
S[{α, 3, 5 + 40}, 0,  x]  + 
          S[{α, 3/2, 15 + 40}, 0,  
            x] + 
S[{α, 3, 5 - 40}, 0,  x]  + 
          S[{α, 3/2, 15 - 40}, 0,  x];

sol  =  
      Block[{α  =  1}, 
        NDSolve[{KdV,  KdVinit,  KdVbc},  ψ,  {t, 0,  10}, {x, 0, 
            40}]]

The problems that show up are because of stiffness.

Do[Plot[First[ψ[t, x]  /.  sol],  {x, 0, 40},  
        PlotRange -> {0,  .1}],  {t, 0, 3.5,  .5}]

The default method for NDSolve in Mathematica switches to a stiff method when stiffness is detected.  The default stiff methods are BDFs of order 1-5 and are adequate for most problems (they are A-α stable).  However, for problems with large eigenvalues on the imginary axis like this, it is usually necessary to use an A-stable method, or one that includes the left hand side of the entire complex plane in its stability region.

BFDs of order 1, 2 are A-stable so one way to remove the stiffness problem is to reduce the order to 2

\(Timing[
      sol12\  = \ 
        Block[{α\  = \ 1}, 
          NDSolve[{KdV, \ KdVinit, \ KdVbc}, \ ψ, \ {t, 0, \ 10}, {x, 0, 
              40}, \ Method -> {"BDF", \ "MaxDifferenceOrder" -> 2}]]]\)

Plot3D[First[ψ[t, x]  /.  sol12],  {x, 0, 40}, {t, 0, 10},  
      Mesh -> False,  PlotPoints -> 200,  PlotRange -> All, 
      ViewPoint -> {0.134,  -2.274,  2.502}]

Linearly implicit extrapolation methods also provide A-stable methods.  The StiffnessSwitching method switches to an A-stable method when stiffness is detected.

Timing[
      solss  =  
        Block[{α  =  1}, 
          NDSolve[{KdV,  KdVinit,  KdVbc},  ψ,  {t, 0,  10}, {x, 0, 
              40},  Method -> "StiffnessSwitching"]]]

For this particular example, it takes longer, but that is because the integration of the problem is fairly easy to resolve, so the limitation of BDFs to second order is not tragic.

With this example, the instability showed up right away.  For some PDEs, it takes much longer to start creeping in.
Also, you may just want to check how the solution is evolving for a set of initial conditions or parameters and stop it when you are not getting what you hoped for.

Using GUIKit, it is not difficult to have Mathematica bring up a window that monitors the solution of the PDE.

Needs["GUIKit`"]

The widget defined below sets this up:

monitorwidget = Widget["Frame", {
  "title" -> "Solution Monitor",
  WidgetSpace[5],
  { 
    WidgetSpace[5],
      Widget["TextField", {"text" -> "", "editable"->False}, Name->"time"],
     WidgetFill[],
     WidgetSpace[5], Widget["Button", {"text" -> "Stop",
      BindEvent["action",
        Script[
          stopRequested = True;
          ]
       ]}],
       WidgetSpace[5]},

  Widget["ImageLabel", {
    "preferredSize" -> Widget["Dimension", {"width"->288, "height"->288}]
    }, WidgetLayout -> {"Stretching" -> {Maximize, Maximize}}, Name -> "imagePanel"],
    
  Script[
    stopRequested = False;
    ]
    
  }, Name -> "monitorFrame"];

In priniciple, the EvaluationMonitor or StepMonitor options of NDSolve could be used to give information to the monitor widget, but they are less than ideal for this purpose because the times at which the solution appears will be unpredictable.

A way to have the monitor show regular frames with spacing corresponding to the evolution of the temporal variable is to set up an event for detection by the EventLocator method.

The command below sets up an action to be taken every time the event Sin[π t/Δt] is zero-- i.e. every interval of 0.25 in t.
The action sends some text, the integration time and the memory used, and an image generated by a plot of the solution.
(Note it is usually desireable to set a PlotRange since otherwise the image will jump around.

Δt  =   .25; 
event  =  
        Sin[π t/Δt];
\(action\  := \ \((
timem\  = \ 
          StringJoin["t = ", \ ToString[t], \ "  RAM in use = ", 
            ToString[ .1 
                 Round[MemoryInUse[]/
                    10^5]], \ " MB"]; 
Sow[{t, \ 
            timem}]; 
monitor\ @\ 
          SetPropertyValue[{"time", \ "text"}, 
            timem]; 
monitor\ @\ 
          InvokeMethod[{"imagePanel", "setData"}, 
ExportString[\ 
              Plot[ψ[t, x], \ {x, 0, 40}, \ PlotRange -> {0, 10}, 
                DisplayFunction -> \ 
                  Identity], \ "GIF"]]; 
If[\ 
          TrueQ[\ monitor\ @\ Script[stopRequested]], \ 
          monitor\ @\ Script[stopRequested\  = \ False]; 
          Throw[Null, \ "StopIntegration"]];
)\)\)

monitor  =  GUIRun[monitorwidget]; 
    sol  =  Block[{α  =  1}, 
        NDSolve[{KdV,  KdVinit,  KdVbc},  ψ,  {t, 0,  10}, {x, 0, 
            40},  Method -> {"EventLocator",  "Event" -> 
                event,  "EventAction" :>  
                action}]]

One of the advantages of MOL is that it makes it relatively easy to change the spatial discretization. In this example, I'll show the effect of changing the asymptotic order of the spatial derivative approximation.

Here is an equation from Stephen's A New Kind of Science:

pde  =  
        D[u[t, x], t, t]  ==  
          D[u[t, x], x, 
              x]  +  (1  -  u[t, x]^2) (1  +  
                  8/3 u[t, x]);
init  =  {u[0, x]  ==  
            Exp[-x^2],  Derivative[1, 0][u][0, x]  ==  
            0};
    bc[
          L_]  :=  {u[t,  -L]  ==  u[t, L]};

A concise way to work with different options is to parameterize the whole problem.  
Here I will parameterize it by the width of the computational domain, L.
The goal is to work up to L = 50 or more, but for trial runs, L = 10 is a good value to start with.

problem[T_,  L_]  :=  
      Apply[Sequence,  {{pde,  init,  bc[L]},  
          u,  {t, 0, T},  {x, -L, L}}]

For now start with a small trial run:

T  =  10;
L  = T;

This solves it using the default settings of Mathematica.  We have set the step method at this point to make later comparison easier.
We choose an explicit method since the problem does not lead to any particular stiffness.  The main reason for this is that the timesteps required to resolve solution features are somewhat smaller than those required by the stability limit.

stepmethod  =  "ExplicitRungeKutta";

This uses the default settings for the method of lines.

Timing[
      sold  =  
        First[ u /.  
            NDSolve[problem[T, L],  Method -> stepmethod]]]

The reason the message is issued is because the solution has more rapidly varying structure than the intial condition on which the grid is based does.

Plot3D[sold[t, x], {t, 0, T}, {x, -L, L},  
        PlotPoints -> 250,  Mesh -> False];

When started from zero initial conditions, the equation has a peroidically varying solution which you can think of as a background.

v[t_]  :=  
        3 Sin[Sqrt[5/6] t]^2/(2  +  
                3 Cos[Sqrt[5/6] t]^2);

Subtracting off the background allows you to see the structural changes in the solution more easily.

A density plot makes it a bit easier to see some of the interesting structure that evolves.

DensityPlot[
        sold[-τ, x]  -  v[-τ], {x, -L, L}, {τ, 
          0, -T},  PlotPoints -> 250,  Mesh -> False,  
        ColorFunction -> Hue];

You can use the information from the message to estimate a better starting grid:
The default method is asymptotically fourth order.
The message reports what the scaled error was.
Thus, if nx1 was the points used to compute this, we should try nx2 = nx1*(15)^(1/4) ≃2*nx1

Load a package of functions for working with the data in InterpolatingFunction objects

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

This is the length of the spatial grid.

nx  = 
      Length[InterpolatingFunctionCoordinates[sold][[2]]]

nx4  =  2*nx  -  1

To keep the commands below more concise, I will set up the options for NDSolve that I want to change in a function that allows me to change the number of grid points and the difference order.
The precision and accuracy goals for the temporal integration as set much higher so that the main source of the error in the comparisons shown below is from the spatial discretization.

\(\(opts[n_, \ do_:  4]\  := \ 
        Apply[Sequence, \ 
{PrecisionGoal -> 8, \ 
            AccuracyGoal -> 
              8, \ 
Method -> {"MethodOfLines", 
Method -> 
                  stepmethod, \ 
"SpatialDiscretization" -> {"TensorProductGrid", 
"MinPoints" -> 
                      n, \ "MaxPoints" -> 
                      n, 
"DifferenceOrder" -> 
                      do, 
AccuracyGoal -> 4, \ 
                    PrecisionGoal -> 4}}}];\)\)

Solve with a specified number of points.

Timing[
      sol4  =  
        First[ u /.  NDSolve[problem[T, L],  opts[nx4]]]]

It turns out the the density plots look the same. The difference of two computed solutions often gives far more information.
It is easier for me to get information from 2 dimensional plots, so I'll be comparing solutions at t = T where the error will typically be the largest.

Plot[sold[T, x]  -  sol4[T, x],  {x, -L, L},  
        PlotRange -> All];

In some situations (not really in this one) , part of what you will be seeing in plots like the one above are the effects of interpolating between grid points.  A way to avoid this is to compare solutions at common grid points: this is why the value nx4 = 2*nx1 - 1 was chosen: all points on the coarser grid lie on the finer grid.

CompareOnGrid[s1_,  s2_,  opts___]  :=  
      Module[{tx1  = InterpolatingFunctionCoordinates[s1],  
          tx2  =  InterpolatingFunctionCoordinates[s2], diff,  grid,  
          Tend}, 
Tend  =  
          Min[Max[tx1[[1]]],  
            Max[tx2[[1]]]]; 
grid  =  
          If[Length[tx1[[2]]]  >  Length[tx2[[2]]],  
            tx2[[2]],  tx1[[2]]]; 
diff  =  
          s1[Tend, grid] -  s2[Tend,  grid]; 
Show[
          Block[{$DisplayFunction =  Identity}, {
ListPlot[
                Transpose[{grid,  diff}],  opts,  PlotJoined -> True, 
                PlotRange -> All], 
ListPlot[
                Transpose[{grid,  diff}],  opts,  
                PlotStyle -> PointSize[0.015], 
                PlotRange -> 
                  All]
}]]; 
Max[
          Abs[diff]]
]

CompareOnGrid[sold,  sol4]

One way to estimate the error is to use Richardson extrapolation:
Let s(h) be the solution computed with time step h
      ex(h) be the exact solution
      p be the method order such that  s(h) ≃ ex(h) + (ch^p)
 Then
     s(h) - s(h/2) ≃ (ch^p)(1  -  (2^(-p)))
 so
    c ≃ ((s (h)  -  s (h/2))/((h^p)(
                1  -  (2^(-p)))))

In the max norm, we have

\(c\  = \ 
      0.0012202601171609484`\/\(\(\((2\ L/nx)\)\^4\) \((15/16)\)\)\)

So the error on the finer grid is extrapolated to be about

c (((2  L/nx4))^4)

This can be verified by make a yet finer grid and doing the same process.

nx42  =  nx4*2  -  1;
Timing[
      sol42  =  
        First[ u /.  NDSolve[problem[T, L],  opts[nx42]]]]

CompareOnGrid[sol4,  sol42,  PlotRange -> All]

The extrapolation predicts c = 6.6.

Sometimes a cheaper way of improving the resolution of the solution is to increase the difference order instead of making the grid finer.

Timing[
      sol6  =  
        First[ u /.  NDSolve[problem[T, L],  opts[nx4,  6]]]]

CompareOnGrid[sol6,  sol42]

Timing[
      sol8  =  
        First[ u /.  NDSolve[problem[T, L],  opts[nx4,  8]]]]

CompareOnGrid[sol8,  sol42]

There are three problems with using arbitrarily high difference orders.

On a periodic grid, however, the limiting difference order can be shown to be equivalent to using a trigonometric basis and the derivative approximations for this basis can be computed using FFTs which effecitively remove the first two problems.

Timing[
      solps  =  
        First[ u /.  
            NDSolve[problem[T, L],  
              opts[nx4,  "Pseudospectral"]]]]

CompareOnGrid[solps,  sol8]

Now set the length of the computational domain to be 50.

T  =  50;
L  =  50;

Because of the increasing structure, the initial error estimate NDSolve computes will be inadequate.
We will compute two solutions and compare them.
The value nn is chosen here to be an efficient length for computing the FFTs required.

nn  =  8*7*9 + 1

\(Timing[
      sol\  = \ First[
          u\  /. \ NDSolve[{pde, \ init, \ bc[L]}, \ 
              u, \ {t, 0, T}, \ {x, \(-L\), L}, \ 
              Method -> 
{MethodOfLines, 
DiscretizedMonitorVariables -> 
                    True, 
Method -> {EventLocator, 
Method -> 
                        stepmethod, \ 
Event :> \ Abs[\(u[t, x]\)[\([1]\)]\  - \ v[t]]\  - \ 
                          0.01, 
EventAction :> 
                        Throw[Tend\  = \ 
                            t, \ "StopIntegration"]}, 
SpatialDiscretization -> 
{TensorProductGrid, 
DifferenceOrder -> "Pseudospectral", \ MinPoints -> nn, \ 
                      MaxPoints -> 
                        nn}
}
]
]]\)

\(Timing[
      sol2\  = \ 
        First[u\  /. \ 
            NDSolve[{pde, \ init, \ bc[L]}, \ 
              u, \ {t, 0, T}, \ {x, \(-L\), L}, \ 
              Method -> 
{MethodOfLines, 
DiscretizedMonitorVariables -> 
                    True, 
Method -> {EventLocator, 
Method -> 
                        stepmethod, \ 
Event :> \ Abs[\(u[t, x]\)[\([1]\)] - v[t]]\  - \ 
                          0.01, 
EventAction :> 
                        Throw[Tend\  = \ 
                            t, \ "StopIntegration"]}, 
SpatialDiscretization -> 
{TensorProductGrid, 
DifferenceOrder -> "Pseudospectral", \ MinPoints -> 2*nn\  + \ 1, \ 
                      MaxPoints -> 
                        2*nn + \ 
                          1}
}
]
]]\)

CompareOnGrid[sol, sol2]

DensityPlot[
        sol2[-τ, x]  -  v[-τ], {x, -L, L}, {τ, 
          0, -T},  PlotPoints -> 250,  Mesh -> False,  
        ColorFunction -> Hue];

This is a problem I was asked about by someone at the Lawrence Livermore National Laboratory in California.  He was trying to solve a PDE for which the solution gave a model for turbulent spectra.  Here is what he wrote:

I am trying to solve a pde for the time dependent energy spectrum E(k,t) of an evolving turbulent flow.  The pde comes from eqn (3.1) of D.C. Besnard et. al., "Spectral Transport Model for Turbulence," Theoretical and Computational Fluid Dynamics, Vol 8, pp. 1-35 (1996).  It looks like this :

(∂_t) E = -(Subscript[c, 1])(∂_k) ((k^(5/2))  (E^(3/2)))+ (Subscript[c, 2]) (∂_k) ((k^(7/2))  (E^(1/
                  2)) (∂_k) E)

where (Subscript[c, 1]) and (Subscript[c, 2]) are constants.  The following plot is the initial condition for the energy spectrum E(k,0) vs. wavenumber k.  It contains energy over a narrow band of wavenumbers.  Figure 1 of the above reference shows how this spectrum evolves over time.  I am hoping to start by reproducing this Figure using Mathematica.

Needs["Graphics`Graphics`"]

z[k_] = Exp[-0.5*(Log[10, k]/0.1)^2];
\(\(p1 = 
        LogLogPlot[z[k], {k, 10\^\(-2\), 10\^2}, 
          GridLines -> 
            None, 
		PlotRange -> {{10\^\(-3\), 10\^3}, {10\^\(-7\), 
                10\^1}}, PlotPoints -> 1000, 
		Frame -> True, 
          FrameLabel -> {\*""\!\(Log\_10\)(k)"", "E(k,0)"}, 
          RotateLabel -> False, AspectRatio -> 2\/3, 
          PlotLabel -> "Energy Spectra"];\)\)

I have been able to use the Mathematica PDE solver if I choose boundary conditions covering a rather narrow range of wavenumber 
0.01 < k < 50 for a relatively short period of time 0 < t < 100.  (What I really want is (10^(-4)) < k < (10^4) )

kmin = 0.01; kmax = 50; npts = 70;
klist = 
        Table[10^logk, {logk, Log[10, kmin], 
            Log[10, kmax], (Log[10, kmax] - Log[10, kmin])/
              npts}];

c1=1;c2=1;
f=NDSolve[{D[y[k, t], t] == -c1*D[(k^2.5)*((y[k, t])^1.5),k]+
c2*D[(k^3.5)*((y[k, t])^0.5)*D[(y[k, t]),k],k],
y[k, 0] == z[k],
y[kmin, t] == z[kmin],
y[kmax, t] == z[kmax]},
y, Join[{k},klist],{t, 0, 100},
MaxSteps->200, Method->StiffnessSwitching]

b  =  y  /.  First[f]

The following plot shows E(k,t) for t = 0, 0.1, 0.2, 0.7, 10 & 100.  This looks roughly correct (though not perfect), but when I try to increase the range of wavenumber over which I solve this pde, Mathematica gets hung up and never returns.   I can't seem to make it work.  Any suggestions?

\(\(LogLogPlot[{z, Re[b\ [k, 0.1]], Re[b\ [k, 0.2]], Re[b\ [k, 0.7]], 

          Re[b\ [k, 10]], Re[b\ [k, 100]]}\ , {k, kmin, 
          kmax}, 
		GridLines -> 
          None, 
		PlotRange -> {{10\^\(-2\), 10\^2}, {10\^\(-7\), 
              10\^1}}, 
        PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], 
            RGBColor[1, 1, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1], 
            RGBColor[1, 0, 1]}, 
        ImageSize -> {400, 300}, 
		Frame -> True, 
        FrameLabel -> {\*""\!\(Log\_10\)(k)"", "E(k,t)"}, 
        RotateLabel -> False, AspectRatio -> 2\/3, 
        PlotLabel -> "Evolution of energy spectra vs. time"];\)\)

To make things a bit easier, I'll define z as a function of k (also making the arguments exact)

Clear[z,  kmin,  kmax,  c1,  c2,  b];

z[k_]  =  Exp[-(1/2) (10 Log[10, k])^2];

eqns  =  
      Hold[{D[y[k, t], 
              t] == -c1*
                D[(k^(5/2))*((y[k, t])^(3/2)), k] + 
              c2*D[(k^(7/2))*((y[k, t])^(1/2))*
                    D[(y[k, t]), k], k], y[k, 0] == z[k], 
          y[kmin, t] == 0, y[kmax, t] == 0}]

The problem with the equation is two-fold.  First, an exponential grid is typically a disaster for finite difference PDE solvers.  Second, even when that is transformed away, there are still some nasty terms lingering.  By using some intuition about what might cancel some of these terms and experimenting with a few transformations, I was able to come up with one that led to a nice equation.

eqns1  =  
      eqns  /.  y[k_, 
            t_]  ->  ((((ϕ[Log[10, k],  t]))^2)/(k^3))

Now that the arguments inside of the differentiation have been replaced, the Hold can be released and we can replace k with a logarithmic variable, ν = (Subscript[log, 10]) k

eqns2  = ( 
          ReleaseHold[
              eqns1]  /.  ϕ[Log[10, k], t] ->    ϕ[ν, 
                t])  /.  k -> (10^ν)

Because of the transformation, the boundary conditions are expressed as squares.  Some needless complication in the solution can be removed by taking square roots.

bcs  =  Rest[eqns2]

bcs  =  Map[Sqrt,  bcs, {2}]

bcs  =  PowerExpand[bcs]

Before attempting a numerical solution, it makes sense to simplify the equations as much as possible.

de  =  
      FullSimplify[
        First[eqns2], {ν > 0,  ϕ[ν, t] > 0}]

de  =  Map[(# (10^(-11 ν))) &,  de]

Here are the parameters for the extended range.

kmin = 10. ^-4; kmax = 10. ^4; c1 = 1; c2 = 1;

With the transformation, and default settings, the solution can be computed in just over half a minute

Timing[
      sol  =  NDSolve[{de,  bcs}, ϕ, {ν,  Log[10. , kmin],  
            Log[10. ,  kmax]}, {t, 0, 100}]]

Tuning the settings of the numerical solver leads to a quicker, better,  solution.

\(Timing[
      sol\  = \ NDSolve[{de, \ bcs}, ϕ, {ν, \ Log[10. , kmin], \ 
            Log[10. , \ kmax]}, {t, 0, 100}, \ 
          Method -> {"MethodOfLines", \ "SpatialDiscretization" -> {"TensorProductGrid", \ "DifferenceOrder" -> 
                    6}}]]\)

This uses the transform to express the solution in terms of the transformed solution.

ψ  =  First[ϕ  /.  sol];
\(\(b[k_, \ t_]\  = \ ψ[Log[10, k], \ t]\^2\/k\^3;\)\)

This makes a plot of the spectra

\(\(LogLogPlot[{z[k], Re[b\ [k, 0.1]], Re[b\ [k, 0.2]], 
          Re[b\ [k, 0.7]], 
Re[b\ [k, 10]], Re[b\ [k, 100]]}\ , {k, kmin, 
          kmax}, 
		GridLines -> 
          None, 
		PlotRange -> {{10\^\(-4\), 10\^4}, {10\^\(-7\), 
              10\^1}}, 
        PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], 
            RGBColor[1, 1, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1], 
            RGBColor[1, 0, 1]}, 
        ImageSize -> {400, 300}, 
		Frame -> True, 
        FrameLabel -> {\*""\!\(Log\_10\)(k)"", "E(k,t)"}, 
        RotateLabel -> False, PlotPoints -> 10000, \ 
        AspectRatio -> 2\/3, 
        PlotLabel -> "Evolution of energy spectra vs. time"];\)\)

Last Slide

The discretization works essentially the same in multiple spatial dimensions as it does in one.

However, there are a few things to keep in mind:

This example comes from Michael Trott, who has been instrumental in pushing the PDE solving capabilities of Mathematica to the limit.

<< PSDirect`

Module[{ΔTC,ξ,η,𝒻,SEqs,packet,nds,ur,ui,makeColoredPolygon,makeAnimationFrame,
  ppξ = 320,ppη = 40,ΔT = 1/8,L = 5,TMax = 7},
  
ΔTC[f_][
              u_, {ξ_, η_}] := (D[
                    u, η, η] - η f[ξ] f''[ξ] D[
                      u, η] + η f'[ξ]^2 (2 D[
                            u, η] + η D[u, η, η]) - 
                  2 η f[ξ] f'[ξ] D[u, ξ, η] + 
                  f[ξ]^2 D[u, ξ, ξ])/f[ξ]^2;
  
𝒻 = 
            Function[ξ, (1 + 1/2 Cos[4  Pi ξ/(4  L)])];
  
SEqs = {-ΔTC[𝒻][ur[ξ, η, t], {ξ, η}] == -D[
                    ui[ξ, η, t], 
                    t], -ΔTC[𝒻][
                    ui[ξ, η, t], {ξ, η}] == +D[
                    ur[ξ, η, t], t]};
  
packet[{ξ_, η_}] := 
            If[Sqrt[(ξ + 2  L)^2 + η^2] < L/2, 
              Cos[Pi/2 Sqrt[(ξ + 2  L)^2 + η^2]/(L/2)]^2, 
              0];
  
state = First[NDSolve`ProcessEquations[
{N[SEqs], 
  
ur[ξ, η, 
                          0]  ==  
                        Cos[4 ξ] packet[{ξ, η}], ui[ξ, η, 0]  ==  
                        Sin[4 ξ] packet[{ξ, η}],
  
ur[-4 L, η, 
                          t] == ur[4  L, η, t],ui[-4 L, η, t] == 
                        ui[4  L, η, t],
 (* 
                          Boundary conditions at the waveguide edge *) 
Derivative[0, 1, 0][
                            ur][ξ, 0, t] == 0,ur[ξ, L, t] == 0,Derivative[0, 1, 0][ui][ξ, 0, 
                          t] == 0,ui[ξ, L, t] == 0},
{  },
{ξ, -4 L, 4  L},{η, 0, L},{t, 0, TMax},
MaxSteps -> ∞,PrecisionGoal -> 3,AccuracyGoal -> 3,Compiled -> True, \(Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> 
                            3, "MaxPoints" -> {ppξ, 
                              ppη}, "MinPoints" -> {ppξ, 
                              ppη}}}\)
]];
  
makeColoredPolygon[points_, maxAbs_] := 
            With[{exp = 1/3, meanAbs = Mean[Abs[Last /@ points]], 
                meanArg = If[# == 0. , 0, Arg[#]] &@
                    Mean[Last /@ points]}, {SurfaceColor[#, #, 2.2] &[
                  Hue[(meanArg + Pi)/(2  
                          Pi),   (meanAbs/maxAbs)^exp, (meanAbs/maxAbs)^exp]], 
                Polygon[{#1, #2, Abs[#3]^2} & @@@ points]}];
  
makeAnimationFrame[ndstate_]:=Block[{points, maxAbs, polys,  τ,  sol},
τ  =  
                    ndstate@"CurrentTime"["Forward"];
sol  =  
                    NDSolve`ProcessSolutions[ndstate,  "Forward"];
  
points = 
                    Table[Evaluate[{ξ, η 𝒻[ξ], 
                            ur[ξ, η, τ] + 
                              I ui[ξ, η, τ]} /. 
                          sol], {ξ, -4.  L, 4.  L, 
                        8.  L/400}, {η, 0. , 1.  L, 1.  L/50}];
maxAbs = 
                    Max[Abs[Map[Last, points, {2}]]];
  
polys = 
                    Table[makeColoredPolygon[{points[[i, j]], 
                          points[[i + 1, j]], points[[i + 1, j + 1]], 
                          points[[i, j + 1]]}, maxAbs], {i, 
                        Length[points] - 1}, {j, 
                        Length[points[[1]]] - 1}];
  
Show[
                    Graphics3D[{EdgeForm[], polys, 
                        polys /. 
                          Polygon[l_] :> 
                            Polygon[{1, -1, 1} # & /@ l]}, 
                      BoxRatios -> {4, 1, 1.3}, 
                      PlotRange -> {All, All, {0, 1.01}}, 
                      ViewPoint -> {1, -3, 2}, Boxed -> False, 
                      Axes -> False]]];
  
monitor  =  GUIRun[textmonitorwidget];
t1  =  t2  =  t3  =  
                0;
Catch[
Do[
  
timem  =  
                      StringJoin[" t = ",  
                        ToString[(k - 
                                1. ) ΔT],  "\n RAM in use \
= ", ToString[ .1  Round[MemoryInUse[]/
                                10^5]],  " MB \n Graphics time ",  
                        ToString[t2  -  t1],  "\n NDSolve time ",  
                        ToString[t3  -  t2]];
monitor @ 
                      SetPropertyValue[{"time",  "text"}, timem];
If[ 
                      TrueQ[ monitor @ Script[stopRequested]],  
                      monitor @ Script[stopRequested  =  False]; 
                      Throw[Null,  "StopIntegration"]];
t1  =  AbsoluteTime[];
makeAnimationFrame[state];
t2  =  AbsoluteTime[];
  
NDSolve`Iterate[state, 
                      k  ΔT];
t3  =  AbsoluteTime[],
{k, TMax/ΔT}], 
"StopIntegration"];
state]

Mathematica's symbolic capabilities make it possible to use the method of lines as a very general method for solving evolution equations.  The method attempts to choose parameters which work well in general, but considering the variety of solutions possible with PDEs, it is not surprising that these parameters can usually be tuned to obtain a solution more quickly or accurately.h