How would you plot a graph which a ball then rolls down (say a y=x^2 graph)

Question

Sena Clarke 2019-8-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/477104-how-would-you-plot-a-graph-which-a-ball-then-rolls-down-say-a-y-x-2-graph

评论： Jack McCarthy 2022-5-12

11 个评论
显示 9更早的评论隐藏 9更早的评论

Kaushik Lakshminarasimhan 2019-8-26

编辑：Kaushik Lakshminarasimhan 2019-8-26

This sounds like physics homework, so you should try your luck at https://physics.stackexchange.com/ first before worrying about how to implement the solution in MATLAB

But here's a generous hint. The velocity of the ball would depend on where you release it so start by assuming an initial condition (x,y position where you release the ball). Now if you draw the free-body diagram, you'll see that at any given moment, the ball should accelerate along the tangent to curve y=x^2 at g*cos(theta) where theta = arctan(1/2x) and g is 10m/s^2.

Sena Clarke 2019-8-27

编辑：Sena Clarke 2019-8-27

Ha, no its not homework, but I could see why you would think that. No, I was just trying to program brachistochrone curves and doing a write-up on the times with different curves, like straight lines vs y=x^2 and y=sqrt(x) etc... with a ball going down and I don't know how to program a model.

请先登录，再进行评论。

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

Answer 1

David K. 2019-8-27

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/477104-how-would-you-plot-a-graph-which-a-ball-then-rolls-down-say-a-y-x-2-graph#answer_389272

ballOnSlope.m

I attached a .m file I created that should work pretty well. It uses a numerical solution by calculating the slope at discrete time steps and determining the acceleration. I only considered a point mass so I ignored rotation. It can also take both static and kinetic friction into account (im slightly worried I messed up here but I think it is right). The function takes in a function handle for what you want the slope to be. I made it so that the part of the function you want the ball on needs to be between x = 0 and 10 but that should be easily changed by changing the function itself.

I did not get around to allowing arguments to be dropped from the function or personally checking arguments and throwing errors.

Let me know how it works for you!

5 个评论
显示 3更早的评论隐藏 3更早的评论

Jim Riggs 2019-8-27

编辑：Jim Riggs 2019-8-27

I see a few problems with your posted answer.

First, the the change in the X position is not correct;

xloc = xloc + velocity * tstep;

The velocity is along the tangent to the curve, but this formula assumes that it is parallel to x.

Second, the friction should be proportional to the force (or acceleration) that is normal to the velociry, not proportional to the velocity.

Also, i recommend using central difference formulas for calculating the derivatives, rather than the backward difference. These are more accurate:

The total normal force (acceleration) is composed of the vector sum of the kinematic acceleration, plus the weight (mg).

where Cf is the coefficient of friction (either static or sliding).

The sliding friction is usually computed as a positive scalar value, and applied in the oposite direction from the velocity. This way, you avoid any sign errors in the friction contribution.

Comments on the Normal Force/Acceleration:

If you are comfortable using the kinematic equations for the acceleration of a point mass, note that the normal force (in this case) is the centrifugal force. It is the reaction force due to the kinematic centripetal acceleration. It, therefore, has the opposite sign from the centripetal acceleration.

David K. 2019-8-28

编辑：David K. 2019-8-28

在 MATLAB Online 中打开

Before I get to my response - I turned my function into a gui and attached it because they are fun and easy to use/make. You will need all 3 files named as is in the same folder to work.

The velocity I am using is parallel to the x axis. I calculated the acceleration tangent to the curve and then multiplied by cosine so that I only have the part of the acceleration that is parallel to the x axis. Then I assume the y position is on the ramp at the given x value.

I chose to model this as a point mass on a ramp problem with different angles so I am ignoring radius of curvatures and all the rotational forces.

Yeah I probably should have used central, that was just a hold over from when I was first throwing stuff together and it was more convenient to do backward. But I figured with a small enough dx it doesn't really matter too much. I made that change in the gui I attached.

And for friction, since I am calculating normal force with a rotated axis as if in the instance it is on a ramp and not a curve my acceleration is simply this:

