{"id":2937,"date":"2024-12-11T09:33:44","date_gmt":"2024-12-11T14:33:44","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=2937"},"modified":"2024-12-11T09:33:44","modified_gmt":"2024-12-11T14:33:44","slug":"solving-higher-order-odes-in-matlab","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2024\/12\/11\/solving-higher-order-odes-in-matlab\/","title":{"rendered":"Solving Higher-Order ODEs in MATLAB"},"content":{"rendered":"
In a comment<\/span><\/a> to last year's introduction to the new ODE solution framework in MATLAB, Ron asked if I could provide an example of using it to solve a 2nd order ODE since most tutorials don't deal with the new syntax. Absolutely I can! <\/span><\/div>
So let's start at the beginning. An <\/span>ordinary differential equation<\/span> (ODE) contains one or more derivatives of a dependent variable, <\/span>y<\/span>, with respect to a single independent variable, <\/span>t<\/span>, usually referred to as time. The notation used here for representing derivatives of <\/span>y<\/span> with respect to <\/span>t<\/span> is <\/span>y<\/span>\u2032 for a first derivative, <\/span>y<\/span>\u2032\u2032 for a second derivative, and so on. The <\/span>order<\/span> of the ODE is equal to the highest-order derivative of <\/span>y<\/span> that appears in the equation.<\/span><\/div>
This equation then is second order<\/span><\/div>
my<\/mi>\u2032<\/mo>\u2032<\/mo>+<\/mo>ky<\/mi>=<\/mo>0<\/mn><\/mrow><\/math><\/span><\/div>
MATLAB's current crop of numerical ODE solvers only solve first order equations so we have work to do before we can proceed. It turns out that <\/span>you can write any higher-order ODE as a system of first-order ODEs<\/span>. <\/span><\/div>

Calculating the system of ODEs by hand<\/span><\/h2>
You first make the following definitions<\/span><\/div>
y<\/mi><\/mrow>1<\/mn><\/mrow><\/msub>=<\/mo>y<\/mi><\/mrow><\/math><\/span><\/div>
y<\/mi><\/mrow>2<\/mn><\/mrow><\/msub>=<\/mo>y<\/mi><\/mrow>\u2032<\/mo><\/mrow><\/msup><\/mrow><\/math><\/span><\/div>
Then we take their derivatives<\/span><\/div>
y<\/mi><\/mrow>1<\/mn><\/mrow><\/msub><\/mrow>\u2032<\/mo><\/mrow><\/msup>=<\/mo>y<\/mi><\/mrow>\u2032<\/mo><\/mrow><\/msup><\/mrow><\/math><\/span><\/div>
y<\/mi><\/mrow>2<\/mn><\/mrow><\/msub><\/mrow>\u2032<\/mo><\/mrow><\/msup>=<\/mo>y<\/mi><\/mrow>\u2032<\/mo>\u2032<\/mo><\/mrow><\/msup><\/mrow><\/math><\/span><\/div>
Do substitutions using our definitions and our original second order equation<\/span><\/div>
y<\/mi><\/mrow>1<\/mn><\/mrow><\/msub><\/mrow>\u2032<\/mo><\/mrow><\/msup>=<\/mo>y<\/mi><\/mrow>2<\/mn><\/mrow><\/msub><\/mrow><\/math><\/span><\/div>
y<\/mi><\/mrow>2<\/mn><\/mrow><\/msub><\/mrow>\u2032<\/mo><\/mrow><\/msup>=<\/mo>-<\/mo>k<\/mi><\/mrow>m<\/mi><\/mrow><\/mfrac>y<\/mi><\/mrow>1<\/mn><\/mrow><\/msub><\/mrow><\/math><\/span><\/div>
and we have our system of first order ODEs. We write this system as the following function<\/span><\/div>
% Define the system of first-order ODEs<\/span><\/span><\/div><\/div>
k = 10; <\/span>% Spring Constant<\/span><\/span><\/div><\/div>
m = 1; <\/span>% Mass<\/span><\/span><\/div><\/div>
odeSystem = @(t, y) [y(2); -k\/m * y(1)];<\/span><\/span><\/div><\/div><\/div>
So far, nothing is new. This is the exact same method you've always had to follow for doing this in MATLAB.<\/span><\/div>
We plug <\/span>odeSystem <\/span>into the new framework as follows, also stating our initial conditions which are <\/span>y=1<\/span> and <\/span>y'=0<\/span> at <\/span>t=0<\/span>. This is the only step that's different from any tutorials you might have seen before.<\/span><\/div>
shm = ode(ODEFcn=odeSystem,InitialTime=0,InitialValue=[1;0]);<\/span><\/span><\/div><\/div>
shmSol = solve(shm,-5,5);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
plot(shmSol.Time,shmSol.Solution,<\/span>'-o'<\/span>);<\/span><\/span><\/div><\/div>
legend([<\/span>\"y\"<\/span>,<\/span>\"y'\"<\/span>])<\/span><\/span><\/div><\/div>
title(<\/span>\"Solution of the 2nd order system and first derivative\"<\/span>)<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"t\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>

Making use of a solutionFcn<\/span><\/h3>
Something to point out here is that our starting point is at <\/span>t=0 <\/span>and that since we asked for a range from <\/span>t=-5<\/span> to <\/span>t=5<\/span> we are integrating in both directions automatically, allowing the solver to chose whatever values of <\/span>t<\/span> it deems are best. As you can see from the plot above, the solver used a lot of points around <\/span>t=0<\/span> in both directions. <\/span><\/div>
We can get more control by defining a <\/span>solutionFcn<\/span><\/a>. Then we can see the solution at a set of uniform points in time<\/span><\/div>
y = solutionFcn(shm,-5,5,OutputVariables=1); <\/span>% Represents the solution of our system<\/span><\/span><\/div><\/div>
t = linspace(-5,5,100);<\/span><\/span><\/div><\/div>
plot(t,y(t),<\/span>'-o'<\/span>)<\/span><\/span><\/div><\/div>
title(<\/span>\"Solution of the 2nd order system\"<\/span>)<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"t\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
To recover the derivative, use the second output of the <\/span>solutionFcn,<\/span> <\/span>y<\/span><\/div>
[yt,ypt] = y(t);<\/span><\/span><\/div><\/div>
plot(t,yt,<\/span>'-o'<\/span>,t,ypt,<\/span>'-o'<\/span>)<\/span><\/span><\/div><\/div>
legend([<\/span>\"y\"<\/span>,<\/span>\"y'\"<\/span>])<\/span><\/span><\/div><\/div>
title(<\/span>\"Solution of the 2nd order system and first derivative\"<\/span>)<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"t\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>

Calculating the system of ODEs using Symbolic Toolbox<\/span><\/h2>
I'm not going to lie to you, I don't trust it when I do mathematics by hand. Even for something as simple as this I make mistakes and for a more complex equation you can bet good money that I'll mess it up. I much prefer it when I can get the computer to do the heavy lifting for me. That's what they are for! <\/span><\/div>
I started off by asking development why we don't 'just' support higher order ODEs directly in the new framework. Once I put my dollar in the 'Just Jar', they told me that its something they've thought about off and on for a while. Could I get more user stories where people would value such functionality? So here I am asking, dear reader, would such direct support be useful for you and, if so, what's your application?<\/span><\/div>
In the meantime, I found a page in the documentation that <\/span>shows how to use symbolic toolbox to do the transformation to a system of ODEs for us<\/span><\/a> using the <\/span>odeToVectorField<\/span><\/a> function. Let's apply it to the problem here. Recall, this is the problem we want to solve<\/span><\/div>
my<\/mi>\u2032<\/mo>\u2032<\/mo>+<\/mo>ky<\/mi>=<\/mo>0<\/mn><\/math><\/span><\/div>
We set up some symbolic variables<\/span><\/div>
syms <\/span>y(t) m k <\/span><\/span><\/div><\/div><\/div>
and pass our equation to <\/span>odeToVectorField<\/span><\/a> using the Symbolic Toolbox function <\/span>diff<\/span> to handle the derivatives<\/span><\/div>
[V] = odeToVectorField(m*diff(y, 2) + k*y ==0)<\/span><\/span><\/div>
V = <\/div><\/span><\/div><\/div><\/div><\/div>
That result looks familiar! It's the same system I derived by hand earlier. We need to convert it to a MATLAB function so I can pass it to the numerical ODE framework<\/span><\/div>
M = matlabFunction(V,<\/span>'vars'<\/span>,{<\/span>'t'<\/span>,<\/span>'Y'<\/span>,<\/span>'m'<\/span>,<\/span>'k'<\/span>})<\/span><\/span><\/div>
M = function_handle with value:<\/span><\/span><\/div>
@(t,Y,m,k)[Y(2);-(k.*Y(1)).\/m]\r\n<\/div><\/div><\/div><\/div><\/div><\/div>
Now I need to sort out the parameters <\/span>m<\/span> and <\/span>k<\/span>. The function above has four input parameters but I want one that only takes two.<\/span><\/div>
k = 10; <\/span>% Spring Constant<\/span><\/span><\/div><\/div>
m = 1; <\/span>% Mass<\/span><\/span><\/div><\/div>
odeSystemFromSymbolic = @(t, y) M(t,y,m,k)<\/span><\/span><\/div>
odeSystemFromSymbolic = function_handle with value:<\/span><\/span><\/div>
@(t,y)M(t,y,m,k)\r\n<\/div><\/div><\/div><\/div><\/div><\/div>
we pass this into the ODE framework as before<\/span><\/div>
shm = ode(ODEFcn=odeSystemFromSymbolic,InitialTime=0,InitialValue=[1;0]);<\/span><\/span><\/div><\/div>
shmSol = solve(shm,-5,5);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
plot(shmSol.Time,shmSol.Solution,<\/span>'-o'<\/span>);<\/span><\/span><\/div><\/div>
legend([<\/span>\"y\"<\/span>,<\/span>\"y'\"<\/span>])<\/span><\/span><\/div><\/div>
title(<\/span>\"Solution of the 2nd order system\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
This gives the identical result but I didn't need to any algebra by hand, which is always a good thing! Of course, we are also free to use any of the other machinery in the new ODE framework including <\/span>solutionFcn<\/span> but I won't show that again for the sake of space.<\/span><\/div>\r\n<\/div>