MATLAB: How to solve first order ODEs to find four components (for velocity and displacement of projectile in X and Y directions)

Question

Hello dear experts in Matlab,

I did a Matlab code to simulate the motion of a projectile in air.

I hope to find TEN reasonable corresponging values (vectors) for the four components (X-velocity, Y-velocity, X-displacement, and Y-displacement) against the Time period (t=[0:0.5:4.5]).

The code that I did fails in executing reasonable values. I don't know what is wrong with it OR with the Four Functions it calls to run.

To simplify reading my codes that I ATTACHED HERE;

xcomp1 file is the functiond of x-velocity (v_x)

xcomp11 file is the functiond of x-displacement (s_x)

ycomp2 file is the functiond of y-velocity (v_y)

ycomp22 file is the functiond of y-displacement (s_y)

My deep thanks for each notice and support!

Answer 1

A little more like the following perhaps (it can all be done in one file):

%Primary Conditions 
tspan=[0:0.5:4.5];     %Period of Time that we will  calculate through it.
% You need an initial velocity and an initial angle or there will be no 
% variation in the y-direction!!
% For example
v0 = 10; theta0 = pi/6;
initial_vx=v0*cos(theta0);      %initial value for v_x  
initial_vy=v0*sin(theta0);       %initial value for v_y
initial_x11=0;     %initial value for s_x
initial_y22=0;     %initial value for s_y
xyv0 = [initial_x11, initial_y22, initial_vx, initial_vy];
[t,xyv]=ode45(@rate,tspan,xyv0);
s_x = xyv(:,1);
s_y = xyv(:,2);
v_x = xyv(:,3);
v_y = xyv(:,4);
plot(t,v_x)  
xlabel('Time (S)') 
ylabel('Vx (m/s)') 
figure 
plot(t,v_y)
xlabel('Time (S)') 
ylabel('Vy (m/s)')
figure
plot(s_x,s_y)
xlabel('Sx (m)') 
ylabel('Sy (m))')
function dxyvdt = rate(t,xyv)
%         s_x = xyv(1);   % Not needed inside this function
%         s_y = xyv(2);   %       ditto
        v_x = xyv(3);
        v_y = xyv(4);
        rho_p=0.9986*10^3;      %density of object
        rho_f=1.2077;           %density of air 
        v=sqrt(v_x^2+v_y^2);    %object velocity
%         T_p0=33;                %initial droplet temperature
%         T_f=18;                 %air temperature 
%         R_H=0.5;                %Relative Humidity
        k1=74*10^-12;           %evaporation constant 
        R0=100*10^-6;           %initial object radius
        Reynolds_no=672;  %%%%% Why isn't Reynolds number a function of velocity?%%%%
        c_d=(21.12/Reynolds_no)+(6.3/sqrt(Reynolds_no))+0.25;  %drag coefficient
        
        % The drag force should be applied to the velocity v and then resolved
        % into x and y components.  You shouldn't use the resolved velocities 
        % to calculate separate drag forces (though, given that your drag 
        % is linear in velocity here, that probably won't matter too much!).
        % 
        drag = -3*c_d*rho_f*v;
        theta = atan(v_y/v_x);
        drag_x = drag*cos(theta);
        drag_y = drag*sin(theta);
        
        dxyvdt = [v_x;
                  v_y;
                  drag_x/(8*sqrt(R0^2-2*k1*t)*rho_p);
                  drag_y/(8*sqrt(R0^2-2*k1*t)*rho_p)];
end

NB If you change your timespan to , say, tspan = 0:0.1:4.5; you will get smoother curves for v_x and v_y.