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

29 views (last 30 days)
Sena Clarke on 22 Aug 2019
Commented: Jack McCarthy on 12 May 2022
How would you plot a graph which a ball then rolls down (say a y=x^2 graph)
Sena Clarke on 27 Aug 2019
Edited: Sena Clarke on 27 Aug 2019
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.

David K. on 27 Aug 2019
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!
Jim Riggs on 29 Aug 2019
Edited: Jim Riggs on 29 Aug 2019
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.

Jim Riggs on 29 Aug 2019
Edited: Jim Riggs on 29 Aug 2019
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');
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.
Jack McCarthy on 12 May 2022
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!

### Categories

Find more on Automotive Applications in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!