Practical Advice for People on the Leading Edge

The new solution framework for Ordinary Differential Equations (ODEs) in MATLAB R2023b

Along with linear algebra, one of the iconic features of MATLAB in my mind is how it handles ordinary differential equations (ODEs). ODEs have been part of MATLAB almost since the very beginning.
One of the features of how MATLAB traditionally allows users to solve ODEs is that it provides a suite of functions. For many years, there were 7 routines in the MATLAB suite, described in this 2014 post by Cleve Moler and in more detail in the 1997 paper, The MATLAB ODE Suite. In 2021b, 2 new high-order methods were added to the suite, ode78 and ode89 based on the algorithms described in the paper Numerically optimal Runge–Kutta pairs with interpolants.
The design of the suite was elegant and has served the community well for over 25 years! However, MATLAB has evolved a lot since it was designed and user-expectations have evolved with it. It was felt that it was time for a fresh look at how to solve ODEs in MATLAB, one that would additionally support our future plans in a modern and elegant manner.
Before I dive into the new interface, I'd like to make it clear that the existing suite of ODE functions are not going anywhere! Millions of lines of code make use of ode45 and friends and we have no plans on doing anything that would break that code. This is about adding functionality, not taking anything away.
Furthermore, the focus of this release is the new interface. There are no new solvers or any fundamentally new functionality....yet! However, some things will be significantly easier to do using the new interface.
With that said, let's take a look at how the new design looks.

Solving ODEs in MATLAB the OOP way

Say we want to solve the ODE
with the initial condition y0=0 over the interval t=[0,10]
With the new interface, we can set up and solve this problem as follows
F = ode(ODEFcn=@(t,y) 2*t,InitialTime=0,InitialValue=0); % Set up the problem by creating an ode object
sol = solve(F,0,10); % Solve it over the interval [0,10]
plot(sol.Time,sol.Solution,"-o") % Plot the solution
An alternative way to proceed would have been to start with an empty ode object and add one property at a time:
F = ode; % Empty ode object called F
F.ODEFcn = @(t,y) 2*t; % Add the function we want to solve to F
F.InitialTime = 0; % Add the initial time to F
F.InitialValue = 0; % Add the initial value to F
Once our problem is set up, we pass it to the solve function.
sol = solve(F,0,10);
The solution is an ODEResults object with two properties: Time and Solution
sol =
ODEResults with properties: Time: [0 0.2500 0.5000 0.7500 1 1.2500 1.5000 1.7500 2 2.2500 2.5000 2.7500 3 3.2500 3.5000 3.7500 4 4.2500 4.5000 4.7500 5 5.2500 5.5000 5.7500 6 6.2500 6.5000 6.7500 7 7.2500 7.5000 7.7500 8 8.2500 8.5000 8.7500 9 9.2500 9.5000 9.7500 10] Solution: [0 0.0625 0.2500 0.5625 1.0000 1.5625 2.2500 3.0625 4 5.0625 6.2500 7.5625 9 10.5625 12.2500 14.0625 16 18.0625 20.2500 22.5625 25 27.5625 30.2500 33.0625 36 39.0625 42.2500 45.5625 49 52.5625 56.2500 60.0625 64 68.0625 … ] (1×41 double)
You access either property with . notation
ans = 1×41
0 0.0625 0.2500 0.5625 1.0000 1.5625 2.2500 3.0625 4.0000 5.0625 6.2500 7.5625 9.0000 10.5625 12.2500 14.0625 16.0000 18.0625 20.2500 22.5625 25.0000 27.5625 30.2500 33.0625 36.0000 39.0625 42.2500 45.5625 49.0000 52.5625 56.2500 60.0625 64.0000 68.0625 72.2500 76.5625 81.0000 85.5625 90.2500 95.0625 100.0000
ans = 1×41
0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500 3.0000 3.2500 3.5000 3.7500 4.0000 4.2500 4.5000 4.7500 5.0000 5.2500 5.5000 5.7500 6.0000 6.2500 6.5000 6.7500 7.0000 7.2500 7.5000 7.7500 8.0000 8.2500 8.5000 8.7500 9.0000 9.2500 9.5000 9.7500 10.0000

Automatic solver selection