N = mgcos(theta) 
ma = mgsin(theta)-coefF*mgcos(theta)
a = g(sin(theta)-coefF*cos(theta)) % acceleration tangent to curve
ax = g(sin(theta)-coefF*cos(theta))*cos(theta) % accel in x axis

However, apparently OP didn't even need friction so by setting coefF to 0 it kind of makes it a moot point anyways.

The way that I at least tried to make it always opposite of velocity was by multiplying it by the signum of velocity (outputs 1 or -1). I think that is what made you think I made it proportional to velocity. I checked my signs for +vel,+ramp; +vel,-ramp; -vel,+ramp; and -vel,-ramp and it appears to all check out with the way I am doing it.

TLDR it looks to me most of the problems with it were either misunderstandings of what I did or a disagreement in the way it should be modeled. If you wish to change it so that instead of being instantaneous ramps it takes into account that it is on a curve, feel free to change it and post another version.

Jim Riggs 2019-8-29

编辑：Jim Riggs 2019-8-29

OK. I see what you are doing. This might work for the frictionless case. (I regret to say that, for security reasons, I do not have Matlab on this machine, so I cannot download and run your code. I would like to try it out, but my Matlab is on an isolated network)

On a side note, I would sugest that you get in the habit of using vector forms whenever possible, avoiding trig forms, as vector calculations are much more robust, and do not have issues when changing quadrants. For example, I once worked on a missile program where the guidance law was coded using a lot of trig expressions. A navigation error in early flight tests caused the missile to believe that it went below the ground level while in terminal homing on the target, resulting in a negative altitude calculation. This caused an unexpected sign change in the guidance law, and the missule started flying away from the target, missing by a significant distance. This error would have been avoided using vector equations, and it would have hit the target even with the navigation error. You might argue that it's reasonable that missiles should't be required to work under the ground, but in the real world, these things happen, and it's always best to hit the target when possible.

For the problem at hand, I would define a unit tangent vector for the curve:

Ut = [Utx, Uty] = [dx, dy] / mag[dx, dy] (where dy is obtained from the numerical first derivative formula).

Now the unit normal vector Un = [Unx, Uny] = [ -Uty, Utx]

Use these two unit vectors and the dot product to obtain the tangent and normal quantities. No trig function required, and it works for all directions with no singularities.

This will also allow you to apply energy conservation principles because you are retaining the full velocity vector. In a frictionless environment the system energy remains constant, and the total velocity should be a function of the Y position, as potential energy is converted to kinetic energy, and vice-versa. You can use this relationship to test the accuracy of your solution, but only if you calculate the total velocity, and not just the x component.

So, I assure you that if you use vectors in place if trig functions, you are more likely to hit the target :)

Oh, and for the record, it does make a difference using central vs. fwd/bkwd difference for the slope calculation. Try it and see how it compares. Use conservation of energy as the accuracy metric and see how the step size effects it.

David K. 2019-8-29

在 MATLAB Online 中打开

Welp, I went about trying to calculate frictionless case for checking the velocity. My velocity function based on conservation of energy was

mghmax-mgh = 1/2 mv^2
v = sqrt(2*g*(hmax-h));

Then I tried to calculate the velocity in the y direction and it went pretty terribly. I couldn't do the acceleration in the same way so I used the difference in y position. And the velocities were super off the expected values at some places. It seems that xvelocity was pretty close if I use the expected x part of the expected velocity but the y part is very off.

I guess to make this a decent answer I will need to change it at some point. Ug, I really thought modeling each time step as a slope would work.

btw, the only reason I would need to be in the habit of vector forms is for making better answers here since I do not do this stuff for work xD.

Jim Riggs 2019-8-29

编辑：Jim Riggs 2019-8-29

Yes, it's very important to have a way to check yourself. It sounds like you understand the problem and you are on the right track.

Another way to test your ramp assumption/approximation is to actually use a linear function for the curve, and see if it works there. If not, there is some implementation issue.

请先登录，再进行评论。

Answer 2

Jim Riggs 2019-8-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/477104-how-would-you-plot-a-graph-which-a-ball-then-rolls-down-say-a-y-x-2-graph#answer_389672

编辑：Jim Riggs 2019-8-29

Here is a model for the kinematics.

