problem figuring out 2 different solutions

Question

Eliraz Nahum 2018-9-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/421469-problem-figuring-out-2-different-solutions

回答： Eliraz Nahum 2018-9-30

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 个评论
显示无隐藏无

madhan ravi 2018-9-30

"I attached the definition of the physical problem as a photo" you didn't

Eliraz Nahum 2018-9-30

thanks for informing me. I attached now.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bruno Luong 2018-9-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/421469-problem-figuring-out-2-different-solutions#answer_339097

编辑：Bruno Luong 2018-9-30

                 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)