One of the first things to note in the above workflow is that we didn't choose a solver. That is, I didn't have to think about if I should use ode45, ode78 etc. All I did was state the problem mathematically and ask MATLAB to solve it. Let's take a closer look at the details.
I set up my problem like this:
F = ode(ODEFcn=@(t,y) 2*t,InitialTime=0,InitialValue=0);
F is an ode object:
ans = 'ode'
I can see the details by evaluating F
F =
ode with properties: Problem definition ODEFcn: @(t,y)2*t InitialTime: 0 InitialValue: 0 Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: auto SelectedSolver: ode45 Show all properties
The object display separates the mathematical problem definition from how we are going to solve it. By default, the Solver property of the class is set to "auto" which means that MATLAB attempts to choose a suitable solver based on various properties of the problem we've asked it to solve. In this case it has chosen ode45 which is a pretty safe bet for many problems.
Ask for a tighter tolerance for this problem, however, and it switches to using ode78.
F.RelativeTolerance = 1e-7
F =
ode with properties: Problem definition ODEFcn: @(t,y)2*t InitialTime: 0 InitialValue: 0 Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-07 Solver: auto SelectedSolver: ode78 Show all properties
At this stage, MATLAB hasn't solved anything. We've just set up our problem and made decisions about how we are going to solve it. We need the solve function to complete the job.
One thing to bear in mind is that the solver that "auto" chooses may change in future releases for a number of reasons. For example, we may improve the heuristics used to select the best solver or maybe add new solvers that do a better job for a given problem.
If you want to override what MATLAB chooses and fix the solver type you can do that as follows
F.Solver = "ode45" % Force the framework to use ode45.
F =
ode with properties: Problem definition ODEFcn: @(t,y)2*t InitialTime: 0 InitialValue: 0 Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-07 Solver: ode45 Show all properties

Trying a stiff ODE problem

I was curious about how it would handle a stiff ODE and so fed it the example given by Cleve in his blog post Ordinary Differential Equations, Stiffness.
delta = 0.0001;
CleveF = ode(ODEFcn=@(t,y) y^2 - y^3);
CleveF.InitialTime = 0;
CleveF.InitialValue = delta;
sol = solve(CleveF,0,2/delta);
Looks reasonable. Following Cleve's suggestion, we zoom in on the steady state that begins at x=1 and see that the solver is working hard to do its job, just as it did when Cleve explicitly chose ode45.
xlim([10006.3 10112.1])
ylim([0.99615 1.00250])
Sure enough, the "auto" option has also chosen ode45 in this case.
ans =
SolverID enumeration ode45
Using the old solver suite, I'd now have to dig into the documentation and read about the stiff solvers available to me before trying those out. With the new interface, however, I just tell it that I think the problem is stiff and it will choose a stiff solver for me.
CleveF.Solver = "stiff"
CleveF =
ode with properties: Problem definition ODEFcn: @(t,y)y^2-y^3 InitialTime: 0 InitialValue: 1.0000e-04 Jacobian: [] Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: stiff SelectedSolver: ode15s Show all properties
This time it has selected ode15s for me. Let's see how that does
stiffSol = solve(CleveF,0,2/delta);
Zooming in at x=1:
xlim([9620.2 9748.2])
ylim([0.99464 1.00232])
Much better behaved!
Being able to tell MATLAB that the problem is "stiff" makes life a little easier than before but I was a disappointed that "auto" didn't realise that my problem was stiff and choose a relevant solver for me.
I reached out to development to ask why it failed me. They told me that the automatic solver selection is a heuristic that operates without peering inside your equations. It reacts to the data you supply to define the problem, e.g. if you've supplied a Jacobian, what the relative tolerances are, etc. There's no attempt to diagnose stiffness at all. That's why there are "auto", "stiff", and "nonstiff" automatic modes.
Of course, this may change in future releases. Now this feature exists, it will be possible to improve on it and its already pretty useful!

More flexible in time

Consider the following classic system of first-order ODEs that describes simple harmonic motion
$ \frac{dy_1}{dt} = y_2 \ \frac{dy_2}{dt} = -y_1 $
subject to the initial conditions y1(0)=1and y2(0)=0
With the classic suite of solvers, you might solve this as follows
y0 = [1; 0]; % initial conditions: y1(0)=1, y2(0)=0
tspan = [0 10]; % time interval from 0 to 10 for the solution
[t,y] = ode45(@(t,y)[y(2);-y(1)], tspan,y0);
So far so simple. Imagine now that I want to go backwards in time from my initial condition as well as forwards. That is, I want my time span to be
tspan2 = [-10 10]; % time interval from 0 to 10 for the solution
This is a little complicated using ode45 since it assumes that the initial condition you set using y0 corresponds to the beginning of the span. That is, if I do
[t,y] = ode45(@(t,y)[y(2);-y(1)], tspan2,y0);
title("Not the boundary condition I wanted!")
it applies the boundary condition y1(-10)=1and y2(-10)=0 which is not the problem I wanted it to solve! Instead I have to make two calls to ode45 -- one that goes forwards in time and the other that goes backwards.
y0 = [1; 0]; % initial conditions: y1(0)=1, y2(0)=0
[tf,yf] = ode45(@(t,y)[y(2);-y(1)], [0 10],y0);
[tb,yb] = ode45(@(t,y)[y(2);-y(1)], [0 -10],y0);
We got there but it took a couple of backflips. With the new interface, I don't need to worry about this as its much more flexible.
shm = ode(ODEFcn=@(t,y)[y(2);-y(1)],InitialTime=0,InitialValue=[1;0],Solver="ode45");
shmSol = solve(shm,-10,10);
title("The same result but we got there more easily")
Unlike the traditional suite, the solution span we ask for doesn't even need to include our initial conditions.
% Request a solution span that doesn't include the initial conditions at t=0
shmSol2 = solve(shm,2*pi,4*pi)
shmSol2 =
ODEResults with properties: Time: [6.2832 6.4381 6.6484 6.8587 7.0690 7.2792 7.4806 7.6820 7.8833 8.0847 8.2869 8.4891 8.6913 8.8935 9.0955 9.2975 9.4995 9.7015 9.8989 10.0962 10.2935 10.4908 10.6933 10.8959 11.0984 11.3009 11.4952 11.6896 11.8839 12.0783 … ] (1×34 double) Solution: [2×34 double]
The closeness of the solution points around t = 0 is the solver starting the two integrations with a small initial step size. If we would rather choose the output points ourselves, we can just specify them.
shmSol3 = solve(shm,linspace(-pi,pi));
We can get all of these results using the traditional suite, of course, it's just that its easier and more intuitive now.

