In[105]:=

This gives the number of spatial and temporal points used, respectively.

In[107]:=

The temporal points are chosen adaptively by the ODE method based on local error control. NDSolve used 49 (50 including the periodic endpoint) spatial points. This choice will be illustrated through the steps that follow.

In the equation processing phase of NDSolve, one of the first things that happen is that equations with second (or higher) order temporal derivatives are replaced with systems with only first order temporal derivatives.

This is a first order system equivalent to the sine-Gordon equation above.

In[108]:=

The next stage is to solve for the temporal derivatives.

This is the solution for the temporal derivatives, with the right hand side of the equations in normal (ODE) form.

In[109]:=

Now the problem is to choose a uniform grid that will approximate the derivative to within local error tolerances as specified by AccuracyGoal and PrecisionGoal. For this illustration, use the default DifferenceOrder (4) and the default AccuracyGoal and PrecisionGoal (both 4 for PDEs). The methods used to integrate the system of ODEs that results from discretization base their own error estimates on the assumption of sufficiently accurate function values. The estimates here have the goal of finding a spatial grid for which (at least with the initial condition) the spatial error is somewhat balanced with the local temporal error.

This sets variables to reflect the default settings for DifferenceOrder, AccuracyGoal, and PrecisionGoal.

In[110]:=

The error estimate is based on Richardson extrapolation. If you know that the error is and you have two approximations and at different values, and of h, then you can, in effect, extrapolate to the limit h->0 to get an error estimate

so the error in is estimated to be

Here and are vectors of different length and is a function, so we need to choose an appropriate norm. If we choose = 2 , then we can simply use a scaled norm on the components common to both vectors, which is all of and every other point of .This is a good choice because it does not require any interpolation between grids.

For a given interval on which we want to set up a uniform grid, we can define a function , where is the length of the interval such that the grid is , , ..., }, where .

This defines functions that return the step size and a list of grid points as a function of for the sine-Gordon equation.

In[113]:=

