MATLAB: Particle Velocity. Trying to solve for v0 with loops/conditional statements

loopsparticle motion

Disclamier: I am really bad at Matlab. Please go easy on me 🙂
Problem: I am given a final distance (8.9m), and I want to solve for initial velocity
Approach: I want to use loops/conditionals to solve for the problem. I know it's not the most elegant approach, but this problem is driving me crazy. I would like to increment velocity during each while loop, then check the final value of distance against my given distance (8.9m). If my distance from the while loop is too low, then I would like to increment v0 by 0.01, and run the loop again. When the distance from the while loop equals my given distance, then I have my v0.
I have a good working while loop to figure out the distance from a given v0; however, I am having a hard time incrementing v0 in a way that doesn't put my computer in an infinite loop or returns poor results.
Here's my code:
% Prepare Workspace
clear;      % Clears variables from the workspace
clc;        % Clears the command window
clf;        % Clears current figure window
% Define variables and constraints
g=9.81;  %gravity
m=80;  % mass
Cd= 0.72;  %coefficent of drag
A = 0.5;   % cross-sectional area
theta=pi/8;  %22.5 deg
v=11.1;  %inital velocity
dt=0.0001;  %change in time / time step
% Control Jump 1 Data
rho1=0.94;  %kg/m^3
D1=.5*Cd*rho1*A;
x1=zeros(1,10000);  %x axis index Jump 1

y1=zeros(1,10000);  %y axis index Jump 1

t1=zeros(1,10000);  %time index Jump 1

    xdot=v*cos(theta);  % problem description

    ydot=v*sin(theta);  % problem description
%While loop Control Jump 1  rho 0.94
i=1;
L1=0;
if L1<8.9
    v=v+0.01;
     x1=zeros(1,10000);  %x axis index Jump 1
    y1=zeros(1,10000);  %y axis index Jump 1
    t1=zeros(1,10000);  %time index Jump 1
while min(y1)>-0.001
    ax1= -(D1/m)*(xdot^2+ydot^2)^0.5*xdot;
    ay1= -g-(D1/m)*(xdot^2+ydot^2)^0.5*ydot;
    xdot=xdot+ax1*dt;
    ydot=ydot+ay1*dt;
    x1(i+1)=x1(i)+xdot*dt+.5*ax1*dt^2;
    y1(i+1)=y1(i)+ydot*dt+.5*ay1*dt^2;
    t1(i+1)=t1(i)+dt;
    i=i+1;
end
    
    L1=max(x1);  %Length of Jump 1
    T1=max(t1);  %Time of Jump 1
    Y1=max(y1);  %Vertical Height of Jump 1
   
i=1;
end
Thanks for any help!!

Best Answer

There is nothing wrong with the approach that you want to use, however there are some details that will help improve the efficiency of the method.
First, I the really smart Matlab people are always reminding me that the variable i is a defined Matlab constant, and is used in complex math notation, therefore, using i as a loop index can be considered a bad habit. Your code will work using i as the counting variable, or loop index, but just realize that you are ghosting a Matlab constant.
Next, it is my opinion that a counting loop should always be structured such that when the loop ends, the counting variable contains the number of loop itterations (and the number of elements in the accumulated vectors). To accomplish this, the counter should be incremented at the start of the loop, and therefore it needs to be initialized to zero outside the loop; i.e.:
i=0
while  (..)
   i=i+1;  % this is the first statement inside the loop
   ...
   ...
end
Now when the loop is completed, i contains the number of elements in x1, y1, and t1.
Next, your algorithm will be very sensitive to the step size because you are using the last loop itteration as the final value for the distance. This means that there will be an error in the final value that is a function of the step size, and the only way to minimize this error is to make the step very small. A better aproach is to use linear interpolation between the last and the next-to-last step to fine the precise point where the vertical trajectory crossed the zero boundary (i.e. the ground plane). Using interpolation will allow you to get a more precise solution while using a much larger step size.
To do this, you need to change the termination condition on the while lop so that it will always terminate on the first step after the Y position passes zero (goes negative). You can do this using
while y1(i) >= 0
    
end
(edit: on the first pass, y1(i) will have a subscript of zero, and this will give an error. Need to change this to avoid the error)
h = 0
 while h>=0
     i=i+1;
      ...
      ...
     h=y1(i);
 end
This will allow you to initialize y1 at zero and the loop will not terminate until y1 goes negative.
After the loop terminates, you want to replace the last value in x1,y1,and t1 with the interpolated value.
First, calculate the interpolation fraction (remember i contains the number of points in the trajectory, and the last point is below the ground):
k = -y1(i-1) / (y1(i) - y1(i-1));
Now calculate the interpolated values:
  xinterp = x(i-1) + k * (x(i) - x(i-1));
  yinterp = 0;
  tinterp = t(i-1) + k * (t(i) - t(i-1));
Now replace the end values with the iterpolated ones:
 x1(i) = xinterp;
 y1(i) = yinterp;
 t1(i) = tinterp;
Now the last point in the trajectory is exactly on the ground.
the length of the jump = x1(i)
the time of the jump = t1(i)
and the height of the jump = max(y1)
After making this change, you should be able to get a pretty good answer using a much biger step size, and this will speed up your code.
Finally, in order to locate the initial velocity for a given distance, use an interpolation method to calculate the next guess and you will converge on the answer much faster than simply taking small steps in the initial velocity.
Related Question