Azzera filtri
Azzera filtri

Solving a differential equation using ode45

6 visualizzazioni (ultimi 30 giorni)
is it possible to solve this equation using ode45?
θ'' - µ*θ'^2 + (g/r)*(µ*cos(θ) - sin(θ)) = 0
µ, g, and r are given
  3 Commenti
John D'Errico
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
function [t,y,N] = lab1(V0,mu,tf)
N=(-W*r*(y(:,2)).^2 + W*g*sin(y(:,1)))/mu;
function g=EOM(t,y,param,mu)
g(2,1)=-mu*y(2).^2 + param*(mu*cos(y(1)) - sin(y(1)));

Accedi per commentare.

Risposta accettata

Sam Chak
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);
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
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)));
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
  5 Commenti
Sam Chak
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)));
% ^ ^


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