ODE Event Function not Properly Configured

2 次查看（过去 30 天）
Mark Munson
Mark Munson 2022-8-18
编辑： Torsten 2022-8-24
I am trying to execute two event functions "reaching L2" and "going beyond L3". But everytime I get this error:
Error using odeevents
The ODE option 'Events' must be set to 'on' | 'off'
Error in ode45 (line 140)
odeevents(odeIsFuncHandle,ode,t0,y0,options,varargin);
Error in Lyapunov_Manifolds_Main (line 195)
[t1u, x1u, te1u, xe1u, ie1u] = ode45('CR3BP', 4*t1span, ..
-----------------------------------------------------------------------------------------------------------------
%% Clear MATLAB Workspace.
clear
clc
close all
format compact
format long
%% Constants.
mu = 0.0121505856; % Mass parameter of the Earth-Moon system.
n = 1; % Mean motion of the Moon's orbit around Earth.
L2x = 1.155682165407869; % x-coordinate for the L2 Lagrange point.
% Set the time periods of L1 and L2 Lyapunov orbits.
T1 = 2.7953;
T2 = 3.4155;
% Set the initial STM.
phi0 = eye(6,6);
%% Integration setup.
% Input setup.
s1_0 = [ 0.8189;
0;
0;
0;
0.17454;
0;
reshape(phi0, 36, 1)];
s2_0 = [ 1.1809;
0;
0;
0;
-0.15587;
0;
reshape(phi0, 36, 1)];
t1span = [0 T1];
t2span = [0 T2];
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
options1 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', 'reaching_L2');
options2 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', 'going_beyond_L3');
%% Integrate the L_1 and L_2 Lyapunov orbits.
[~, x1] = ode45('CR3BP_STM', t1span, s1_0, options);
[~, x2] = ode45('CR3BP_STM', t2span, s2_0, options);
%% Start plotting the orbits in a figure, later to overlay with the
% manifolds in this same figure.
plot3(x1(:,1), x1(:,2), x1(:,3), 'r', 'DisplayName', 'L_1 Lyapunov Orbit',...
'LineWidth', 3);
hold on;
plot3(x2(:,1), x2(:,2), x2(:,3), 'r', 'LineWidth', 3, ...
'DisplayName', 'L_2 Lyapunov Orbit');
scatter3(1 - mu, 0, 0, 'b*', 'LineWidth', 4, 'DisplayName', 'Moon');
scatter3(-mu, 0, 0, 'b*', 'LineWidth', 10, 'DisplayName', 'Earth');
grid on;
axis equal;
xlabel('x [NDU]');
ylabel('y [NDU]');
zlabel('z [NDU]');
%% Select 100 points along each of the L1 and L2 Lyapunov orbits
% respectively.
points = 100;
x1Points = zeros(points, 6);
x2Points = zeros(points,6);
step1 = floor(size(x1, 1)/points);
step2 = floor(size(x2, 1)/points);
for ii = 1:points
x1Points(ii,:) = x1((ii - 1)*step1 + 1, 1:6);
x2Points(ii,:) = x2((ii - 1)*step2 + 1, 1:6);
end
%% Initialize plotting control variables for stable and unstable manifolds
% of the L1 and L2 Lyapunov orbits.
magenta = 0;
yellow = 0;
red = 0;
cyan = 0;
%% Iterate over the 100 points along the L1 and L2 Lyapunov orbits.
% First 100 iterations go over in the positive direction of perturbations
% for both stable and unstable manifolds for L1 orbit and in the negative
% direction of perturbations for both stable and unstable manifolds for L2
% orbit.
% Second 100 iterations go over in the negative direction of perturbations
% for both stable and unstable manifolds for L1 orbit and in the positive
% direction of perturbations for both stable and unstable manifolds for L2
% orbit.
% In both the 100 iterations, same 100 points along the original orbits are
% used.
% L_1.
for ii = 1:2*points
% Set the correct index for initial points to perturb (especially
% when ii greater than 100.
k = ii;
if ii > points
k = ii - points;
end
end
%% Integration setup.
s1_0 = [x1Points(k,:)'; reshape(phi0, 36, 1)];
s2_0 = [x2Points(k,:)'; reshape(phi0, 36, 1)];
%% Integrate the Lyapunov L1 and L2 orbits.
[t1, x1] = ode45('CR3BP_STM', t1span, s1_0, options);
[t2, x2] = ode45('CR3BP_STM', t2span, s2_0, options);
%% Compute the monodromy matrix for L1 and L2 orbits.
monodromy1 = reshape(x1(end, 7:42), 6, 6);
monodromy2 = reshape(x2(end, 7:42), 6, 6);
%% Compute the eigenvectors and digonal matrices of eigenvalues for L1 and
%% L2 orbits.
[eigVectors1, eigVDiag1] = eig(monodromy1);
[eigVectors2, eigVDiag2] = eig(monodromy2);
%% Extract eigenvalues from the diagonal matrices of eigenvalues.
eigValues1 = diag(eigVDiag1);
eigValues2 = diag(eigVDiag2);
%% Extract the indices of the unstable and stable eigenvalues.
[~, uIndex1] = max(eigValues1);
[~, uIndex2] = max(eigValues2);
[~, sIndex1] = min(eigValues1);
[~, sIndex2] = min(eigValues2);
%% Set the perturbation value, approximately 38 kilometers.
d = 0.0001;
%% Normalize the unstable and stable eigenvectors using the position
% components of the corresponding eigenvectors.
normUEigVec1 = eigVectors1(:, uIndex1)/norm(eigVectors1(1:3, uIndex1));
normUEigVec2 = eigVectors2(:, uIndex2)/norm(eigVectors2(1:3, uIndex2));
normSEigVec1 = eigVectors1(:, sIndex1)/norm(eigVectors1(1:3, sIndex1));
normSEigVec2 = eigVectors1(:, sIndex1)/norm(eigVectors1(1:3, sIndex2));
%% Introduce perturbations into the initial conditions using the unstable
% and stable eigenvectors and perturbation value.
if ii < points + 1
x1InitUnstable = x1Points(ii, :)' + d*normUEigVec1;
x2InitUnstable = x2Points(ii, :)' - d*normUEigVec2;
x1InitStable = x1Points(ii, :)' + d*normSEigVec1;
x2InitStable = x2Points(ii, :)' - d*normSEigVec2;
else
k = ii - points;
x1InitUnstable = x1Points(k, :)' - d*normUEigVec1;
x2InitUnstable = x2Points(k, :)' + d*normUEigVec2;
x1InitStable = x1Points(k, :)' - d*normSEigVec1;
x2InitStable = x2Points(k, :)' + d*normSEigVec2;
end
%% Integrate the perturbed initial conditions forward to get the unstable
% manifolds.
[t1u, x1u, te1u, xe1u, ie1u] = ode45('CR3BP', 4*t1span, ...
x1InitUnstable, options1);
[t2u, x2u, te2u, xe2u, ie2u] = ode45('CR3BP', 4*t1span, ...
x2InitUnstable, options2);
%% For unstable L1 manifolds, plot only those manifolds that do not go to
% L2 point, but traverse towards the Earth and the interior region, this
% is for Condition - 1 from the problem statement.
if size(ie1u, 1) == 0
p1u = plot3(x1u(:,1), x1u(:,2), x1u(:,3), 'm', 'DisplayName', ...
'L_1 Lyapunov Unstable');
% Conditions for suppressing legend entries in the manifolds plot for
% every manifold except the first one.
if magenta == 1
set(get(get(p1u, 'Annotation'), 'LegendInformation'),...
'IconDisplayStyle', 'off');
else
magenta = 1;
end
end
% For Condition - 2 of the problem statement, the following logic is used to
% enforce plotting only those unstable manifolds that go beyond Moon, L3
% and cross the x-axis (y=0).
indicesUnstable1 = find(ie2u == 1);
indicesUnstable2 = find(ie2u == 2);
indicesUnstable3 = find(ie2u == 3);
startUnstable = 1;
if size(indicesUnstable2, 1) > 0
startUnstable = indicesUnstable2(1);
end
indicesUnstable3 = find(ie2u(startUnstable:end) == 3);
if size(indicesUnstable2, 1) > 0
indexUnstable = find(abs(t2u - te2u(indicesUnstable2(1))) < 0.005);
if size(indicesUnstable3, 1)>0
indexUnstable = find(abs(t2u - te2u(startUnstable-1+...
indicesUnstable3(1))) < 0.008);
end
p2u = plot3(x2u(1 : indexUnstable, 1), x2u(1 : indexUnstable, 2), ...
x2u(1 : indexUnstable,3), 'r', 'DisplayName', 'L_2 Lyapunov Unstable');
% Conditions for suppressing legend entries in the manifolds plot for
% every manifold except the first one.
if red == 1
set(get(get(p2u, 'Annotation'), 'LegendInformation'),...
'IconDisplayStyle', 'off');
else
red = 1;
end
end
% Integrate the perturbed initial conditions backwards for stable
% manifolds.
[t1s, x1s, te1s, xe1s, ie1s] = ode45('CR3BP', -4*t1span, ...
x1InitStable, options1);
[t2s, x2s, te2s, xe2s, ie2s] = ode45('CR3BP', -4*t2span, ...
x2InitStable, options2);
% For stable L1 manifolds, plot only those manifolds that do not go to
% L2 point, but traverse towards the Earth and the interior region, this
% is for Condition - 1 from the problem statement.
if size(ie1s, 1) == 0
p1s = plot3(x1s(:,1), x1s(:,2) ,x1s(:,3), 'y', 'DisplayName', ...
'L_1 Lyapunov Stable');
% Conditions for suppressing legend entries in the manifolds plot for
% every manifold except the first one.
if yellow == 1
set(get(get(p1s, 'Annotation'), 'LegendInformation'),...
'IconDisplayStyle', 'off');
else
yellow = 1;
end
end
% For Condition - 2 of the problem statement, the following logic is used to
% enforce plotting only those stable manifolds that go beyond Moon, L3
% and cross the x-axis (y=0).
indicesStable1 = find(ie2s == 1);
indicesStable2 = find(ie2s == 2);
start = 1;
if size(indicesStable2, 1) > 0
start = indicesStable2(1);
end
indicesStable3 = find(ie2s(start : end) == 3);
if size(indicesStable2, 1) > 0
indexStable = find(abs(t2s - te2s(indicesStable2(1))) < 0.005);
if size(indicesStable3, 1)>0
indexStable=find(abs(t2s - te2s(start - 1 + indicesStable3(1))) < 0.005);
end
p2s = plot3(x2s(1 : indexStable, 1),x2s(1 : indexStable, 2),...
x2s(1 : indexStable,3), 'c', 'DisplayName', 'L_2 Lyapunov Stable');
% Conditions for suppressing legend entries in the manifolds plot for
% every manifold except the first one.
if cyan == 1
set(get(get(p2s, 'Annotation'), 'LegendInformation'),...
'IconDisplayStyle', 'off');
else
cyan = 1;
end
end
legend('Location', 'northeast', 'FontWeight', 'bold');
fig.Color = 'white';

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-8-18
Always work with function handles instead of strings:
options1 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', @reaching_L2);
options2 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', @going_beyond_L3);
[~, x1] = ode45(@CR3BP_STM, t1span, s1_0, options);
[~, x2] = ode45(@CR3BP_STM, t2span, s2_0, options);
[t1, x1] = ode45(@CR3BP_STM, t1span, s1_0, options);
[t2, x2] = ode45(@CR3BP_STM, t2span, s2_0, options);
[t1u, x1u, te1u, xe1u, ie1u] = ode45(@CR3BP, 4*t1span, ...
x1InitUnstable, options1);
[t2u, x2u, te2u, xe2u, ie2u] = ode45(@CR3BP, 4*t1span, ...
x2InitUnstable, options2);
[t1s, x1s, te1s, xe1s, ie1s] = ode45(@CR3BP, -4*t1span, ...
x1InitStable, options1);
[t2s, x2s, te2s, xe2s, ie2s] = ode45(@CR3BP, -4*t2span, ...
x2InitStable, options2);
If this does not solve your problem, you will have to include the functions to run your code.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Computations 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by