 # ODE Test Problems: a MATLAB suite of initial value problems

ODE Test Problems (OTP) is an object-oriented MATLAB package offering a broad range of initial value problems which can be used to test numerical methods such as time integration methods and data assimilation (DA) methods. It includes problems that are linear and nonlinear, homogeneous and nonhomogeneous, autonomous and nonautonomous, scalar and high-dimensional, stiff and nonstiff, and chaotic and nonchaotic. Many are real-world problems from fields such as chemistry, astrophysics, meteorology, and electrical engineering. OTP also supports partitioned ODEs for testing IMEX methods, multirate methods, and other multimethods. Functions for plotting solutions and creating movies are available for all problems, and exact solutions are provided when available. OTP is desgined for ease of use-meaning that working with and modifying problems is simple and intuitive.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

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

1model = csl.otp.lorenz96.presets.Canonical;
2[t, y] = ode45(model.RHS.F, model.TimeSpan, model.Y0);
3model.movie(t, y);

A non-trivial example showcasing the versatility of the framework is

1model = csl.otp.qgso.presets.GC( ’small’);
2model.TimeSpan=[0,5000];
3[~,y1]=ode45(model.RHS.F,model.TimeSpan,model.Y0);
4model.Parameters.linearsolver=’multigrid’;
5model.Parameters.linearsolvertol=1e-8;
6[~,y2]=ode45(model.RHS.F,model.TimeSpan,model.Y0);
7disp(norm(y1(end,:)-y2(end,:))/norm(y1(end,:)));
8␣␣␣0.821952174839369

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.

## Ii Background

We consider the classic initial value problem

 d\*y(t)dt=\*f(t,\*y(t)),t∈[t0,tf]\*y(t0)=\*y0, (1)

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

### iii.1 Presets

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:

1model = csl.otp.{problemname}.presets.{Presetname};

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:

1help csl.otp.{problemname}.presets.{Presetname}

### iii.2 Usage

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:

1% define the model
2model = csl.otp.lorenz96.presets.Canonical;
3% accessing the right hand side---a function of time and state
4model.RHS.F;
5% accessing the timespan---a column vector of initial time and final time
6model.TimeSpan;
7% accessing the initial condition---a column vector
8model.Y0;
9% getting the value of the right hand side function at the initial condition
10model.RHS.F(model.TimeSpan(1), model.Y0);

Accessing additional properties of the model (if it has them defined) is just as easy:

1% defining a random column vector of the size of the model state
2rv = randn(model.NumVars, 1);
3% evaluating the adjoint model at the initial condition on the random vector

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:

1options = odeset( ’Jacobian’,model.RHS.Jacobian);
2[t,y]=ode15s(model.RHS.F,model.TimeSpan,model.Y0,options);
3model.movie(t,y);

### 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.

1model = csl.otp.lorenz63.presets.Canonical;
2model.Parameters.rho = -1;
3    Error using csl.odeutils.StructParser/checkField (line 26)
4    The field rho does not satisfy nonnegative
5
6model.Parameters.rho = [1, 1];
7    Error using csl.odeutils.StructParser/checkField (line 26)
8    The field rho does not satisfy scalar

### 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.

1>>> import matlab.engine
2>>> me = matlab.engine.start_matlab()
3>>> model = me.csl.otp.lorenz63.presets.Canonical()
4>>> me.workspace[ "model"] = model
5>>> rval = me.eval( "model.RHS.F(model.TimeSpan(1),model.Y0)")
6>>> print(rval)
710.

### 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.

1model = csl.otp.lorenz63.presets.Canonical;
2% evolve the system initially
3[~, y] = ode45(model.RHS.F, model.TimeSpan, model.Y0);
4model.Y0 = y(end, :). ’;
5%makethenexttimespanstartattheendofthepreviousandthe
6%spanitselfbe500timeunits
7model.TimeSpan=model.TimeSpan(end)+[0;500];
8%setthetolerancestothelowestallowablebyMATLABandsolve
9options=odeset(’AbsTol’,100*eps,’RelTol’,100*eps);
10[t,y]=ode45(model.RHS.F,model.TimeSpan,model.Y0,options);
11dts=diff(t);
12Q=eye(model.NumVars);
13les=zeros(model.NumVars,1);
14J=model.RHS.Jacobian(t(1),y(1,:).’);
15fori=1:(numel(dts)-1)
16␣␣␣␣%computethefulllinearlizationof
17␣␣␣␣dt=dts(i);
18␣␣␣␣%Doatrapezpidalrulestep
19␣␣␣␣W1=J*Q;
20␣␣␣␣J=model.RHS.Jacobian(t(i+1),y(i+1,:).’);
21␣␣␣␣W2=J*Q;
22␣␣␣␣%wewantourQRtobeincanonicalform.
23␣␣␣␣[Q,R]=qr(Q+(dt/2)*(W1+W2));D=diag(sign(diag(R)));Q=Q*D;R=D*R;
24␣␣␣␣rd=diag(R);les=les+log(rd);
25end
26les=les/(t(end)-t(1));
27k=sum(arrayfun(@(ki)sum(les(1:ki)),1:numel(les))>0);
28fracd=k+sum(les(1:k))/abs(les(k+1));
29disp(les);
30␣␣␣␣0.9102
31␣␣␣␣0.0027
32␣␣-14.5891
33
34disp(fracd);
35␣␣␣␣2.0626

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.

### iv.1 Linear

Path: csl.otp.linear

One of the simplest and most fundamental ODEs is the linear, first order ODE:

 dydt=A(t)\*y. (2)

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 .

