Unexpected results of MATLAB's ode45

Question

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))

dadt1 = a(2);
dadt2 = -2*a(1);  
dadt3 = a(4);
dadt4 = -2*a(3) - **Fy**;

dadt = [dadt1;dadt2;dadt3;dadt4];

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)}}$

Answer 1

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.

Answer 2

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.

Answer 3

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.