# Solving a differential equation using ode45

3 views (last 30 days)
Melhem on 7 Feb 2023
Commented: Melhem on 8 Feb 2023
is it possible to solve this equation using ode45?
θ'' - µ*θ'^2 + (g/r)*(µ*cos(θ) - sin(θ)) = 0
µ, g, and r are given
Melhem on 7 Feb 2023
yes I've tried it but the values don't make sense
here's the code:
In here I'm also trying to find the value of the normal force N
v0=10
mu=0.2
tf=1
function [t,y,N] = lab1(V0,mu,tf)
W=1;
g=32.2;
r=5;
tspan=linspace(0,tf,1000);
y0=[0;V0/r];
param=g/r;
[t,y]=ode45(@(t,y)EOM(t,y,param,mu),tspan,y0);
N=(-W*r*(y(:,2)).^2 + W*g*sin(y(:,1)))/mu;
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)));

Sam Chak on 7 Feb 2023
Edited: Sam Chak on 8 Feb 2023
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
Melhem on 8 Feb 2023
Thank you so much you're a life saviour.

Jan on 7 Feb 2023
Moved: Jan on 7 Feb 2023
θ'' - µ*θ'^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)));
% ^ ^