Event detection

Event detection has always been part of the MATLAB ODE suite and so, of course, this is also possible with the new framework. The documentation for ode demonstrates how to solve the classic bouncing ball problem using event detection. Here, I'll demonstrate how to use it to find maxima and minima of the solution to an ODE problem. The original problem comes from a 2011 blog post by John Kitchin.
First, let's solve the ODE y=e-0.05tsin(t)
myODE = ode(ODEFcn=@(t,y) exp(-0.05*t)*sin(t),InitialTime=0,InitialValue=0);
sol = solve(myODE,0,20,Refine=10);
We want to find the maxima and minima of this solution. We know from calculus that these occur when our original equation y=e-0.05tsin(t)is equal to zero.
odeEvent objects define events that the solver detects while solving an ordinary differential equation. An event occurs when one of the event functions you specify crosses zero.
So, we create an odeEvent object with an event function equal to our original ODE.
E = odeEvent(EventFcn = myODE.ODEFcn)
E =
odeEvent with properties: EventFcn: @(t,y)exp(-0.05*t)*sin(t) Direction: both Response: proceed CallbackFcn: []
This time, let's solve the ODE with this event definition
myODE = ode(ODEFcn=@(t,y) exp(-0.05*t)*sin(t),InitialTime=0,InitialValue=0,EventDefinition=E);
sol = solve(myODE,0,20,Refine=10)
sol =
ODEResults with properties: Time: [0 0.2000 0.4000 0.6000 0.8000 1 1.2000 1.4000 1.6000 1.8000 2 2.2000 2.4000 2.6000 2.8000 3 3.2000 3.4000 3.6000 3.8000 4 4.2000 4.4000 4.6000 4.8000 5 5.2000 5.4000 5.6000 5.8000 6 6.2000 6.4000 6.6000 6.8000 7 … ] (1×101 double) Solution: [0 0.0205 0.0797 0.1733 0.2969 0.4451 0.6122 0.7917 0.9766 1.1594 1.3319 1.4872 1.6201 1.7256 1.8002 1.8417 1.8489 1.8223 1.7635 1.6754 1.5622 1.4285 1.2800 1.1238 0.9664 0.8139 0.6717 0.5448 0.4377 0.3544 0.2983 … ] (1×101 double) EventTime: [2.8422e-14 3.1416 6.2832 9.4248 12.5664 15.7080 18.8496] EventSolution: [4.2572e-28 1.8503 0.2696 1.6201 0.4659 1.4521 0.6088] EventIndex: [1 1 1 1 1 1 1]
The solution includes all of the times the event occurs. Let's plot them
hold on
hold off
Exactly what we were looking for.

Obtaining the ODE solution as a function

Until now, we have been working with discrete representations of the solutions to our ODE problems and this can present limitations on further analysis. For example, what if we wanted to integrate the solution? What we need is a function that represents the solution and we can get one using solutionFcn, passing in the ode object we are aiming to solve.
myODE = ode(ODEFcn=@(t,y) exp(-0.05*t)*sin(t),InitialTime=0,InitialValue=0); % Define the ODE problem
solFcn = solutionFcn(myODE,0,20); % Use it to create a solutionFcn
I can now evaluate the solution at any point I like
ans = 1.3404
or pass it to MATLAB's integral function
ans = 19.5314
Again, you could have done something similar with the old framework using deval and friends but I much prefer the new way!

Mass Matrices, Jacobians, Parameters and everything else I haven't covered here

This is just an introduction to the new framework and I encourage you to read the documentation to explore everything else it offers. Let me know what you think in the comments section.
  • print


To leave a comment, please click here to sign in to your MathWorks Account or create a new one.