Not able to run
1 view (last 30 days)
Show older comments
y0 = [r0 Vn0];
t0 = 0;
nrevs = 100;
tmax = nrevs*Period;
h0 = 0;
x0 = [r0; h0; V0];
r= @(x) [sqrt(x(1)^2 + x(2)^2 + x(3)^2)];
f = @(t,x) [x(3); x(2)*x(3)^2/x(1) -mu/x(1)^2; -x(2)*x(3)/x(1)];
settings1=odeset('RelTol',1e-4,'AbsTol',1e-4);
settings2=odeset('RelTol',1e-8,'AbsTol',1e-8);
settings3=odeset('RelTol',1e-11,'AbsTol',1e-11);
[T1,x1]=ode45(f,[0 tmax],x0,settings1)
[T2,x2]=ode45(f,[0 tmax],x0,settings2)
[T3,x3]=ode45(f,[0 tmax],x0,settings3)
% Plot the trajectories for each tolerance
% Plot the trajectories for each tolerance
figure;
plot(x1(:,1).*cos(x1(:,2)), x1(:,1).*sin(y1(:,2)), '-', ...
x2(:,1).*cos(x2(:,2)), x2(:,1).*sin(x2(:,2)), '-', ...
x3(:,1).*cos(x3(:,2)), x3(:,1).*sin(x3(:,2)), '-');
legend('tol=1e-4', 'tol=1e-8', 'tol=1e-11');
title('Orbit Trajectories');
0 Comments
Answers (1)
Sulaymon Eshkabilov
on 20 Feb 2023
Here is the corrected code (Note that the Initial Conditions were missing and some arbitrary values are taken. Therefore they need to be changed. Also, note that in figure 1 three solution sets overlap each other. Hence, they are plotted in separate) :
t0 = 0;
nrevs = 100;
Period = 2;
tmax = nrevs*Period;
mu = 2.5;
h0 = 0; %
r0 = 1; % use the given value
V0 = 0.5; % use the given value
x0 = [r0; h0; V0];
f = @(t,x) ([x(3); x(2).*x(3).^2./x(1)-mu./x(1).^2; -x(2).*x(3)./x(1)]);
settings1=odeset('RelTol',1e-4,'AbsTol',1e-4);
settings2=odeset('RelTol',1e-8,'AbsTol',1e-8);
settings3=odeset('RelTol',1e-11,'AbsTol',1e-11);
[T1,y1]=ode45(f,[0 tmax],x0,settings1);
[T2,y2]=ode45(f,[0 tmax],x0,settings2);
[T3,y3]=ode45(f,[0 tmax],x0,settings3);
% Plot the trajectories for each tolerance
figure;
plot(y1(:,1).*cos(y1(:,2)), y1(:,1).*sin(y1(:,2)), 'r-', ...
y2(:,1).*cos(y2(:,2)), y2(:,1).*sin(y2(:,2)), 'b--', ...
y3(:,1).*cos(y3(:,2)), y3(:,1).*sin(y3(:,2)), 'g-.');
legend('tol=1e-4', 'tol=1e-8', 'tol=1e-11');
title('Orbit Trajectories');
axis equal
figure
plot(y1(:,1).*cos(y1(:,2)), y1(:,1).*sin(y1(:,2)), 'r-')
legend('tol=1e-4');
title('Orbit Trajectories');
axis equal
figure;
plot(y2(:,1).*cos(y2(:,2)), y2(:,1).*sin(y2(:,2)), 'b--');
legend('tol=1e-8');
title('Orbit Trajectories');
axis equal
figure;
plot(y3(:,1).*cos(y3(:,2)), y3(:,1).*sin(y3(:,2)), 'g-.');
legend('tol=1e-11');
title('Orbit Trajectories');
axis equal
0 Comments
See Also
Categories
Find more on Array Geometries and Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!