# Unexpected results of MATLAB's ode45

Whilst working with MATLAB recently I encountered something odd that I cannot explain. I was using the ode45 solver to solve a system of two coupled second order ODEs. I wasn't convinced about the results, so I tried something easier, just to see where it goes wrong. I uncoupled the equations and I swapped one of them to be a simple harmonic oscillator.

Quite surprisingly, I did not get what I expected. However when I removed one of the terms from the second equation (which is now uncoupled!) and ran it again, it produced a correct answer. This blows my mind, since there is no link at all between the two equations now. Here is my code after alteration:

Main function:

function Ord2ODE

t=0:0.0001:20;
%x,xdot,y,ydot
ainit = [0; 1; 0; 0];

[t,a] = ode45(@rhs,t,ainit);

figure;
plot(t, a(:,1));

end


and the function rhs:

function dadt = rhs(t,a)

mu = 0;

x = a(1);
xdot = a(2);
y = a(3);
ydot = a(4);

**Fy = stuff in terms of x, y (a(1), a(2))**
Fx = stuff in terms of x, y (a(1), a(2))

end


with the prblematic bit marked with ** ** (definition and occurance in computation of dadt4). You can see that neither a(3) nor a(4) are present in the calculation of x. The results with and without that term are posted below. Does anyone have any idea why would a term in an uncoupled equation cause this kind of divergence in a solution of the other equation?

Fy = -mu*y - (1-mu)*y + mu*y/(sqrt(((x+1-mu)^2 + y^2)^3)) + (1-mu)*y/(sqrt(((x-mu)^2 + y^2)^3));

Fx = -mu*(x+1-mu) - (1-mu)*(x-mu) + mu*(x+1-mu)/(sqrt(((x+1-mu)^2 + y^2)^3)) + (1-mu)*(x-mu)/(sqrt(((x-mu)^2 + y^2)^3));


$Fy = -\mu*y - (1-\mu)*y + \frac{\mu*y}{(\sqrt{((x+1-\mu)^2 + y^2)^3)}} + \frac{(1-\mu)*y}{(\sqrt{((x-\mu)^2 + y^2)^3)}}$

$Fx = -\mu*(x+1-\mu) - (1-\mu)*(x-\mu) + \frac{\mu*(x+1-\mu)}{(\sqrt{((x+1-\mu)^2 + y^2)^3)}} + \frac{(1-\mu)*(x-\mu)}{(\sqrt{((x-\mu)^2 + y^2)^3)}}$

• Unless there is a good reason for excluding the information, I recommend you include the full expressions for Fx and Fy in your code. Otherwise, there's really no way to know what is going wrong. – Geoff Oxberry Sep 16 '15 at 19:25
• @GeoffOxberry done now, they were just a bit longish. Does that help? – Glo Sep 16 '15 at 19:34
• I don't see Fx used at all in your rhs ODE function. – horchler Sep 16 '15 at 23:02
• @horchler that is because after I obtained something dubious whilst working, I simplified down my equations, as mentioned in the original post. The point is that the first equation is just a SHM, and the second one is whatever but the important thing is it's uncoupled from the first one. So the set of equations: ($\frac{d^2x}{dt^2} = -2x$ ; $\frac{d^2y}{dt^2} = -2\frac{dy}{dt} - Fy$) gives different result for x from what ($\frac{d^2x}{dt^2} = -2x$ ; $\frac{d^2y}{dt^2} = -2\frac{dy}{dt}$) gives. See my point? The change in 2nd equation should not influence the solution for x. – Glo Sep 17 '15 at 7:40

From what I understand from your question, the first plot has become unstable, so perhaps if you did a plot with the two methods for the first time period (so up to around 4 seconds) your two solutions will match (looking at your plots they seem to match). If that was the case, then try solving the problem with a different solver, or making your time steps smaller to ensure that the problem is not related to the time stepping scheme.

• it would not be a problem if the equation simply become unstable. The thing that worries me is that there is no coupling between equation one and two (or, after I transform into first order equations - between equations 1, 2 and equations 3, 4). How can an equation that is not coupled to the other equation affect the result? It's as if I ran my program on the first equation, then ran it again on another one and two independent runs affected each other. Note that dadt4 has no link to dadt1 or dadt2! – Glo Sep 16 '15 at 14:35
• the structure of the bit of the code u've put here seems fine. Another thought, this may not be a problem in you case, but when using ode solvers you are putting all the equations in the form dy/dt = f(y,t), so essentially they are being put into one system of equations, even if the actual equations are not coupled, they form a single system, that may be unstable. I can't think of anything else other than try other solvers and rule out different possibilities, or go back to the code and look for mistakes. good luck. – Hooman Sep 16 '15 at 16:03
• @Glo Continuing on Hooman's comment, the ode45 function in MATLAB uses dynamic time steps, so if you solve two uncoupled ODE at the same time, then both will affect the step size of the time steps. – fibonatic Sep 18 '15 at 2:22
• @fibonatic but I specified t = 0:0.001:10. Doesn't that fix the timestep? I tried much smaller timesteps, too. – Glo Sep 18 '15 at 10:26
• @Glo That only fixes some points of the time vector, depending on the problem there might be additional time steps in between. You can test this by looking at the length of the vectors that get returned by MATLAB. – fibonatic Sep 18 '15 at 10:30

Here is a guess as to what might be going on. ode45 (like all the MATLAB ode solvers) adjusts the step size based on it's estimate of the solution error. All equations are considered in this calculation whether they are coupled or not.

The solution errors are compared to two parameters, AbsTol and RelTol. Both of these have default values but sometimes these defaults don't produce an accurate solution. I suggest experimenting with these two parameters (e.g. make them an order of magnitude smaller) and see how it affects your solution.

• I just tried that - even very small values (1e-10) give the same result. As far as I can tell changing the stepsize helps a little, but it's not worth it - execution takes 10^4 times longer and the result barely changes. Does my code seem correct? Perhaps there's something obvious I'm missing and the equations somehow mix? – Glo Sep 16 '15 at 16:28

Although both equations are uncoupled in the absence of Fy, they are both being solved together and if there is a problem solving one equation this could be a problem solving the other.

I copied your code and ran it. The equation for y is unstable and quickly returns NaN. However, for me, the result for x is the same whether I include Fy or not.

• Really? Could it be that it's MATLAB's fault, i.e. the ODE solver got updated at some point? I'm using R2013a version of MATLAB. Did you make any alterations to the code (error tolerance etc.) and are you sure you plotted the same thing each time (a(1) against t) or that you cleared the variables and not plotted the old variable? Sorry, I'm just trying to think of all the possibilities! – Glo Sep 17 '15 at 12:48
• I'm not sure what the problem is but it cannot be the MALTAB version. – Hooman Sep 17 '15 at 14:05
• I've copied it exactly as you have it above (substituting in your expression for Fy) and I definitely have not made any mistake like plotting previous variables. I've run it in R2008b and R2007b and it plots exactly what you have shown in your second plot above. – jamble Sep 17 '15 at 15:59
• I have now tried it again in R2012a and it plots what you have shown in your first plot. So it must be something to do with the version of Matlab. – jamble Sep 17 '15 at 19:56
• It's strange that an older version is producing better results! Out of interest, have u guys tried using other solvers than ode45, for instance ode15s, would that change the results? – Hooman Sep 17 '15 at 20:11