I can't figure out how to write equations of motion into Matlab.

9 views (last 30 days)
I'm attempting to code a simulation for a hockey puck on a playground carousel which I have found the equations of motion for.
The equations are:
With m being the mass, omega (and theta dot) being the angular velocity, x being the distance from the puck to the center of the carousel, and v being the tangental velocity. Cf is the centrifugal velocity and Cp is the centripetal velocity.
I've been trying to utilize the ode45 function and directly writing down the equations to no avail. Can I get a nudge into the right direction? Any help would be greatly appreciated.

Accepted Answer

William Rose
William Rose on 21 Nov 2022
Please read the examples in ode45() help pages carefully, They are very helpful.
You have two angular velocities: ω and . Is ω a constant, the carousel rotation rate? Is Θ the time-varying angular position of the puck, in an inertial, i.e. non-rotating, reference frame?
You say v = tangential velocity. Is that in the non-rotating frame? Can you express v in terms of x and Θ, or in terms of and ?
You say Cp=centripetal velocity and Cf=centrifugal velocity. I think you meant to say Cp=centripetal force, etc.
You will need to write the equations of motin as a set of first order differential equations. I suspect you will need four variables: two equations for x and two for Θ. The vector y will have elements . I assume that you can write the velocity v as a function of x and theta.
You will need to create a derivative function which returns a 1x4 vector whose elements are . For example, your function could look like this:
function dy=carousel(t,y)
m=1; w=1; %define constants
v=...; %define v as a constant, or in terms of components of y.
dy(4)=...; %fill in as appropriate
You will pass the name carousel to ode45(). Try it. Good luck.
William Rose
William Rose on 22 Nov 2022
@Alejandro Nava, the plot looks great!
Cpngratulations on getting it working.
I see that the state vector now has six elements, rather than four. The code to convert from polar to Cartesian coordinates suggests that y(:,1)=radius and y(:,2)=theta. I believe that theta is measured in the rotating frame, because the equaitons of motion in your original posting appear to be a form of the Coriolis equaitons for movement relative to a rotating frame of reference.
The code for the differential equations, in carousel.m, specifies that , and , , . I do not understand why these are the correct derivatives.
Finally, carousel.m indicates that , which means . Likewise, carousel.m indicates that , which means .
The specified equatons for y5 and y6 in carousel.m suggest that you consider y3 to equal dr/dt, and that you consider y4 to equal dTheta/dt. However, as noted above, the equaitons as written also suggest y5 equals dr/dt, and y6 equuals dTheta/dt. Thus there appears to be an inconsistency.
In the equaiton for y5, you use y4/R in place of v. Is this correct? I suppose it depends in part on the definition of v. Is v the puck velocity? If so, measured in what frame? Why do you use a constant value of R in the equaiton for v? The equation gives the velocity v in the NON-rotating frame of an object moving on a circular path, at a constant distance r from the center. The puck does not move on a circular path with constant r, so does not describe the puck's velocity in the inertial frame.
The constant w, defined in carousel.m, does not appear to be used. Why not? I am afraid that its non-use reflects confusion about the equivalence and/or role of ω and .
carousel.m defines R=3, in line 2. This means the value of R that is passed from the calling program is not used. Perhaps that is a mistake.
Does the frame rotate at a constant rate? In your original Coriolis equations of motion, ω is the angular velocity of the frame, relative to an inertial frame. Do not confuse ω with . They are not equal.
I am sorry if all my comments seem overly critical. The comments only show that I do not understand your code.

Sign in to comment.

More Answers (0)




Community Treasure Hunt

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

Start Hunting!

Translated by