Not able to run
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposte (1)
Sulaymon Eshkabilov
il 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 Commenti
Vedere anche
Categorie
Scopri di più su Array Geometries and Analysis in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!