For a given grid, the equation can be discretized using finite differences. This is easily done using NDSolve`FiniteDifferenceDerivative.

This defines a symbolic discretization of the right hand side of the sine-Gordon equation as a function of a grid. It returns a function of u and v, which gives the approximate values for and in a list. (Note that in principle this works for any grid, uniform or not, though in the following, only uniform grids will be used.)

In[116]:=

For a given step size and grid, we can also discretize the initial conditions for u and v.

This defines a function that discretizes the initial conditions for u and v. The last grid point is dropped because, by periodic continuation, it is considered the same as the first.

In[117]:=

The quantity of interest is the approximation of the right hand side for a particular value of with this initial condition.

This defines a function that returns a vector consisting of the approximation of the right hand side of the equation for the initial condition for a given step size and grid. The vector is flattened to make subsequent analysis of it simpler.

In[118]:=

Starting with a particular value of , we can obtain the error estimate by generating the right hand side for and points.

This gives an example of the right hand side approximation vector for a grid with 10 points.

In[119]:=

This gives an example of the right hand side approximation vector for a grid with 20 points.

In[120]:=

As mentioned above, every other point on the grid with points lies on the grid with points. Thus, for simplicity, we can use a norm that only compares point common to both grids. Because the goal is to ultimately satisfy absolute and relative tolerance criteria, it is appropriate to use a scaled norm. In addition to taking into account the size of the right hand side for the scaling, it is also important to include the size of the corresponding components of u and v on the grid since error in the right hand side is ultimately included in u and v.

This defines a norm function for the difference in the approximation of the right hand side.

In[121]:=

This applies the norm function to the two approximations found above.

In[122]:=

To get the error estimate form the distance, according to the Richardson extrapolation formula (3), this simply needs to be divided by .

This computes the error estimate for n 10. Since this is based on a scaled norm, the tolerance criteria are satisfied if the result is less than 1.

In[123]:=

This makes a function that combines the above functions to give an error estimate as a function of n.

In[124]:=

The goal is to find the minimum value of n, such that the error estimate is less than or equal to 1 (since it is based on a scaled norm). In principle, it would be possible to use a root finding algorithm on this, but since n can only be an integer, this would be overkill and adjustments would have to be made to the stopping conditions. An easier solution is simply to use the simple Richardson extrapolation formula to predict what value of n would be appropriate and repeat the prediction process until the appropriate n is found.

The condition to satisfy is

and we have estimated that

so we can project that

or in terms of n, which is proportional to the reciprocal of h,

This computes the predicted optimal value of n based on the error estimate for n 10 computed above.

In[125]:=

This computes the error estimate for the new value of n.

In[126]:=

Often the case that a prediction based on a very coarse grid will be inadequate. A coarse grid may completely miss some solution features that contribute to the error on a finer grid. Also, the error estimate is based on an asymptotic formula, so for coarse spacings, the estimate itself may not be very good, even when all the solution features are resolved to some extent.

In practice, it can be fairly expensive to compute these error estimates. Also, it is not necessary to find the very optimal n, but one that satisfies the error estimate. Remember, everything can change as the PDE evolves, so it is simply not worth a lot of extra effort to find an optimal spacing for just the initial time. A simple solution is to include an extra factor greater than 1 in the prediction formula for the new n. Even with an extra factor, it may still take a few iterations to get to an acceptable value. It does, however, typically make the process faster.

This defines a function that gives a predicted value for the number of grid points, which should satisfy the error estimate.

In[127]:=

This iterates the predictions until a value is found that satisfies the error estimate.

In[129]:=

The prediction formula can be modified to use n instead of p in a similar way.

This modifies the function predicting an appropriate value of n to work with pseudospectral derivatives. This formulation does not try to pick an efficient FFT length.

In[130]:=

When finalizing the choice of n for a pseudospectral method, an additional consideration is to choose a value that not only satisfies the tolerance conditions, but is also an efficient length for computing FFTs. In Mathematica, an efficient FFT does not require a power of two length since the Fourier command has a prime factor algorithm built in.

Typically, the difference order has a profound effect on the number of points required to satisfy the error estimate.

This makes a table of the number of points required to satisfy the a-priori error estimate as a function of the difference order.

In[131]:=

A table of the number of points required as a function of difference order goes a long way towards explaining why the default setting for the method of lines is DifferenceOrder->4: the improvement from 2 to 4 is usually most dramatic and in the default tolerance range, fourth order differences do not tend to produce large roundoff errors, which can be the case with higher orders. Pseudospectral differences are often a good choice, particularly with periodic boundary conditions, but it is not a good choice for the default because they lead to full Jacobian matrices, which can be expensive to generate and solve if needed for stiff equations.

For non-periodic grids, the error estimate is done using only interior points. The reason is that the error coefficients for the derivatives near the boundary are different. This may miss features that are near the boundary, but the main idea is to keep the estimate simple and inexpensive since the evolution of the PDE may change everything anyway.

For multiple spatial dimensions, the determination is made one dimension at a time. Since better resolution in one dimension may change the requirements for another, the process is repeated in reverse order to improve the choice.

A-posteriori Error Estimates

When the solution of a PDE is computed with NDSolve, a final step is to do a spatial error estimate on the evolved solution and issue a warning message if this is excessively large.

These error estimates are done in a manner similar to the a-priori estimates described previously. The only real difference is that, instead of using grids with n and 2n points to estimate the error, grids with n/2 and n points are used. This is because, while there is no way to generate the values on a grid of 2n points without using interpolation, which would introduce its own errors, values are readily available on a grid of n/2 points simply by taking every other value. This is easily done in the Richardson extrapolation formula by using 2 , which gives

This defines a function (based on functions defined in the previous section) that can compute an error estimate on the solution of the sine-Gordon equation from solutions for u and v expressed as vectors. The function has been defined to be a function of the grid since this is applied to a grid already constructed. (Note, as defined here, this only works for grids of even length. It is not difficult to handle odd length, but it makes the function somewhat more complicated.)

In[132]:=

This solves the sine-Gordon equation with a Gaussian initial condition.

In[134]:=

This makes a function that gives the a-posteriori error estimate at a particular numerical value of t.

In[136]:=

This makes a plot of the a-posteriori error estimate as a function of t.

In[138]:=

This shows a plot of the solution at t 4. It is apparent that the warning message is appropriate because the oscillations near the peak are not physical.

All of the options for tensor product grid discretization can be given as a list with length equal to the number of spatial dimensions, in which case the parameter for each spatial dimension is determined by the corresponding component of the list.

With the exception of pseudospectral methods on non-periodic problems, discretization is done with uniform grids, so when solving a problem on interval length , there is a direct correspondence between the Points and StepSize options:

The StepSize options are effectively converted to the equivalent Points values. They are simply provided for convenience since sometimes it is more natural to relate problem specification to step size rather then the number of discretization points. When values other then Automatic are specified for both the Points and corresponding StepSize option, generally, the more stringent restriction is used.

In the previous section an example was shown where the solution was not resolved sufficiently because the solution steepened as it evolved. The examples below will show some different ways of modifying the grid parameters so that the near shock is better resolved.

One way to avoid the oscillations that showed up in the solution as the profile steepened is to make sure that we use sufficient points to resolve the profile at its steepest. In the one-hump solution of Burgers` equation,

it can be shown [W76] that the width of the shock profile is proportional to as ->0. More than 95% of the change is included in a layer of width 10 . Thus, if we pick a maximum step size of half the profile width, there will always be a point somewhere in the steep part of the profile, and there is a hope of resolving it without significant oscillation.

This computes the solution to Burgers' equation, such that there are sufficient points to resolve the shock profile.

In[140]:=

Note that resolving the profile alone is by no means sufficient to meet the default tolerances of NDSolve, which requires an accuracy of . However, once we have sufficient point to resolve the basic profile, it is not unreasonable to project based on the a-posteriori error estimate shown in the NDSolve::eerr message (with an extra 10% since, after all, it is just a projection).

This computes the solution to Burgers' equation with the maximum step size chosen so that it should be small enough to meet the default error tolerances based on a projection from the error of the previous calculation.

In[142]:=

To compare solutions like this, it is useful to look at a plot of the solution only at the spatial grid points. Because the grid points are stored as a part of the InterpolatingFunction, it is fairly simple to define a function that does this.

