Testing Methods: Solving Initial Value Problems¶
In addition to directly analyzing solver properties, NodePy also facilitates the testing of solvers through application to problems of interest. Furthermore, NodePy includes routines for automatically running sets of tests to compare the practical convergence or efficiency of various solvers for given problem(s).
Contents
Initial Value Problems¶
The principal objects in NodePy are ODE solvers. The object upon which a solver acts is an initial value problem. Mathematically, an initial value problem (IVP) consists of one or more ordinary differential equations and an initial condition:
\begin{align*} u’(t) & = F(u) & u(0) & = u_0. \end{align*}
In NodePy, an initial value problem is an object with the following properties:
rhs(): The right-hand-side function; i.e. F where \(u(t)'=F(u)\).
u0: The initial condition.
T: The (default) final time of solution.
- Optionally an IVP may possess the following:
exact(): a function that takes one argument (t) and returns the exact solution (Should we make this a function of u0 as well?)
dt0: The default initial timestep when a variable step size integrator is used.
Any other problem-specific parameters.
The module ivp contains functions for loading a variety of initial value problems. For instance, the van der Pol oscillator problem can be loaded as follows:
>> from NodePy import ivp
>> myivp = ivp.load_ivp('vdp')
Instantiation¶
-
ivp.detest()¶ Load problems from the non-stiff DETEST problem set. The set consists of six groups of problems, as follows:
A1-A5 – Scalar problems
B1-B5 – Small systems (2-3 equations)
C1-C5 – Moderate size systems (10-50 equations)
D1-D5 – Orbit equations with varying eccentricities
E1-E5 – Second order equations
F1-F5 – Problems with discontinuities
SB1-SB3 – Periodic Orbit problem from Shampine Baca paper pg.11,13
Note
Although this set of problems was not intended to become a standard, and although there are certain dangers in accepting any particular set of problems as a universal standard, it is nevertheless sometimes useful to try a new method on this test set due to the availability of published results for many existing methods.
Reference: [EP87].
Solving Initial Value Problems¶
Any ODE solver object in NodePy can be used to solve an initial value problem simply by calling the solver with an initial value problem object as argument:
>> t,u = rk44(my_ivp)
-
ODESolver.__call__(ivp, t0=0, N=5000, dt=None, errtol=None, controllertype='P', x=None, diagnostics=False, use_butcher=False, max_steps=7500)[source]¶ Calling an ODESolver numerically integrates the ODE u’(t) = f(t,u(t)) with initial value u(0)=u0 from time t0 up to time T using the solver.
The timestep can be controlled in any of three ways:
By specifying N, the total number of steps to take. Then dt = (T-t0)/N.
By specifying dt directly.
For methods with an error estimate (e.g., RK pairs), by specifying an error tolerance. Then the step size is adjusted using a PI-controller to achieve the requested tolerance. In this case, dt should also be specified and determines the value of the initial timestep.
The argument x is used to pass any additional arguments required for the RHS function f.
- Input:
ivp – An IVP instance (the initial value problem to be solved)
t0 – The initial time from which to integrate
N – The # of steps to take (using a fixed step size)
dt – The step size to use
errtol – The local error tolerance to be observed (using adaptive stepping). This requires that the method have an error estimator, such as an embedded Runge-Kutta method.
controllerType – The type of adaptive step size control to be used; available options are ‘P’ and ‘PI’. See [HNorW93] for details.
diagnostics – if True, return the number of rejected steps and a list of step sizes used, in addition to the solution values and times.
- Output:
t – A list of solution times
u – A list of solution values
If ivp.T is a scalar, the solution is integrated to that time and all output step values are returned. If ivp.T is a list/array, the solution is integrated to T[-1] and the solution is returned for the times specified in T.
TODO:
Implement an option to not keep all output (for efficiency).
Option to keep error estimate history
Convergence Testing¶
-
convergence.ctest(ivp, grids=[20, 40, 80, 160, 320, 640], verbosity=0, parallel=False)¶ Runs a convergence test, integrating a single initial value problem using a sequence of fixed step sizes and a set of methods. Creates a plot of the resulting errors versus step size for each method.
- Inputs:
methods – a list of ODEsolver instances
ivp – an IVP instance
- grids – a list of grid sizes for integration.
optional; defaults to [20,40,80,160,320,640]
parallel – to exploit possible parallelization (optional)
Example:
>>> import nodepy.runge_kutta_method as rk >>> from nodepy.ivp import load_ivp >>> rk44=rk.loadRKM('RK44') >>> myivp=load_ivp('nlsin') >>> work, err=ctest(rk44,myivp)
- TODO:
Option to plot versus f-evals or dt
Performance testing with automatic step-size control¶
-
convergence.ptest(ivps, tols=[0.1, 0.01, 0.0001, 1e-06], verbosity=0, parallel=False, controllertype='P')¶ Runs a performance test, integrating a set of problems with a set of methods using a sequence of error tolerances. Creates a plot of the error achieved versus the amount of work done (number of function evaluations) for each method.
- Input:
- methods – a list of ODEsolver instances
Note that all methods must have error estimators.
ivps – a list of IVP instances
tols – a specified list of error tolerances (optional)
parallel – to exploit possible parallelization (optional)
Example:
>>> import nodepy.runge_kutta_method as rk >>> from nodepy.ivp import load_ivp >>> bs5=rk.loadRKM('BS5') >>> myivp=load_ivp('nlsin') >>> work,err=ptest(bs5,myivp)