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

Question

Alejandro Nava ，2022-11-21

0
链接

Direct link to this question

https://ww2.mathworks.cn/matlabcentral/answers/1857933-i-can-t-figure-out-how-to-write-equations-of-motion-into-matlab

评论： William Rose ，2022-11-22

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.

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

请先登录，再进行评论。

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

Answer 1

William Rose 2022-11-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1857933-i-can-t-figure-out-how-to-write-equations-of-motion-into-matlab#answer_1106743

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

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

Alejandro Nava 2022-11-21

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 2022-11-22

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

请先登录，再进行评论。