problem figuring out 2 different solutions
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Eliraz Nahum
le 30 Sep 2018
Réponse apportée : Eliraz Nahum
le 30 Sep 2018
hello everyone,
my friend and I have to hand out our HW today, and can't understand why after going over and over again on our codes - we can't get identical plots. the main difference is with x values. I attached the definition of the physical problem as a photo.
In addition, I am attaching the first code:
clear all
close all
clc
t1=32.2;
w0=6169;
w11=2695;
t2=138.2;
w12=2082;
w21=397;
w22=100;
T1=26950;
T2=3973;
x0=0;
y0=2.0926*(10^7);
xdot0=34.7;
ydot0=196.9;
g=32.15;
GM=1.4077*(10^16);
t_delta=1;
counter=1;
time=0;
v_time=[];
x=[];
y=[];
xdot=[];
ydot=[];
x2dot=[];
y2dot=[];
x(counter)=x0;
y(counter)=y0;
xdot(counter)=xdot0;
ydot(counter)=ydot0;
division1=T1/w0;
x2dot(counter)=(-GM*(x(1)/((x(1)^2 + y(1)^2)^1.5)))+(g*division1)*(xdot(1)/((xdot(1)^2 + ydot(1)^2)^0.5));
y2dot(counter)=(-GM*(y(1)/((x(1)^2 + y(1)^2)^1.5)))+(g*division1)*(ydot(1)/((xdot(1)^2 + ydot(1)^2)^0.5));
r=sqrt(x0^2+y0^2);
size=r;
while (size>=r)
time=time+t_delta;
if time<t1
W=w0+(((w11-w0)/t1)*time);
T=T1;
elseif (time>=t1)&&(time<t2)
W=w12+(((w21-w12)/(t2-t1))*(time-t1));
T=T2;
else
W=w22;
T=0;
end
TvsW=T/W;
counter=counter+1;
x(counter)=x(counter-1)+xdot(counter-1)*t_delta;
y(counter)=y(counter-1)+ydot(counter-1)*t_delta;
xdot(counter)=xdot(counter-1)+x2dot(counter-1)*t_delta;
ydot(counter)=ydot(counter-1)+y2dot(counter-1)*t_delta;
x2dot(counter)=(-GM*(x(counter)/((x(counter)^2 + y(counter)^2)^1.5)))+(g*TvsW)*(xdot(counter)/((xdot(counter)^2 + ydot(counter)^2)^0.5));
y2dot(counter)=(-GM*(y(counter)/((x(counter)^2 + y(counter)^2)^1.5)))+(g*TvsW)*(ydot(counter)/((xdot(counter)^2 + ydot(counter)^2)^0.5));
v_time(counter)=time;
size=sqrt(x(counter)^2+y(counter)^2);
end
plot(x,y,'.b')
xlim([0 (10*10^6)])
ylim([0 (3*10^7)])
figure
plot(v_time,x,'.r')
hold on
plot(v_time,y,'.g')
figure
plot(v_time,xdot,'.r')
hold on
plot(v_time,ydot,'.g')
figure(1)
hold on
for counter=0:0.001*pi:2*pi
plot(r*cos(counter),r*sin(counter),'.m')
hold on
end
hold off
and here is the second code:
close all;
clear all;
clc;
%Basic data
GM = 1.4077 * 10^16;
v0 = 200;
gama = deg2rad(80);
g = 32.17;
dt=1;
Time=100;
%Starting conditions
x0 = 0;
y0 = 2.0926 * 10^7;
xdot0 = v0 * cos(gama);
ydot0 = v0 * sin(gama);
for n=1:Time
T = 26950;
w = 6169;
t=0;
s1(1,n) = 30.59 + rand * 2 * 1.61;
s2(1,n) = 131.29 + rand * 2* 6.91;
x_val(1,1) = x0;
x_dot_val(1,1) = xdot0; %=x1dot
y_val(1,1) = y0; % =R world
y_dot_val(1,1) = ydot0; %=y1dot
while sqrt(x_val ^ 2 + y_val ^ 2) >= y0
v = ( x_dot_val ^ 2 + y_dot_val ^ 2 ) ^0.5;
aT = ( g * T ) / w;
x_2dot_val = ( -GM * x_val) / ((x_val ^ 2 + y_val ^ 2)^1.5) + (aT * x_dot_val) / v;
y_2dot_val = ( -GM * y_val) / ((x_val ^ 2 + y_val ^ 2)^1.5) + (aT * y_dot_val) / v;
x_val = x_val + dt * x_dot_val;
x_dot_val = x_dot_val + dt * x_2dot_val;
y_val = y_val + dt * y_dot_val;
y_dot_val = y_dot_val + dt * y_2dot_val;
if t<=s1(1,n)
w = ((2695-6169)/(32.2-0))*dt + 6169;
elseif t>s2(1,n)
w = 100;
T = 0;
else
w = ((397 - 2082) / (138.2 - 32.2)) * dt + 2082;
T = 3973;
end
t=t+dt;
x_vector(1,t) = x_val;
y_vector(1,t) = y_val;
x_dot_vector(1,t) = x_val;
y_dot_vector(1,t) = y_val;
x_2dot_vector(1,t) = x_val;
y_2dot_vector(1,t) = y_val;
end
length_arc(1,n) = atan (x_val/y_val) * y0; %חישוב הקשת במעגל
figure(1);
plot(x_vector,y_vector);
hold on;
end
% g = viscircles(x,y,0,0,y0);
% plot(g);
th = 0:pi/50:2*pi;
xunit = y0 * cos(th) ;
yunit = y0 * sin(th);
h = plot(xunit, yunit);
hold off
figure(2);
hist(length_arc,20);
2 commentaires
madhan ravi
le 30 Sep 2018
"I attached the definition of the physical problem as a photo" you didn't
Réponse acceptée
Bruno Luong
le 30 Sep 2018
Modifié(e) : Bruno Luong
le 30 Sep 2018
if t<=s1(1,n) w = ((2695-6169)/(32.2-0))*dt + 6169; elseif t>s2(1,n) w = 100; T = 0; else w = ((397 - 2082) / (138.2 - 32.2)) * dt + 2082; T = 3973; end
The dt coded above is timestep, which is wrong if should be (t-s1) or (t-s0)
0 commentaires
Plus de réponses (1)
Voir également
Catégories
En savoir plus sur Logical dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!