Solving a differential equation using ode45

6 visualizzazioni (ultimi 30 giorni)
Melhem il 7 Feb 2023
Commentato: Melhem il 8 Feb 2023
is it possible to solve this equation using ode45?
θ'' - µ*θ'^2 + (g/r)*(µ*cos(θ) - sin(θ)) = 0
µ, g, and r are given
3 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
John D'Errico il 7 Feb 2023
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 il 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)));

Accedi per commentare.

Risposta accettata

Sam Chak il 7 Feb 2023
Modificato: Sam Chak il 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
5 CommentiMostra 3 commenti meno recentiNascondi 3 commenti meno recenti
Sam Chak il 8 Feb 2023
Hi @Melhem, I have edited the code in my Answer to capture the event by the block hits the ground.
Melhem il 8 Feb 2023
Thank you so much you're a life saviour.

Accedi per commentare.

Più risposte (1)

Jan il 7 Feb 2023
Spostato: Jan il 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)));
% ^ ^
0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by