I am running a code that runs an ode in which the equations of motion are exactly the same for both cases (except the line where I am making the little change). With the same input, I obtain different output.
Here's a simplified working version of the main code. The only difference between solver_ex1 and solver_ex2 is that in the calculation of r23, I use (x(1)+mu-1)^2 versus (x(1)-1+mu)^2. mu is a constant.
mu = 1.21506683e-2;x0 = [0.8 0.3 0.1 0 0.2 0];tspan = linspace(0,200,1000);options = odeset('RelTol',1e-12,'AbsTol',1e-12);[~,y1] = ode113(@(t,y) solver_ex1(t,y,mu),tspan,x0,options);[~,y2] = ode113(@(t,y) solver_ex2(t,y,mu),tspan,x0,options);figure; hold on;plot(y1(:,1),y1(:,2),'g-','LineWidth',0.5);plot(y2(:,1),y2(:,2),'m--','LineWidth',0.5);function xout = solver_ex1(~,x,mu)% Distances (in r^(3/2) format)
r13 = ( (x(1)+mu-0)^2 + x(2)^2 + x(3)^2 )^1.5;r23 = ( (x(1)+mu-1)^2 + x(2)^2 + x(3)^2 )^1.5;% New state (nondimensional equations of motion)
dx = [ x(4) x(5) x(6) x(1) + 2*x(5) - (1-mu)*(x(1)+mu)/r13 - mu*(x(1)-1+mu)/r23 x(2) - 2*x(4) - (1-mu)*x(2)/r13 - mu*x(2)/r23 0 - (1-mu)*x(3)/r13 - mu*x(3)/r23];xout(1:6) = dx(1:6);xout = xout';endfunction xout = solver_ex2(~,x,mu)% Distances (in r^(3/2) format)r13 = ( (x(1)+mu-0)^2 + x(2)^2 + x(3)^2 )^1.5;r23 = ( (x(1)-1+mu)^2 + x(2)^2 + x(3)^2 )^1.5;% New state
dx = [ x(4) x(5) x(6) x(1) + 2*x(5) - (1-mu)*(x(1)+mu)/r13 - mu*(x(1)-1+mu)/r23 x(2) - 2*x(4) - (1-mu)*x(2)/r13 - mu*x(2)/r23 0 - (1-mu)*x(3)/r13 - mu*x(3)/r23];xout(1:6) = dx(1:6);xout = xout';end
The output differs in each case. Plotting them as in the code also gives different visual trajectories after a while. That is, they are similar in the beginning, but start to diverge later on.
I get that it might be the tolerance, and the fact that it goes near the singularity (for r23), but shouldn't MATLAB deal with both cases the same way?
I am a bit confused at the moment. Could you give me a hint on what's going on?
Best Answer