The initial value problem hairer1993solving ; hairer1996solving is very often an important sub-problem in many different applications. The aim of this project is to provide not only a useful collection of initial value problems, but an easy-to-use interface for working with them. Unlike other projects of this type, our aim is to have setup and use be as simple as possible, with almost no boilerplate code, but with the extendability that comes with more sophisticated, in-depth tools. The project was born out of necessity: disjoint collections of problems with their own setup scripts and uses existed internally. Their standardization in a unified framework written in MATLAB is simply the natural evolution, and a sign of things to come. All problems are implemented purely in MATLAB without any external dependencies or mex files.
Using OTP is as simple as
A non-trivial example showcasing the versatility of the framework is
This showcases that not only is it trivial to use complex models, but that modifying fundamental properties of the model is as easy as typing a few lines of code.
We consider the classic initial value problem
where . With all possible trajectories defining the space .
The three things that uniquely define an initial value problem are the function , commonly referred to as the ‘right hand side’ or RHS for short, the time span of the problem (from an initial to a final ), and the initial value, , defined to be the known value of the solution at the initial time.
Additional things that are helpful but not necessary in the solution of an initial value problem numerically might be the knowledge of higher order derivatives of , such as the Jacobian, spatial relations between the variables in terms of a distance function, splittings of the RHS into physical processes that occur at different temporal scales, etc.
Iii ODE Test Problems API
By far, the most useful part of OTP are presets. A preset provides problem parameters, the time span, and initial conditions in order to simplify the initialization of a model. The calling convention of a preset could not be simpler:
Because of MATLAB’s default support of tab-completion, it is trivial to cycle through all possible problems and all possible presets. All presets are callable without any arguments, with sensible defaults being applied, and all problems have a default preset named Canonical, which represents the problem specification that most commonly appears in literature, or a default designed by us. Arguments to all presets are passed with Name-Value pairs, the documentation for which appears whenever the help command is issued on the preset:
As before, the three things required for defining an initial value problem are the right hand side function, time span, and initial condition. These are very easy to access in OTP:
Accessing additional properties of the model (if it has them defined) is just as easy:
Finally, a practical example. In order to efficiently solve the problem with an implicit time integration method like ode15s built into MATLAB, one needs to pass the Jacobian to the time integrator:
iii.3 Parameter Validation
An important design principle of OTP is to protect the user. If the user does something outside of the design parameters of the problem, like inputting in impossible initial conditions or impossible parameter values, OTP is smart enough to error out, and give useful errors.
iii.4 Python Interface
Python has an interface with MATLAB. This interface does amount to writing matlab code in python, but it does make it possible to interface.
iii.5 Complex Example
We will compute the Lyapunov exponents of the Lorenz ’63 system. We will be using a large range of the APIs found in OTP.
This is really close to the best known empirical results viswanath1998lyapunov !
Iv List of Test problems
This section is a comprehensive list of and user guide to all the test problems that are included in the OTP test suite.
Each listing gives the path in MATLAB, the corresponding number of variables that the problem contains, and a short description of all the presets with some of their properties.
One of the simplest and most fundamental ODEs is the linear, first order ODE:
It is well known that the exact solution is . When eq. 2 is a scalar, constant coefficient ODE, the problem is often called to Dahlquist test problem and is commonly used to access the stability of ODE solvers dahlquist1963special .
|The Dahlquist scalar linear test problem.|
iv.2 Bouncing Ball
The bouncing ball test problem tufillaro1992experimental is meant to test event function support in time integration methods. It describes the mechanics of an inelastic collision in a perfect conservative system. The ball is modeled as a point-mass in a friction-less 2D environment.
The effect of a point mass falling to the Earth can be described by the second order ODE
where is the horizontal spacial dimension, is the vertical spacial dimension, and is the acceleration due to gravity. We can rewrite this into a 4 variable first-order ODE as follows:
It is immediately evident that this is a simple linear problem. The uniqueness of this problem comes from the event function
where is the function defining the height of the ground at the current horizontal spatial location. When the event function is triggered, the time integration method terminates, and the velocities are transformed by
where represents the angle of the ground at the current horizontal position. This boils down to rotating the velocities into a reference frame where the velocity is purely vertical, inverting it, and rotating back into the original reference frame.
This, of course, can be rewritten in the efficient form,
An immediate condition as a result of this is that the ground function is differentiable everywhere, and that there are no vertical walls anywhere in the space.
An example of the problem in action is
|This simply calls the Flat Terrain example.|
|The ball bounces on a parabola, starting slightly off-center as to create an interesting trajectory.|
|The ball has no horizontal velocity, and bounces on a perfectly flat surface.|
The Brusselator is a model for an autocatalytic reaction characterized by the reactions
This chemical reaction corresponds to the dimensionless system of first order differential equations
where are positive, real constants and are the concentrations of the two reactants ault2003dynamics .
An interesting splitting of this problem is
This is an additive splitting of the right hand side into linear and nonlinear parts, which can be used in partitioned or exponential time integrators.
|Rapid descent to a fixed point|
|A periodic cycle.|
|Rapid decay into a fixed orbit|
iv.4 Double Pendulum
The double pendulum problem models the trajectory of two pendulums connected end to end by massless rods. The first pendulum has a rod of length and a bob of mass . The second pendulum is connected to the first pendulum’s bob and has rod length and bob mass . This dynamical system is very sensitive to initial conditions and is governed by
where is the acceleration due to gravity, , and . The angles are measured counterclockwise from the negative y-axis.
This HIRES (High Irradiance Responses) problem is a mildly stiff system of eight first order differential equations. It corresponds to a chemical reaction that models how light affects morphogenesis in a plant. The equations are given by
with as the initial condition schafer1975new .
|The canonical initial conditions as used in the literature|
iv.6 Lorenz ’63 (3-variable)
The three variable Lorenz problem lorenz1963deterministic ; strogatz2018nonlinear (named here as Lorenz ’63 to contrast it with the -variable Lorenz ’96) is the first concrete ODE described to encompass chaotic dynamics.
with the canonical initial conditions outside of the trapping region:
Lorenz ’63 is a chaotic problem and, as such, is sensitive to initial conditions and parameters. In this example, we will propagate forward the same initial conditions through a problem with one of the parameters being slightly perturbed.
|This is the original problem presented in the literature. The initial condition is purposefully outside of the trapping region, but converges to it quite quickly.|
|A non-chaotic set of parameters for the problem to showcase potential periodic behavior.|
iv.7 Lorenz ’96
The Lorenz ’96 system lorenz1996predictability can be defined as follows
where with and where is a forcing function usually defined as . The canonical amount of variables of the system is .
In this discretization, one time unit corresponds to 5 days. The canonical time interval for solution is 0.05 units which corresponds to 6 hours.
Lorenz ’96 is one of the simplest -dimensional chaotic test problems. In the canonical (, ) case, the system has 13 positive Lyapunov exponents with a fractal dimension of about (popov2018bayesian, ).
Note that is a critical point, and is thus likely to be inside a trapping region. A slight perturbation from that is therefore a good choice of initial condition. Canonically, the 40 variable case,
|This is the original problem presented in the literature. The initial condition is a slight perturbation of a critical point.|
|Used in (Popov and Sandu, 2019)|
iv.8 Quasi-Geostrophic 1.5 Layer Model (QGSO)
The quasi-geostrophic model of Sakov and Oke (QGSO) sakov2008deterministic is defined by the PDE
by non-dimensionalizing potential vorticity. Equation 8 is discretized in terms of the stream function variable on the spatial domain using a -order central finite difference discretization, which by default is . The constants are taken to be , , and , with being the most variable, and potentially an order of magnitude smaller in some instances. Homogeneous Dirichlet boundary conditions are applied, and therefore the zero-valued boundary points are not explicitly stored in the computational state-space variable. The dimension is taken to be by default. The value of is calculated by the Arakawa approximation (arakawa1966computational, ; ferguson2008numerical, ). The operator is computed offline as the cube of the sparse Laplacian.
We provide several different ways of solving the Helmhotlz equation in eq. 8, which by default is solved through an offline pivoted sparse Cholesky decomposition. A multigrid solver and the ability to solve via GMRES are also provided.
Additionally, as several practical Data Assimilation applications require the computational of physical distance between state-space variables. We define the distance function as
where and are the and components of two distinct state space variables when the state space is transformed into a two dimensional grid.
This particular implementation of QGSO was written specifically for (popov2018bayesian, ).
The problem also has both Jacobian vector products and Adjoint vector products configured to work with the same linear system solvers as that of the right hand side. A Jacobian approximation is also provided.
|Uses the Gaspari-Cohn compactly supported function to build an interesting initial condition|
The Gray–Scott model gray1983autocatalytic is a reaction–diffusion problem that simulates the chemical reaction
The corresponding PDEs are given by
with periodic boundary conditions. and are diffusion rates and and are reaction rates. This problem is discretized on a uniform 2D grid with second order finite differences.
|Random initial conditions, useful for stiff integrator testing|
We provide a new MATLAB
-based framework of the creation, manipulation, and usage of problems posed as first order ordinary differential equations. OTP has an extensible, easy-to-use interface. In the future, additional problems and functionality will be added. Of special interest are more PDEs and large-scale problems.
Acknowledgements.Special thanks to Mahesh Narayanamurthi, S. Ross Glandon, Arash Sarshar, Bibek Regmi, and the rest of the Computational Science Lab at Virginia Tech for their patience and support of this project.
-  Akio Arakawa. Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I. Journal of Computational Physics, 1(1):119–143, 1966.
-  Shaun Ault and Erik Holmgreen. Dynamics of the brusselator. Math 715 Projects (Autumn 2002), page 2, 2003.
-  Germund G Dahlquist. A special stability problem for linear multistep methods. BIT Numerical Mathematics, 3(1):27–43, 1963.
-  James Ferguson. A numerical solution for the barotropic vorticity equation forced by an equatorially trapped wave. Master’s thesis, University of Victoria, 2008.
-  P Gray and SK Scott. Autocatalytic reactions in the isothermal, continuous stirred tank reactor: isolas and other forms of multistability. Chemical Engineering Science, 38(1):29–43, 1983.
-  Ernst Hairer, Syvert P Nørsett, and Gerhard Wanner. Solving ordinary differential equations. I, volume 8 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin,, 1993.
-  Ernst Hairer and Gerhard Wanner. Solving ordinary differential equations ii: Stiff and differential-algebraic problems second revised edition with 137 figures, volume 14. Springer-Verlag, 1996.
-  Edward N Lorenz. Deterministic nonperiodic flow. Journal of the atmospheric sciences, 20(2):130–141, 1963.
-  Edward N Lorenz. Predictability: A problem partly solved. In Proc. Seminar on predictability, volume 1, 1996.
-  Andrey A Popov and Adrian Sandu. A bayesian approach to multivariate adaptive localization in ensemble-based data assimilation with time-dependent extensions. arXiv preprint arXiv:1809.08984, 2018.
Pavel Sakov and Peter R Oke.
A deterministic formulation of the ensemble kalman filter: an alternative to ensemble square root filters.Tellus A: Dynamic Meteorology and Oceanography, 60(2):361–371, 2008.
-  E Schäfer. A new approach to explain the “high irradiance responses” of photomorphogenesis on the basis of phytochrome. Journal of Mathematical Biology, 2(1):41–56, 1975.
-  Steven H Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. CRC Press, 2018.
-  Nicholas B Tufillaro, Tyler Abbott, and Jeremiah Reilly. An experimental approach to nonlinear dynamics and chaos. Addison-Wesley Redwood City, CA, 1992.
-  Divakar Viswanath. Lyapunov exponents from random fibonacci sequences to the lorenzequations. Technical report, Cornell University, 1998.