Name Path Dahlquist №Vars Stiff Chaotic 1 X X

### iv.2 Bouncing Ball

Path: csl.otp.bouncingball

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

 d2xdt2=0,d2ydt2=−g, (3)

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:

 dx1dt =x2, dy1dt =y2, dx2dt =0, dy2dt =−g.

It is immediately evident that this is a simple linear problem. The uniqueness of this problem comes from the event function

 e(x1,y1,x2,y2)=y1−h(x1),

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

 [x1y1]new=[cos(θ)sin(θ)−sin(θ)cos(θ)]−1[100−1][cos(θ)sin(θ)−sin(θ)cos(θ)][x1y1],

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,

 [x1y1]new=⎡⎢ ⎢ ⎢ ⎢⎣(1−h′(x2)2)x1+2h′(x2)y11+h′(x2)22h′(x2)x1−(1−h′(x2)2)y1+1+h′(x2)2⎤⎥ ⎥ ⎥ ⎥⎦.

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

1model = csl.odetestproblems.bouncingball.presets.RandomTerrain;
2options = odeset( ’AbsTol’,1e-6,’RelTol’,1e-6,’MaxStep’,2^-4);
3[t,y]=csl.odeutils.solveeventproblem(@ode45,model,options);
4model.phaseSpacePlot(t,y); Figure 1: Modeling a perfect bouncing ball in a vacuum, on randomly generated sinusoidal terrain.
Name Path Canonical №Vars Stiff Chaotic 4 X X 4 X X X X 4 X X

### iv.3 Brusselator

Path: csl.otp.brusselator

The Brusselator is a model for an autocatalytic reaction characterized by the reactions

 A →X, 2X+Y →3X, B+X →Y+C, X →D.

This chemical reaction corresponds to the dimensionless system of first order differential equations

 dxdt=1−(b+1)x+ax2y,dydt=bx−ax2y, (4)

where are positive, real constants and are the concentrations of the two reactants ault2003dynamics .

An interesting splitting of this problem is

 f=[1−(b+1)xbx]linear+[ax2y−ax2y]nonlinear.

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.

Name Path Decay №Vars Stiff Chaotic 2 X X 2 X X 2 X X

### iv.4 Double Pendulum

Path: csl.otp.doublependulum

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

 (5)

where is the acceleration due to gravity, , and . The angles are measured counterclockwise from the negative y-axis.

### iv.5 Hires

Path: csl.otp.hires

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

 dy1dt=−1.71y1+0.43y2+8.32y3+0.0007,dy2dt=1.71y1−8.75y2,dy3dt=−10.03y3+0.43y4+0.035y5,dy4dt=8.32y2+1.71y3−1.12y4,dy5dt=−1.745y5+0.43y6+0.43y7,dy6dt=−280y6y8+0.69y4+1.71y5−0.43y6+0.69y7,dy7dt=280y6y8−1.81y7,dy8dt=−280y6y8+1.81y7, (6)

with as the initial condition schafer1975new .

Name Path Canonical №Vars Stiff Chaotic 8 ✓ X

### iv.6 Lorenz ’63 (3-variable)

Path: csl.otp.lorenz63

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.

 dxdt=σ(y−x),dydt=x(ρ−z)−y,dzdt=xy−βz, (7)

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.

1model = csl.otp.lorenz63.presets.Canonical;
2[t1, y1] = ode45(model.RHS.F, model.TimeSpan, model.Y0);
3pert = sqrt(eps(model.Parameters.sigma));
4model.Parameters.sigma = model.Parameters.sigma + pert;
5[~, y2] = ode45(model.F, t1, model.Y0);
6model.plot(t1, y1 - y2); Figure 2: Perturbations in model runs of the canonical Lorenz ’63 test problem resulting from a tiny perturbation in the parameter σ.
Name Path Canonical №Vars Stiff Chaotic 3 X ✓ 3 X X 3 X X

### iv.7 Lorenz ’96

Path: csl.otp.lorenz96

The Lorenz ’96 system lorenz1996predictability can be defined as follows

 fi=dxidt=(xi+1−xi−2)xi−1−xi+F(t),

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,

 xi={8.008i=208sonst.
Name Path Canonical №Vars Stiff Chaotic 40 X ✓ 40 X ✓

### iv.8 Quasi-Geostrophic 1.5 Layer Model (QGSO)

Path: csl.otp.qgso

The quasi-geostrophic model of Sakov and Oke (QGSO) sakov2008deterministic is defined by the PDE

 ∂q∂t=−ψx−εJ(ψ,q)−AΔ3ψ+2πsin(2πy),Δψ−Fψ=q,J(ψ,q)≡ψxqy−ψyqx, (8)

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

 d(i,j)=√(ix−jx)2+(iy−jy)2,

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. Figure 3: A typical model state of the 1.5 layer quasi-geostrophic model (QGSO). The left graph represents the stream function while the right represents vorticity.
Name Path Gaspari-Cohn №Vars Stiff Chaotic 16,129 X ✓

### iv.9 Gray–Scott

Path: csl.otp.grayscott

The Gray–Scott model gray1983autocatalytic is a reaction–diffusion problem that simulates the chemical reaction

 U+2V →3V, V →P.

The corresponding PDEs are given by

 ∂u∂t=ε1Δu−uv2+f(1−u),∂v∂t=ε2Δv+uv2−(f+k)v, (9)

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.

Name Path Random №Vars Stiff Chaotic 5,000 ✓ X

## V Conclusions

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.

## References

•  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.