problem figuring out 2 different solutions

1 visualizzazione (ultimi 30 giorni)
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 Commenti
madhan ravi
madhan ravi il 30 Set 2018
"I attached the definition of the physical problem as a photo" you didn't
Eliraz Nahum
Eliraz Nahum il 30 Set 2018
thanks for informing me. I attached now.

Accedi per commentare.

Risposta accettata

Bruno Luong
Bruno Luong il 30 Set 2018
Modificato: Bruno Luong il 30 Set 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)

Più risposte (1)

Eliraz Nahum
Eliraz Nahum il 30 Set 2018
thank you very much!!! you have no idea how you helped us.

Prodotti


Release

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by