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

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!!

Answer 1

[Jim Riggs]

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.

[threex5x7]

Jim,

Thanks for the detailed answer.

I am still confused.

I have a vauge understanding of the interpolation (Newton method?), my goal is to be able to increment velocity only.

How would you go about that the 'ugly' way? I would like to refine and optimize this code after I can be sure my numbers are good.

My desired approach is:

1. run while loop with current value of v in the code
2. compare distance obtained from while loop above with given distance (8.9m).
3. increment v
4. re-run while loop
5. compare distance obtained from while loop with given distance (8.9m).
6. end loop when v corresponds to given distance (8.9m).

in the steps above, I am confident that step 1 above runs correctly.

thank you very much for the help, and sorry that I'm having problems following along with your help.

Thanks again for your patience and understanding!

threex5x7

[Jim Riggs]

1) Run loop to calculate distance L1 given initial velocity v1.

2) Compare to desired value: Lerror = Ldesired - L1 (Positive error means we need more V, negative error means we need less V)

Vloopcount = 1

3) For the first correction (Vloopcount=1), simply pick an increment for V, say deltaV = 0.5*sign(Lerror)

set Lold = L1, Vold = V

V = V + deltaV

4) Re-run the trajectory solution to get a new L1, based on the new V.

5) Vloopcount = Vloopcount + 1

Lerror = Ldesired - L1.

if abs(Lerror) > (Desired Accuracy) then calculate the update to V using interpolation:

K = (Ldesired - Lold ) / (L1 - Lold)

deltaV = K * (V - Vold)

Lold = L1

Vold = V

V = V + deltaV

(go back to 4)

Keep doing 4) and 5) until the error is sufficiently small.

Vloopcount = how many iterations it takes to converge on the answer. It also helps you do something different on the first update (i.e. when Vloopcount =1)

[threex5x7]

Thanks again Jim for the help!

I'm going to work with this, and see how it goes.

I'll reply back soon.....

[Jim Riggs]

LINEAR INTERPOLATION

Suppose I have two measured coordinates, p1 and p2 (i know x1,y1 and x2,y2). In your case, this would be two solutions where you have a velocity and a distance solution, v1,L1 and v2, L2 If x is the value that I know (the desired distance) and Y is the value that I want t solve for (the initial velocity) we can use linear interpolation to find the value.

In the figure, X represents the value that I want (Ldesired) , and you can see that the ratio (x-x1) / (x2-x1) is the fraction of the distance that the desired point is between x1 and x2. (this is k, above). Let's assume that k has a value of .75 We can find the Y value that corresponds to X by going .75 of the distance from Y1 to Y2, Y = Y1 + .75(Y2 - Y1) .

This formula works even if X is not between X1 and X2. If X is to the left of X1, then k is negative. If X is to the right of X2, then k is greater than 1.

[threex5x7]

Thank you Jim for all the help!!

Makes alot more sense now visualizing linear interpolation.

I appreciate you helping me with this!!

Robert