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

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:

or

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.

Answer 1

[William Rose]

@Alejandro Nava,

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(1)=y(2);
dy(2)=m^2*(w^2*y(1)^2-v^2)/y(1);
dy(3)=y(4);
dy(4)=...;  %fill in as appropriate
end

You will pass the name carousel to ode45(). Try it. Good luck.

[Alejandro Nava]

Thank you for your reply. It really gave me a great push. I've been trying to get the code working and unfortunately it has given some errors and I am unsure of what is wrong. Can you take a look at my code?

Function:

function dy=carousel(t,y)

m=10; w=15; %define constants

v=5; %define v as a constant, or in terms of components of y.

dy(1)=y(2);

dy(2)=m^2*(w^2*y(1)^2-v^2)/y(1) + y(1)*w^2; %X DD Eq

dy(3)=y(4);

dy(4)=(-2*(y(2)*w))/y(1); %Theta DD Eq

dy = [dy(1),dy(2),dy(3),dy(4)];

end

Main file:

% Hockey Puck on a Carousel Simulation %

% Clear previous run %

clear, close all

% Set Conditions %

t=0:0.1:30; %time array

y=15; % initial condition y0

[t,y]= ode45(@carousel,t,y);

Error:

>> carousel_sim

Index exceeds the number of array elements. Index must not exceed 1.

Error in carousel (line 4)

dy(1)=y(2);

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in carousel_sim (line 8)

[t,y]= ode45(@carousel,t,y);

Thanks in advance and thanks again for the help earlier.

[William Rose]

@Alejandro Nava,

See attached code. This code runs.

The errors you encountered are due to:

1. Use a different variable name for the time span than for the time returned by ode45(). Also, it is not necessary to specify the intermediate time values when specifying the time span. Just specify the start and end times. ode45() will adjust the step size in order to mee the default accuracy requirements. In my experience, the default level of accuracy is usually good enough.
2. Your initial condition vector was not satisfactory, because it did not have the correct number of elements. You can adjust the values which I used.

I added code to plot the results. The results do not look very realistic. But I believe they are an accurate result of the equations of motion which we have defined.

[Alejandro Nava]

Thank you for your code. I managed to get a running code earlier which shows what seems to be a more accurate plot.

Let me know what you think.

[William Rose]

@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.