clear data  % I use data to save values in the time loop
func = @(x) x.^2;  % this is the user function
% set model parameters and initial conditions
dx = 0.001;  % step used to compute numerical derivatives
dt = .01;    % integration time step
tstart = 0;    % starting value of time
tstop = 12;  % time to stop
x = 0.1;     % starting value of X
y = func(x)  % starting value of Y
speed = 40;  % initial speed
grav = 9.806;   % acceleration due to gravity
G = [0, -grav]; % gravity vector
nstep = (tstop-tstart)/dt;  % numer of calculation steps
cnt = 1;  % itteration counter
% initial energy state (per unit mass)
Ep = gravity*y;   % potential energy
Ek = 0.5*speed^2; % kinetic energy
Etot = Ep + Ek;   % total system energy
% initialize saved data table
% (you can modify this to save the parameters that you want)
data = zeros(nstep,8);
data(1,1) = tstart;
data(1,2) = x;
data(1,3) = func(x);  % truth value of y
data(1,4) = y;        % numerical approximation of y
data(1,5) = speed;
data(1,6) = Ep;
data(1,7) = Ek;
data(1,8) = Etot;
% calculation time loop
cnt = 1;
while cnt <= nstep
    time = cnt*dt + tstart;
    dy = (func(x+dx/2) - func(x-dx/2)) / dx;  % first derivative
    deltax = dx;     % step change in X value
    deltay = dy*dx;  % corresponding change in Y value
    mag = sqrt(deltax^2+deltay^2); % magnitude of step change
    
% compute the unit tangent vector
    Tx = deltax/mag;
    Ty = deltay/mag
    T = [Tx, Ty];  % unit tangent vector
    
% compute accelerations
    At = dot(G, T);  % acceration in the tangential direction
% update states (numerical integration)
    speed = speed + At*dt;
    delta = speed*dt;   % dstance traveled along curve
    x = x + delta*Tx;   % updated X position
    y = y + delta*Ty;   % updated Y position
% update energy states
    Ep = gravity*y;
    Ek = 0.5*speed^2;
    Etot = Ep + Ek;
    
    cnt=cnt+1
% save data
    data(cnt,1) = time;
    data(cnt,2) = x;         % X position
    data(cnt,3) = func(x);   % truth Y position
    data(cnt,4) = y;         % calculated Y position
    data(cnt,5) = speed;
    data(cnt,6) = Ep;     % potential energy
    data(cnt,7) = Ek;     % kinetic energy
    data(cnt,8) = Etot;   % total energy
end  % end of time loop
% extract data vectors from data table
time  = data(:,1);
X     = data(:,2);
fX    = data(:,3);
Y     = data(:,4);
speed = data(:,5);
Ep    = data(:,6);
Ek    = data(:,7);
Etot  = data(:,8);
% Plot data
figure;
plot(X,fX,'r');  % truth Y vs. X
hold on;
plot(X,Y,'ob');   % calculated Y vs X
grid on;
legend('fX','Y');
%.... (add additional plots as desired)

Note that for the function Y = X^2, with X starting near zero there must be a nonzero initial velocity or it won't do anything.

Also note the use of unit vectors and dot product (no trig functions).

Speed is a signed quantity. Positive speed causes x to increase. As the point rises along the curve, the speed drops to zero, and then goes negative as it falls backward.

Numerical integration errors will result in errors in the conservation of energy. Etot will not be constant, but will drift due to integration errors. Try different values for dt and you will see different ammounts of error in the energy conservation. (smaller dt should produce smaller errros). The plot of Y vs X will also become closer to the plot of fX vs X as the time step becomes smaller.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jack McCarthy 2022-5-12

If you wouldn't mind, could you provide a brief explanation as to how this would work in 3 dimensions? I am unsure as to how I could implement the vector method above. Thanks!

请先登录，再进行评论。

How would you plot a graph which a ball then rolls down (say a y=x^2 graph)

11 个评论
显示 9更早的评论隐藏 9更早的评论

回答（2 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

How would you plot a graph which a ball then rolls down (say a y=x^2 graph)

11 个评论 显示 9更早的评论隐藏 9更早的评论

回答（2 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

11 个评论
显示 9更早的评论隐藏 9更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论