problem figuring out 2 different solutions
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Eliraz Nahum
am 30 Sep. 2018
Beantwortet: Eliraz Nahum
am 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 Kommentare
madhan ravi
am 30 Sep. 2018
"I attached the definition of the physical problem as a photo" you didn't
Akzeptierte Antwort
Bruno Luong
am 30 Sep. 2018
Bearbeitet: Bruno Luong
am 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 Kommentare
Weitere Antworten (1)
Siehe auch
Kategorien
Mehr zu Text Files finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!