This defines a function that plots a solution only at the spatial grid points at a time t.

In[144]:=

This makes a plot comparing the three solutions found at t = 4.

In[146]:=

This makes a plot of the values of the solution at the assigned spatial grid points.

In[150]:=

This shows a surface plot of the leading edge of the solution a t = 2.

In[155]:=

It is important to keep in mind that the a-posteriori spatial error estimates are simply estimates of the local error in computing spatial derivatives and may not reflect the actual accumulated spatial error for a given solution. One way to get an estimate on the actual spatial error is to compute the solution to very stringent tolerances in time for different spatial grids. To show how this works, consider again the simpler 1-dimensional Burgers` equation.

This computes a list of solutions using {33, 65, ..., 4097} spatial grid points to compute the solution to Burgers` equation for difference orders 2, 4, 6, and pseudospectral. The temporal accuracy and precision tolerances are set very high so that essentially all of the error comes from the spatial discretization. Note that by specifying {t, 4, 4} in NDSolve, only the solution at t = 4 is saved. Without this precaution, some of the solutions for the finer grids (which take many more time steps) could exhaust available memory. Even so, the list of solutions takes a substantial amount of time to compute.

In[158]:=

Given two solutions, a comparison needs to be done between the two. To keep out any sources of error except for that in the solutions themselves, it is best to use the data that is interpolated to make the InterpolatingFunction. This can be done by using points common to the two solutions.

This defines a function to estimate error by comparing two different solutions at the points common to both. The arguments coarse and fine should be the solutions on the coarser and finer grids, respectively. This works for the solutions generated above with grid spacing varying by powers of 2.

In[160]:=

To get an indication of the general trend in error (in cases of instability, solutions do not converge, so this does not assume that), you can compare the difference of successive pairs of solutions.

This defines a function that will plot a sequence of error estimates for the successive solutions found for a given difference order and uses it to make a logarithmic plot of the estimated error as a function of the number of grid points.

The upper left part of the plot are the grids where the profile is not adequately resolved, so differences are simply of magnitude order 1 (it would be a lot worse if there was instability) However, once there are a sufficient number of points to resolve the profile without oscillation, convergence becomes quite rapid. Not surprisingly, the slope of the logarithmic line is -4, which corresponds to the difference order NDSolve uses by default. If your grid is fine enough to be in the asymptotically converging part, a simpler error estimate could be effected by using Richardson extrapolation as in the previous two sections, but on the overall solution rather than the local error. On the other hand, computing more values and viewing a plot gives a better indication of whether you are in the asymptotic regime or not.

It is fairly clear from the plot that the best solution computed is the pseudospectral one with 2049 points (the one with more points was not computed because its spatial accuracy far exceeds the temporal tolerances that were set). This solution can, in effect, be treated almost as an exact solution, at least up to error tolerances of or so.

To get a perspective of how best to solve the problem, it is useful to do the following: for each solution found that was at least a reasonable approximation, recompute it with the temporal accuracy tolerance set to be comparable to the possible spatial accuracy of the solution and plot the resulting accuracy as a function of solution time. The following (somewhat complicated) commands do this.

This identifies the "best" solution that will be used, in effect, as an exact solution in the computations that follow. It is dropped from the list of solutions to compare it to since the comparison would be meaningless.

In[165]:=

This defines a function that, given a difference order, do, and a solution, sol, computed with that difference order, recomputes it with local temporal tolerance slightly more stringent than the actual spatial accuracy achieved if that accuracy is sufficient. The function output is a list of {number of grid points, difference order, time to compute in seconds, actual error of the recomputed solution}.

In[167]:=

This applies the function to each of the previously computed solutions. (With the appropriate difference order!)

In[168]:=

This removes the cases that were not recomputed and makes a logarithmic plot of accuracy as a function of computation time.

The resulting graph demonstrates quite forcefully that, when they work, as in this case, periodic pseudospectral approximations are incredibly efficient. Otherwise, up to a point, the higher the difference order, the better the approximation will generally be. These are all features of smooth problems, which, this particular instance of Burgers` equation is. However, the higher order solutions would generally be quite poor if we went towards the limit ->0.

In[2]:=

The NDSolve:ibcinc message is issued, in this case, completely to the larger discretization error at the right boundary. For this particular example, the extra error is not a problem because it gets damped out due to the nature of the equation. However, it is possible to eliminate the message by using just a few more spatial points.

This computes the solution to the same equation as above, but using a minimum of 50 points in the x-direction.

In[8]:=

One other case where error problems at the boundary can affect the discretization unexpectedly is when periodic boundary conditions are given with a function which is not truly periodic, so that an unintended discontinuity is introduced into the computation.

This begins the computation of the solution to the sine-Gordon equation with a Gaussian initial condition and periodic boundary conditions. The NDSolve command is wrapped with TimeConstrained since solving the given problem can take a very long time and a large amount of system memory.

In[7]:=

The problem here is that the initial condition is effectively discontinuous when the periodic continuation is taken into account.

This shows a plot of the initial condition over the extent of three full periods.