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

[Torsten]

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]

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]

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

Answer 1

[Sam Chak]

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

[Melhem]

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]

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.

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.

Similar to mass-damper-spring system (stable)

,

when the ODE is rearranged in this form

, you can see that the sign is 'minus'.

[Melhem]

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!

[Sam Chak]

Hi @Melhem, I have edited the code in my Answer to capture the event by the block hits the ground.

[Melhem]

Thank you so much you're a life saviour.

Answer 2

[Jan]

θ'' - µ*θ'^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)));
%        ^           ^