ConvergenceTesting_1.gif > ConvergenceTesting_2.gif

SIMULATIONTOOLS TUTORIAL ConvergenceTesting_3.gif

Convergence Testing

Many simulations involve the solution of partial differential equations using finite difference methods on a coordinate grid.  Several simulations, identical apart from the grid spacing, are then compared to see how the numerical solution changes with the grid spacing.  This is done both to show that the solution converges to an exact solution, and to assess the numerical error.

SimulationTools contains several functions to help in this process.

In the examples below, we will use a set of three DataTables constructed to simulate a second-order accurate numerical solution to some partial differential equation.

The grid spacings are

In[1]:=

ConvergenceTesting_4.gif

The “numerical solutions” are

In[2]:=

ConvergenceTesting_5.gif

Out[2]=

ConvergenceTesting_6.gif

In[3]:=

ConvergenceTesting_7.gif

Out[3]=

ConvergenceTesting_8.gif

In[4]:=

ConvergenceTesting_9.gif

Out[4]=

ConvergenceTesting_10.gif

We can plot the DataTables directly:

In[5]:=

ConvergenceTesting_11.gif

Out[5]=

ConvergenceTesting_12.gif

Ratio between solution differences

Assuming convergence of the numerical solutions at a given order, the ratio between the differences of numerical solutions, ConvergenceTesting_13.gif and ConvergenceTesting_14.gif is given by

    ConvergenceTesting_15.gif

This is accurate to order p+1.  

In[6]:=

ConvergenceTesting_16.gif

Out[6]=

ConvergenceTesting_17.gif

For the example above, we can confirm that the data satisfies this property by plotting the two differences, with the second one rescaled by the convergence multiplier:

In[7]:=

ConvergenceTesting_18.gif

Out[7]=

ConvergenceTesting_19.gif

In[8]:=

ConvergenceTesting_20.gif

Out[8]=

ConvergenceTesting_21.gif

The two curves lie on top of each other, since the convergence is exactly 2nd order in this case.

The WithResampling function is necessary in order for subtraction of DataTables on different coordinate grids to use interpolation automatically (the interpolation is 8th order by default).

Convergence rate

You can also measure the convergence rate:

In[9]:=

ConvergenceTesting_22.gif

Out[9]=

ConvergenceTesting_23.gif

Note that the convergence rate cannot be determined at the points where ConvergenceTesting_24.gif

For general ConvergenceTesting_25.gif, ConvergenceTesting_26.gif and ConvergenceTesting_27.gif, the convergence rate cannot be solved for exactly; Mathematica’ FindRoot function is used to solve the algebraic equation.  As such, in certain cases a solution might not be found.

RichardsonExtrapolant

Assuming a given order of convergence, you can use two of the solutions to estimate a higher-order approximation of the exact solution.  This is called Richardson Extrapolation and is useful for providing an error estimate in one of the numerical solutions.  The RichardsonExtrapolant function gives this estimate:

In[10]:=

ConvergenceTesting_28.gif

Out[10]=

ConvergenceTesting_29.gif

An error estimate (accurate to 3rd order) for f3 can be determined by subtracting the Richardson extrapolant:

In[11]:=

ConvergenceTesting_30.gif

Out[11]=

ConvergenceTesting_31.gif

In[12]:=

ConvergenceTesting_32.gif

Out[12]=

ConvergenceTesting_33.gif

Spikey Created with Wolfram Mathematica 9.0