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)));
% ^ ^

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by