MATLAB Answers

Extra independent component in ode integration affects other components

조회 수: 1(최근 30일)
Ferran
Ferran 6 Apr 2020
댓글: Ferran 7 Apr 2020
I am trying to understand how the ode functions in MATLAB work.
In this case, I am running an ode45 or ode113 with a state that contains 6 components. Secondly, I run the exact same problem, but in the function that contains the linear integration, I am adding a 7th component. This component ("extra_parameter" in the code below) is not included in any equation of the other 6 components.
Nonetheless, the resulting trajectories end up diverging from each other given enough time. This can be seen in the plot below the code. In this plot, the 7th component is also not represented in anyway.
Since this extra component is not included in the computation of the rest of the components, why does it affect their values throughout the integration?
Here's the minimum working example:
x0=[0.8 0.3 0.1 0 0.2 0]; % initial state
mu = 1.21506683e-2; % some constant I pass to the integration functions
tspan = linspace(0,200,10000); % integration time
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
% Integration
[~,y2] = ode45(@(t,y) solver1(t,y,mu),tspan,x0,options);
[~,y1] = ode45(@(t,y) solver2(t,y,mu,1),tspan,[x0 0],options); % initial value of the 7th component is 0 in this example
% Plot
figure; hold on;
plot(y2(:,1),y2(:,2),'b','LineWidth',1);
plot(y1(:,1),y1(:,2),'r--','LineWidth',1);
hold off;
% Integration functions
function dx = solver1(~,x,mu)
% Distance
r1 = sqrt((x(1)+mu)^2+x(2)^2+x(3)^2);
r2 = sqrt((x(1)+mu-1)^2+x(2)^2+x(3)^2);
% State
dx = [ x(4)
x(5)
x(6)
x(1) + 2*x(5) - (1-mu)*(x(1)+mu)/r1^3 - mu*(x(1)-1+mu)/r2^3
x(2) - 2*x(4) - (1-mu)*x(2)/r1^3 - mu*x(2)/r2^3
0 - (1-mu)*x(3)/r1^3 - mu*x(3)/r2^3];
end
function dx = solver2(~,x,mu,a)
% Distance (same as in solver 1)
r1 = sqrt((x(1)+mu)^2+x(2)^2+x(3)^2);
r2 = sqrt((x(1)+mu-1)^2+x(2)^2+x(3)^2);
% State (same as in solver 1, with an additional 7th component)
dx = [ x(4)
x(5)
x(6)
x(1) + 2*x(5) - (1-mu)*(x(1)+mu)/(r1^3) - mu*(x(1)-1+mu)/(r2^3)
x(2) - 2*x(4) - (1-mu)*(x(2))/(r1^3) - mu*(x(2))/(r2^3)
0 - (1-mu)*(x(3))/(r1^3) - mu*(x(3))/(r2^3);
a]; % extra parameter; so that a_new = a_previous + a * time throughout the integration
end
Plot:

  댓글 수: 5

표시 이전 댓글 수: 2
Ferran
Ferran 6 Apr 2020
The difference starts to appear when I integrate for longer times. In the figure I posted, both plots look the same until at some point they are start diverging and no longer overlap.
darova
darova 6 Apr 2020
Indeed, there is a difference. ode45 has some adaptive algorithms for choosing step of integration (dt). Maybe it chooses slightly different steps for these functions
I can't explain it honestly. It's weird. Maybe this is some of those cases when people say "tolerance" or "machine precision"
Ferran
Ferran 7 Apr 2020
Seeing James' answer below, I understand that it is in the nature of adaptive step ode solvers. I tried other odes and obtained similar results. At least it's good to now know this and be aware of that when using the output of any integration. Thank you, darova!

로그인 to comment.

채택된 답변

James Tursa
James Tursa 6 Apr 2020
The integrators do not internally step each element separately. They still step the entire state vector as a whole. So even though the derivative equations are independent, they will all affect the overall internal stepping that the integrator does to arrive at a solution. And if you change the stepping, you will change the result for all elements, even if only slightly.

  댓글 수: 3

Ferran
Ferran 7 Apr 2020
Thank you, James. Do you have any suggestion as to manually fixing the timestep or some other solution in that direction, or would a non-variable timestep produce even worse results? Not sure if this makes sense, I am new to this.
James Tursa
James Tursa 7 Apr 2020
The ode integrators do not give the user the option of setting the step sizes directly. All of the step sizing is adaptive and done internally based on the tolerances in effect ('RelTol' and 'AbsTol'). Slight changes in the conditions of the integration should be expected to produce slight changes in the results, but within the tolerances. For long duration integrations, all methods may eventually diverge from each other. This is particularly true of "orbit" type of ode's, where the rectangular coordinate systems used for the "circular'ish" trajectories tend to produce systematic integration errors that build up over time (often downtrack errors). To get consistency for your ode's being solved, tighter tolerances or different / higher order methods may be appropriate, but you will still get some divergence.

로그인 to comment.

추가 답변(0개)

이 질문에 답변하려면 로그인을(를) 수행하십시오.

태그

제품


릴리스

R2018b

Translated by