Solving a differential equation using ode45
显示 更早的评论
is it possible to solve this equation using ode45?
θ'' - µ*θ'^2 + (g/r)*(µ*cos(θ) - sin(θ)) = 0
µ, g, and r are given
3 个评论
Torsten
2023-2-7
After you have written the equation as a system of two first-order ODEs (the first for theta, the second for theta'), you can use ode45 to solve.
John D'Errico
2023-2-7
Why not? Did you try it? Do you have initial values? You will need two of them, of course, typically theta(0) and theta'(0), or at some point. We can't really offer too much help, as you have not provided any specifics. What are the values of those parameters? What are the initial values?
Read the help docs for ODE45, where it is explicitly described how to convert the problem into a pair of first order differential equations.
Melhem
2023-2-7
采纳的回答
Hi @Melhem
Edit: The code is revised to capture the event of the falling block. As mentioned in the problem, the symbol μ is related to the friction, which should dampen the falling motion at the beginning. The simulation stops when the block hits ground, that is when 
.
.Please check the derivation of the equations of motion again. Not sure if the signs are correct or not.
V0 = 10;
mu = 0.6;
tf = 10;
W = 1;
g = 32.2;
r = 5;
tspan = linspace(0, tf, 1001);
y0 = [0; V0/r];
param = g/r;
options = odeset('Events', @BlockHitsGroundEventFcn);
[t, y, te, ye, ie] = ode45(@(t,y) EOM(t, y, param, mu), tspan, y0, options);
figure(1)
yyaxis left
plot(t, y(:,1)*180/pi), ylabel({'$\theta$, deg'}, 'Interpreter', 'latex')
yyaxis right
plot(t, y(:,2)), ylabel({'$\dot{\theta}$, rad/s'}, 'Interpreter', 'latex')
legend('y_1', 'y_2', 'location', 'best')
xlabel('t, sec'), grid on
figure(2)
N = (- W*r*y(:,2).^2 + W*g*sin(y(:,1)))/mu;
plot(t, N), grid on
xlabel('t'), ylabel('N')
function g = EOM(t, y, param, mu)
g(1,1) = y(2);
g(2,1) = - mu*y(2).^2 - param*(mu*cos(y(1)) - sin(y(1)));
end
function [position, isterminal, direction] = BlockHitsGroundEventFcn(t, y)
position = y(1) - pi/2; % When theta = 90 deg
isterminal = 1; % Halt integration
direction = 1; % When theta is increasing from 0 to 90 deg
end
5 个评论
Melhem
2023-2-8
Hello, Thanks for your reply.
I wrote the equations wrong and they should be:
g(2,1) = mu*y(2).^2 - param*(mu*cos(y(1)) - sin(y(1)))
N = - W*r*y(:,2).^2 + W*g*cos(y(:,1))
Now I'm getting this error for mu = 0.6:
Warning: Failure at t=8.949691e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at
time t.
if I take tf=0.7, the error won't appear
I'm taking mu=0:0.1:0.6
Sam Chak
2023-2-8
Hi @Melhem
The solution of your motion system (original ODE) blows up in finite time due to the nature of the differential equation.
Positive μ means that the kinetic energy (from the angular velocity 
) is being added to the motion system. Also this unbounded term 
dominates the bounded term 
, causing the acceleration to keep increasing, and eventually θ and 
blow up in finite time.
) is being added to the motion system. Also this unbounded term dominates the bounded term , causing the acceleration to keep increasing, and eventually θ and blow up in finite time.However, I suspect that this quadratic term 
is related to the drag force, which should dampen the motion system over time. That's why I simulated the way it was in my Answer. Please check the Equation of Motion (EoM) again.
is related to the drag force, which should dampen the motion system over time. That's why I simulated the way it was in my Answer. Please check the Equation of Motion (EoM) again.Similar to mass-damper-spring system (stable)
,when the ODE is rearranged in this form
, you can see that the sign is 'minus'.
Melhem
2023-2-8
I'm trying to find the angle θd at which the block leaves the cylindrical surface, So I found N=0 and used it to plot θd with respect to time
This is the lab assignment:
Thank you anyways, you helped a lot!
Melhem
2023-2-8
更多回答（1 个）
θ'' - µ*θ'^2 + (g/r)*(µ*cos(θ) - sin(θ)) = 0 means:
θ'' = µ*θ'^2 - (g/r)*(µ*cos(θ) - sin(θ))
This does not match:
g(2,1) = -mu*y(2).^2 + param*(mu*cos(y(1)) - sin(y(1)));
% ^ ^
类别
在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Translated by 
