gusucode.com > demos工具箱matlab源码程序 > demos/orbitode.m
function orbitode %ORBITODE Restricted three body problem. % This is a standard test problem for non-stiff solvers stated in Shampine % and Gordon, p. 246 ff. The first two solution components are % coordinates of the body of infinitesimal mass, so plotting one against % the other gives the orbit of the body around the other two bodies. The % initial conditions have been chosen so as to make the orbit periodic. % Moderately stringent tolerances are necessary to reproduce the % qualitative behavior of the orbit. Suitable values are 1e-5 for RelTol % and 1e-4 for AbsTol. % % ORBITODE runs a demo of event location where the ability to % specify the direction of the zero crossing is critical. Both % the point of return to the initial point and the point of % maximum distance have the same event function value, and the % direction of the crossing is used to distinguish them. % % The orbit of the third body is plotted using the output function % ODEPHAS2. % % L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary % Differential Equations, W.H. Freeman & Co., 1975. % % See also ODE45, ODE23, ODE113, ODESET, ODEPHAS2, FUNCTION_HANDLE. % Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94 % Copyright 1984-2014 The MathWorks, Inc. fprintf('\nThis is an example of event location where the ability to\n'); fprintf('specify the direction of the zero crossing is critical. Both\n'); fprintf('the point of return to the initial point and the point of\n'); fprintf('maximum distance have the same event function value, and the\n'); fprintf('direction of the crossing is used to distinguish them.\n\n'); fprintf('Calling ODE45 with event functions active...\n\n'); fprintf('Note that the step sizes used by the integrator are NOT\n'); fprintf('determined by the location of the events, and the events are\n'); fprintf('still located accurately.\n\n'); % Problem parameters mu = 1 / 82.45; mustar = 1 - mu; y0 = [1.2; 0; 0; -1.04935750983031990726]; tspan = [0 7]; options = odeset('RelTol',1e-5,'AbsTol',1e-4,'OutputFcn',@odephas2,... 'Events',@events); figure; [t,y,te,ye,ie] = ode45(@f,tspan,y0,options); figure; plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o'); title('Restricted three body problem'); ylabel('y(t)'); xlabel('x(t)'); % ----------------------------------------------------------------------- % Nested functions -- problem parameters provided by the outer function. % function dydt = f(t,y) % Derivative function -- mu and mustar shared with the outer function. r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5; r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5; dydt = [ y(3) y(4) 2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - mu*((y(1)-mustar)/r23) -2*y(3) + y(2) - mustar*(y(2)/r13) - mu*(y(2)/r23) ]; end % ----------------------------------------------------------------------- function [value,isterminal,direction] = events(t,y) % Event function -- y0 shared with the outer function. % Locate the time when the object returns closest to the initial point y0 % and starts to move away, and stop integration. Also locate the time when % the object is farthest from the initial point y0 and starts to move closer. % % The current distance of the body is % % DSQ = (y(1)-y0(1))^2 + (y(2)-y0(2))^2 = <y(1:2)-y0(1:2),y(1:2)-y0(1:2)> % % A local minimum of DSQ occurs when d/dt DSQ crosses zero heading in % the positive direction. We can compute d/dt DSQ as % % d/dt DSQ = 2*(y(1:2)-y0)'*dy(1:2)/dt = 2*(y(1:2)-y0)'*y(3:4) % % y0 is shared with the outer function. dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4)); value = [dDSQdt; dDSQdt]; isterminal = [1; 0]; % stop at local minimum direction = [1; -1]; % [local minimum, local maximum] end % ----------------------------------------------------------------------- end % orbitode