Numerical Solution of Partial Differential Equations
The Numerical Method of Lines
Introduction
The numerical method of lines is a technique for solving partial differential equations discretizing in all but one dimension, and then integrating the semi-discrete problem as a system of ODEs or DAEs. A significant advantage of the method is that it allows the solution to take advantage of the sophisticated general-purpose methods and software that have been developed for numerically integrating ODEs and DAEs. For the PDEs to which the method of lines is applicable, the method typically proves to be quite efficient.
It is necessary that the PDE problem be well-posed as an initial value (Cauchy) problem in at least one dimension, since the ODE and DAE integrators used are initial value problem solvers. This rules out purely elliptic equations such as Laplace's equation, but leaves a large class of evolution equations that can be solved quite efficiently.
A simple example illustrates better than mere words the fundamental idea of the method. Consider the following problem (a simple model for seasonal variation of heat in soil).
This is a candidate for the method of lines since you have the initial value
u (x, 0)
0.
Problem (1) will be discretized with respect to the variable
x using second-order finite differences, in particular using the approximation
Even though finite difference discretizations are the most common, there is certainly no requirement that discretizations for the method of lines be done with finite differences; finite volume or even finite element discretizations can also be used.
To use the discretization shown, choose a uniform grid
xi, 0≤i≤n with spacing
h1/n such that
xiih. Let
ui[t] be the value of
u (xi, t). For the purposes of illustrating the problem setup, a particular value of
n is chosen.
This defines a particular value of n and the corresponding value of h used in the subsequent commands. This can be changed to make a finer or coarser spatial approximation. |
This defines the vector of ui.
| Out[2]= |  |
|
For
1≤i≤9, you can use the centered difference formula (2) to obtain a system of ODEs. However, before doing this, it is useful to incorporate the boundary conditions first.
The Dirichlet boundary condition at
x0 can easily be handled by simply defining
u0as a function of
t. An alternative option is to differentiate the boundary condition with respect to time and use the corresponding differential equation. In this example, the latter method will be used.
The Neumann boundary condition at
x1 is a little more difficult. With second-order differences, one way to handle it is with reflection: imagine that you are solving the problem on the interval
0≤x≤2 with the same boundary conditions at
x0 and
x2. Since the initial condition and boundary conditions are symmetric with respect to
x, the solution should be symmetric with respect to
x for all time, and so symmetry is equivalent to the Neumann boundary condition at
x
1. Thus,
u (1+h, t)=u (1-h, t), so
un+1[t]=un-1[t].
This uses ListCorrelate to apply the difference formula. The padding {un-1[t]} implements the Neumann boundary condition.
| Out[3]= |  |
|
This sets up the zero initial condition.
| Out[4]= |  |
|
Now the PDE has been partially discretized into an ODE initial value problem that can be solved by the ODE integrators in
NDSolve.
This solves the ODE initial value problem.
| Out[5]= |  |
|
This shows the solutions u (xi, t) plotted as a function of x and t.
| Out[6]= |  |
|
The plot indicates why this technique is called the numerical "method of lines".
The solution in between lines can be found by interpolation. When
NDSolve computes the solution for the PDE, the result is a two-dimensional
InterpolatingFunction.
This uses NDSolve to compute the solution of the heat equation (1) directly.
| Out[7]= |  |
|
This creates a surface plot of the solution.
| Out[8]= |  |
|
The setting
n10 used did not give a very accurate solution. When
NDSolve computes the solution, it uses spatial error estimates on the initial condition to determine what the grid spacing should be. The error in the temporal (or at least time-like) variable is handled by the adaptive ODE integrator.
In the example (1), the distinction between time and space was quite clear from the problem context. Even when the distinction is not explicit, this document will refer to "spatial" and "temporal" variables. The "spatial" variables are those to which the discretization is done. The "temporal" variable is the one left in the ODE system to be integrated.
| | |
| TemporalVariable | Automatic | what variable to keep derivatives with respect to the derived ODE or DAE system |
| Method | Automatic | what method to use for integrating the ODEs or DAEs |
| SpatialDiscretization | TensorProductGrid | what method to use for spatial discretization |
| DifferentiateBoundaryConditions | True | whether to differentiate the boundary conditions with respect to the temporal variable |
| ExpandFunctionSymbolically | False | whether to expand the effective function symbolically or not |
| DiscretizedMonitorVariables | False | whether to interpret dependent variables given in monitors like StepMonitor or in method options for methods like EventLocator and Projection as functions of the spatial variables or vectors representing the spatially discretized values |
Options for NDSolve`MethodOfLines.
Use of some of these options requires further knowledge of how the method of lines works and will be explained in the sections that follow.
Currently, the only method implemented for spatial discretization is the
TensorProductGrid method, which uses discretization methods for one spatial dimension and uses an outer tensor product to derive methods for multiple spatial dimensions on rectangular regions.
TensorProductGrid has its own set of options that you can use to control the grid selection process. The following sections give sufficient background information so that you will be able to use these options if necessary.
Spatial Derivative Approximations
Finite Differences
The essence of the concept of finite differences is embodied in the standard definition of the derivative
where instead of passing to the limit as
h approaches zero, the finite spacing to the next adjacent point,
xi+1=xi+h, is used so that you get an approximation.
The difference formula can also be derived from Taylor's formula,
which is more useful since it provides an error estimate (assuming sufficient smoothness)
An important aspect of this formula is that
i must lie between
xi and
xi+1 so that the error is local to the interval enclosing the sampling points. It is generally true for finite difference formulas that the error is local to the stencil, or set of sample points. Typically, for convergence and other analysis, the error is expressed in asymptotic form:
This formula is most commonly referred to as the first-order forward difference. The backward difference would use
xi-1.
Taylor's formula can easily be used to derive higher-order approximations. For example, subtracting
and solving for
f' (xi) gives the second-order centered difference formula for the first derivative,
If the Taylor formulas shown are expanded out one order farther and added and then combined with the formula just given, it is not difficult to derive a centered formula for the second derivative.
Note that the while having a uniform step size
h between points makes it convenient to write out the formulas, it is certainly not a requirement. For example, the approximation to the second derivative is in general
where
h corresponds to the maximum local grid spacing. Note that the asymptotic order of the three-point formula has dropped to first order; that it was second order on a uniform grid is due to fortuitous cancellations.
In general, formulas for any given derivative with asymptotic error of any chosen order can be derived from the Taylor formulas as long as a sufficient number of sample points are used. However, this method becomes cumbersome and inefficient beyond the simple examples shown. An alternate formulation is based on polynomial interpolation: since the Taylor formulas are exact (no error term) for polynomials of sufficiently low order, so are the finite difference formulas. It is not difficult to show that the finite difference formulas are equivalent to the derivatives of interpolating polynomials. For example, a simple way of deriving the formula just shown for the second derivative is to interpolate a quadratic and find its second derivative (which is essentially just the leading coefficient).
This finds the three-point finite difference formula for the second derivative by differentiating the polynomial interpolating the three points (xi-1, f (xi-1)), (xi, f (xi)), and (xi+1, f (xi+1)).
| Out[9]= |  |
|
In this form of the formula, it is easy to see that it is effectively a difference of the forward and backward first-order derivative approximations. Sometimes it is advantageous to use finite differences in this way, particularly for terms with coefficients inside of derivatives, such as
(a (x)ux)x, which commonly appear in PDEs.
Another property made apparent by considering interpolation formulas is that the point at which you get the derivative approximation need not be on the grid. A common use of this is with staggered grids where the derivative may be wanted at the midpoints between grid points.
This generates a fourth-order approximation for the first derivative on a uniform staggered grid, xi, where the main grid points xi+k/2 are at xi+h k/2, for odd k.
| Out[10]= |  |
|
The fourth-order error coefficient for this formula is
versus
for the standard fourth-order formula derived next. Much of the reduced error can be attributed to the reduced stencil size.
This generates a fourth-order approximation for the first derivative at a point on a uniform grid.
| Out[11]= |  |
|
In general, a finite difference formula using n points will be exact for functions that are polynomials of degree
n-1 and have asymptotic order at least
n-m. On uniform grids, you can expect higher asymptotic order, especially for centered differences.
Using efficient polynomial interpolation techniques is a reasonable way to generate coefficients, but B. Fornberg has developed a fast algorithm for finite difference weight generation [
F92], [
F98], which is substantially faster.
In [
F98], Fornberg presents a one-line
Mathematica formula for explicit finite differences.
This is the simple formula of Fornberg for generating weights on a uniform grid. Here it has been modified slightly by making it a function definition. |
Here
m is the order of the derivative,
n is the number of grid intervals enclosed in the stencil, and
s is the number of grid intervals between the point at which the derivative is approximated and the leftmost edge of the stencil. There is no requirement that
s be an integer; noninteger values simply lead to staggered grid approximations. Setting
s to be
n/2 always generates a centered formula.
This uses the Fornberg formula to generate the weights for a staggered fourth-order approximation to the first derivative. This is the same one computed earlier with InterpolatingPolynomial.
| Out[13]= |  |
|
A table of some commonly used finite difference formulas follows for reference.
Finite difference formulas on uniform grids for the first derivative.
Finite difference formulas on uniform grids for the second derivative.
One thing to notice from this table is that the farther the formulas get from centered, the larger the error term coefficient, sometimes by factors of hundreds. For this reason, sometimes where one-sided derivative formulas are required (such as at boundaries), formulas of higher order are used to offset the extra error.
NDSolve`FiniteDifferenceDerivative
Fornberg [
F92], [
F98] also gives an algorithm that, though not quite so elegant and simple, is more general and, in particular, is applicable to nonuniform grids. It is not difficult to program in
Mathematica, but to make it as efficient as possible, a new kernel function has been provided as a simpler interface (along with some additional features).
| NDSolve`FiniteDifferenceDerivative[Derivative[m],grid,values] |
| approximate the mth-order derivative for the function that takes on values on the grid |
| NDSolve`FiniteDifferenceDerivative[Derivative[m1,m2,...,mn],{grid1,grid2,...,gridn},values] |
| approximate the partial derivative of order (m1, m2, ..., mn) for the function of n variables which takes on values on the tensor product grid defined by the outer product of (grid1, grid2, ..., gridn) |
| NDSolve`FiniteDifferenceDerivative[Derivative[m1,m2,...,mn],{grid1,grid2,...,gridn}] |
| compute the finite difference weights needed to approximate the partial derivative of order (m1, m2, ..., mn) for the function of n variables on the tensor product grid defined by the outer product of (grid1, grid2, ..., gridn); the result is returned as an NDSolve`FiniteDifferenceDerivativeFunction, which can be repeatedly applied to values on the grid |
Finding finite difference approximations to derivatives.
This defines a uniform grid with points spaced apart by a symbolic distance h.
| Out[14]= |  |
|
This gives the first derivative formulas on the grid for a symbolic function f.
| Out[15]= |  |
|
The derivatives at the endpoints are computed using one-sided formulas. The formulas shown in the previous example are fourth-order accurate, which is the default. In general, when you use symbolic grid and/or data, you will get symbolic formulas. This is often useful for doing analysis on the methods; however, for actual numerical grids, it is usually faster and more accurate to give the numerical grid to
NDSolve`FiniteDifferenceDerivative rather than using the symbolic formulas.
This defines a randomly spaced grid between 0 and 2  .
| Out[16]= |  |
|
This approximates the derivative of the sine function at each point on the grid.
| Out[17]= |  |
|
This shows the error in the approximations.
| Out[18]= |  |
|
In multiple dimensions,
NDSolve`FiniteDifferenceDerivative works on tensor product grids, and you only need to specify the grid points for each dimension.
This defines grids xgrid and ygrid for the x and y direction, gives an approximation for the mixed xy partial derivative of the Gaussian on the tensor product of xgrid and ygrid, and makes a surface plot of the error.
| Out[23]= |  |
|
Note that the values need to be given in a matrix corresponding to the outer product of the grid coordinates.
NDSolve`FiniteDifferenceDerivative does not compute weights for sums of derivatives. This means that for common operators like the Laplacian, you need to combine two approximations.
This makes a function that approximates the Laplacian operator on a tensor product grid. |
This uses the function to approximate the Laplacian for the same grid and Gaussian function used in the previous example.
| Out[25]= |  |
|
| | |
| DifferenceOrder | 4 | asymptotic order of the error |
| PeriodicInterpolation | False | whether to consider the values as those of a periodic function with the period equal to the interval enclosed by the grid |
Options for NDSolve`FiniteDifferenceDerivative.
This approximates the derivatives for the sine function on the random grid defined earlier, assuming that the function repeats periodically.
| Out[26]= |  |
|
When using
PeriodicInterpolation->True, you can omit the last point in the values since it should always be the same as the first. This feature is useful when solving a PDE with periodic boundary conditions.
This generates second-order finite difference formulas for the first derivative of a symbolic function.
| Out[27]= |  |
|
Fourth-order differences typically provide a good balance between truncation (approximation) error and roundoff error for machine precision. However, there are some applications where fourth-order differences produce excessive oscillation (Gibb's phenomena), so second-order differences are better. Also, for high-precision, higher-order differences may be appropriate. Even values of
DifferenceOrder use centered formulas, which typically have smaller error coefficients than noncentered formulas, so even values are recommended when appropriate.
NDSolve`FiniteDifferenceDerivativeFunction
When computing the solution to a PDE, it is common to repeatedly apply the same finite difference approximation to different values on the same grid. A significant savings can be made by storing the necessary weight computations and applying them to the changing data. When you omit the (third) argument with function values in
NDSolve`FiniteDifferenceDerivative, the result will be an
NDSolve`FiniteDifferenceDerivativeFunction, which is a data object that stores the weight computations in a efficient form for future repeated use.
| NDSolve`FiniteDifferenceDerivative[{m1,m2,...},{grid1,grid2,...}] |
| compute the finite difference weights needed to approximate the partial derivative of order (m1, m2, ...) for the function of n variables on the tensor product grid defined by the outer product of (grid1, grid2, ...); the result is returned as an NDSolve`FiniteDifferenceDerivativeFunction object |
| NDSolve`FiniteDifferenceDerivativeFunction[Derivative[m],data] |
| a data object that contains the weights and other data needed to quickly approximate the mth-order derivative of a function; in the standard output form, only the Derivative[m] operator it approximates is shown |
| NDSolve`FiniteDifferenceDerivativeFunction[data][values] |
| approximate the derivative of the function that takes on values on the grid used to determine data |
Computing finite difference weights for repeated use.
This defines a uniform grid with 25 points on the unit interval and evaluates the sine function with one period on the grid.
| Out[4]= |  |
|
This defines an NDSolve`FiniteDifferenceDerivativeFunction, which can be repeatedly applied to different values on the grid to approximate the second derivative.
| Out[5]= |  |
|
Note that the standard output form is abbreviated and only shows the derivative operators that are approximated.
This computes the approximation to the second derivative of the sine function.
| Out[6]= |  |
|
This function is only applicable for values defined on the particular grid used to construct it. If your problem requires changing the grid, you will need to use
NDSolve`FiniteDifferenceDerivative to generate weights each time the grid changes. However, when you can use
NDSolve`FiniteDifferenceDerivativeFunction objects, evaluation will be substantially faster.
This compares timings for computing the Laplacian with the function just defined and with the definition of the previous section. A loop is used to repeat the calculation in each case because it is too fast for the differences to show up with Timing.
| Out[10]= |  |
|
An
NDSolve`FiniteDifferenceDerivativeFunction can be used repeatedly in many situations. As a simple example, consider a collocation method for solving the boundary value problem
on the unit interval. (This simple method is not necessarily the best way to solve this particular problem, but it is useful as an example.)
This defines a function that will have all components zero at an approximate solution of the boundary value problem. Using the intermediate vector v and setting its endpoints (parts {1,-1}) to 0 is a fast and simple trick to enforce the boundary conditions. Evaluation is prevented except for numbers  because this would not work otherwise. (Also, because Times is Listable, Sin[2 Pi grid] u would thread componentwise.) |
This uses FindRoot to find an approximate eigenfunction using the constant coefficient case for a starting value and shows a plot of the eigenfunction.
| Out[14]= |  |
|
With the setup for this problem is so simple, it is easy to compare various alternatives. For example, to compare the solution above, which used the default fourth-order differences, to the usual use of second-order differences, all that needs to be changed is the
DifferenceOrder.
This solves the boundary value problem using second-order differences and shows a plot of the difference between it and the fourth-order solution.
| Out[41]= |  |
|
One way to determine which is the better solution is to study the convergence as the grid is refined. This will be considered to some extent in the section on differentiation matrices below.
While the most vital information about an
NDSolve`FiniteDifferenceDerivativeFunction object, the derivative order, is displayed in its output form, sometimes it is useful to extract this and other information from an
NDSolve`FiniteDifferenceDerivativeFunction, say for use in a program. The structure of the way the data is stored may change between versions of
Mathematica, so extracting the information by using parts of the expression is not recommended. A better alternative is to use any of the several method functions provided for this purpose.
Let
FDDF represent an
NDSolve`FiniteDifferenceDerivativeFunction[data] object.
| FDDF@"DerivativeOrder" | get the derivative order which FDDF approximates |
| FDDF@"DifferenceOrder" | get the list with the difference order used for the approximation in each dimension |
| FDDF@"PeriodicInterpolation" | get the list with elements True or False indicating whether periodic interpolation is used for each dimension |
| FDDF@"Coordinates" | get the list with the grid coordinates in each dimension |
| FDDF@"Grid" | form the tensor of the grid points; this is the outer product of the grid coordinates |
| FDDF@"DifferentiationMatrix" | compute the sparse differentiation matrix mat such that mat.Flatten[values] is equivalent to Flatten[FDDF[values]] |
Method functions for exacting information from an NDSolve`FiniteDifferenceDerivativeFunction[data] object.
Any of the method functions that return a list with an element for each of the dimensions can be used with an integer argument
dim, which will return only the value for that particular dimension such that
FDDF@method[dim]=(FDDF@method)[[dim]].
The following examples show how you might use some of these methods.
Here is an NDSolve`FiniteDifferenceDerivativeFunction object created with random grids having between 10 and 16 points in each dimension.
| Out[15]= |  |
|
This shows the dimensions of the outer product grid.
| Out[20]= |  |
|
Note that the rank of the grid point tensor is one more than the dimensionality of the tensor product. This is because the three coordinates defining each point are in a list themselves. If you have a function that depends on the grid variables, you can use
Apply[f, fddf["Grid"], {n}] where
n=Length[fddf["DerivativeOrder"]] is the dimensionality of the space in which you are approximating the derivative.
This defines a Gaussian function of 3 variables and applies it to the grid on which the NDSolve`FiniteDifferenceDerivativeFunction is defined. |
This shows a 3-dimensional plot of the grid points colored according to the scaled value of the derivative.
| Out[23]= |  |
|
For a moderate-sized tensor product grid like the example here, using
Apply is reasonably fast. However, as the grid size gets larger, this approach may not be the fastest because
Apply can only be used in limited ways with the
Mathematica compiler and hence, with packed arrays. If you can define your function so you can use
Map instead of
Apply, you may be able to use a
CompiledFunction since
Map has greater applicability within the
Mathematica compiler than does
Apply.
This defines a CompiledFunction which uses Map to get the values on the grid. (If the first grid dimension is greater than the system option "MapCompileLength", then you do not need to construct the CompiledFunction since the compilation is done automatically when grid is a packed array.)
| Out[24]= |  |
|
An even better approach, when possible, is to take advantage of listability when your function consists of operations and functions which have the
Listable attribute. The trick is to separate the
x,
y, and
z values at each of the points on the tensor product grid. The fastest way to do this is using
Transpose[fddf["Grid"]], RotateLeft[Range[n+1]]], where
n=Length[fddf["DerivativeOrder"]] is the dimensionality of the space in which you are approximating the derivative. This will return a list of length
n, which has the values on the grid for each of the component dimensions separately. With the
Listable attribute, functions applied to this will thread over the grid.
This defines a function that takes advantage of the fact that Exp has the Listable attribute to find the values on the grid. |
This compares timings for the three methods. The commands are repeated several times to get more accurate timings.
| Out[26]= |  |
|
The example timings show that using the
CompiledFunction is typically much faster than using
Apply and taking advantage of listability is a little faster yet.
Pseudospectral Derivatives
The maximum value the difference order can take on is determined by the number of points in the grid. If you exceed this, a warning message will be given and the order reduced automatically.
This uses maximal order to approximate the first derivative of the sine function on a random grid.
| Out[50]= |  |
|
Using a limiting order is commonly referred to as a
pseudospectral derivative. A common problem with these is that artificial oscillations (Runge's phenomena) can be extreme. However, there are two instances where this is not the case: a uniform grid with periodic repetition and a grid with points at the zeros of the Chebyshev polynomials,
Tn, or Chebyshev-Gauss-Lobatto points [
F96a], [
QV94]. The computation in both of these cases can be done using a fast Fourier transform, which is efficient and minimizes roundoff error.
| DifferenceOrder->n | use nth-order finite differences to approximate the derivative |
| DifferenceOrder->Length[grid] | use the highest possible order finite differences to approximate derivative on grid (not generally recommended.) |
| DifferenceOrder->"Pseudospectral" | use a pseudospectral derivative approximation; only applicable when the grid points are spaced corresponding to the Chebyshev-Gauss-Lobatto points or when the grid is uniform with PeriodicInterpolation ->True |
| DifferenceOrder->{n1,n2,...} | use difference orders n1, n2, ... in dimensions 1, 2, ... respectively |
Settings for the DifferenceOrder option.
This gives a pseudospectral approximation for the second derivative of the sine function on a uniform grid.
| Out[28]= |  |
|
This computes the error at each point. The approximation is accurate to roundoff because the effective basis for the pseudospectral derivative on a uniform grid for a periodic function are the trigonometric functions.
| Out[29]= |  |
|
The Chebyshev-Gauss-Lobatto points are the zeros of
(1 - x2)Tn' (x). Using the property
Tn (x) = Tn (cos (
))cos (n ), these can be shown to be at
.
This defines a simple function that generates a grid of n points with leftmost point at x0 and interval length L having the spacing of the Chebyshev-Gauss-Lobatto points |
This computes the pseudospectral derivative for a Gaussian function.
| Out[31]= |  |
|
This shows a plot of the approximation and the exact values.
| Out[32]= |  |
|
This shows a plot of the derivative computed using a uniform grid with the same number of points with maximal difference order.
| Out[36]= |  |
|
Even though the approximation is somewhat better in the center (because the points are more closely spaced there in the uniform grid), the plot clearly shows the disastrous oscillation typical of overly high-order finite difference approximations. Using the Chebyshev-Gauss-Lobatto spacing has minimized this.
This shows a plot of the pseudospectral derivative approximation computed using a uniform grid with periodic repetition.
| Out[71]= |  |
|
With the assumption of periodicity, the approximation is significantly improved. The accuracy of the periodic pseudospectral approximations is sufficiently high to justify, in some cases, using a larger computational domain to simulate periodicity, say for a pulse like the example. Despite the great accuracy of these approximations, they are not without pitfalls: one of the worst is probably aliasing error, whereby an oscillatory function component with too great a frequency can be misapproximated or disappear entirely.
Accuracy and Convergence of Finite Difference Approximations
When using finite differences, it is important to keep in mind that the truncation error, or the asymptotic approximation error induced by cutting off the Taylor series approximation, is not the only source of error. There are two other sources of error in applying finite difference formulas; condition error and roundoff error [
GMW81]. Roundoff error comes from roundoff in the arithmetic computations required. Condition error comes from magnification of any errors in the function values, typically from the division by a power of the step size, and so grows with decreasing step size. This means that in practice, even though the truncation error approaches zero as
h does, the actual error will start growing beyond some point. The following figures demonstrate the typical behavior as
h becomes small for a smooth function.
A logarithmic plot of the maximum error for approximating the first derivative of the Gaussian f (x)=
- (15 (x-1/2))2 at points on a grid covering the interval [0, 1] as a function of the number of grid points, n, using machine precision. Finite differences of order 2, 4, 6, and 8 on a uniform grid are shown in red, green, blue, and magenta respectively. Pseudospectral derivatives with uniform (periodic) and Chebyshev spacing are shown in black and gray respectively.
A logarithmic plot of the truncation error (dotted) and the condition and roundoff error (solid line) for approximating the first derivative of the Gaussian f (x)=- (15 (x-1/2))2 at points on a grid covering the interval [0, 1] as a function of the number of grid points, n. Finite differences of order 2, 4, 6, and 8 on a uniform grid are shown in red, green, blue, and magenta, respectively. Pseudospectral derivatives with uniform (periodic) and Chebyshev spacing are shown in black and gray, respectively. The truncation error was computed by computing the approximations with very high precision. The roundoff and condition error was estimated by subtracting the machine-precision approximation from the high-precision approximation. The roundoff and condition error tends to increase linearly (because of the 1/h factor common to finite difference formulas for the first derivative) and tends to be a little bit higher for higher-order derivatives. The pseudospectral derivatives show more variations because the error of the FFT computations vary with length. Note that the truncation error for the uniform (periodic) pseudospectral derivative does not decrease below about 10-22. This is because, mathematically, the Gaussian is not a periodic function; this error in essence gives the deviation from periodicity.
A semilogarithmic plot of the error for approximating the first derivative of the Gaussian f (x)=- (x-1/2)2 as a function of x at points on a 45 point grid covering the interval [0, 1]. Finite differences of order 2, 4, 6, and 8 on a uniform grid are shown in red, green, blue, and magenta, respectively. Pseudospectral derivatives with uniform (periodic) and Chebyshev spacing are shown in black and gray, respectively. All but the pseudospectral derivative with Chebyshev spacing were computed using uniform spacing 1/45. It is apparent that the error for the pseudospectral derivatives is not so localized; not surprising since the approximation at any point is based on the values over the whole grid. The error for the finite difference approximations are localized and the magnitude of the errors follows the size of the Gaussian (which is parabolic on a semilogarithmic plot).
From the second plot, it is apparent that there is a size for which the best possible derivative approximation is found; for larger
h, the truncation error dominates, and for smaller
h, the condition and roundoff error dominate. The optimal
h tends to give better approximations for higher-order differences. This is not typically an issue for spatial discretization of PDEs because computing to that level of accuracy would be prohibitively expensive. However, this error balance is a vitally important issue when using low-order differences to approximate, for example, Jacobian matrices. To avoid extra function evaluations, first-order forward differences are usually used, and the error balance is proportional to the square root of unit roundoff, so picking a good value of
h is important [
GMW81].
The plots showed the situation typical for smooth functions where there were no real boundary effects. If the parameter in the Gaussian is changed so the function is flatter, boundary effects begin to appear.
A semilogarithmic plot of the error for approximating the first derivative of the Gaussian f (x)=- (15 (x-1/2))2 as a function of x at points on a 45-point grid covering the interval [0, 1]. Finite differences of order 2, 4, 6, and 8 on a uniform grid are shown in red, green, blue, and magenta, respectively. Pseudospectral derivatives with uniform (nonperiodic) and Chebyshev spacing are shown in black and gray, respectively. All but the pseudospectral derivative with Chebyshev spacing were computed using uniform spacing 1/45. The error for the finite difference approximations are localized, and the magnitude of the errors follows the magnitude of the first derivative of the Gaussian. The error near the boundary for the uniform spacing pseudospectral (order-45
