清空筛选条件
清空筛选条件

Solving a differential equation using ode45

1 次查看（过去 30 天）
Melhem
Melhem 2023-2-7
评论： Melhem 2023-2-8
is it possible to solve this equation using ode45?
θ'' - µ*θ'^2 + (g/r)*(µ*cos(θ) - sin(θ)) = 0
µ, g, and r are given
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2023-2-7
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
Melhem 2023-2-7
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
Sam Chak 2023-2-7
编辑：Sam Chak 2023-2-8
Ran in:
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
  5 个评论
显示 3更早的评论
Sam Chak
Sam Chak 2023-2-8
Hi @Melhem, I have edited the code in my Answer to capture the event by the block hits the ground.

请先登录，再进行评论。

更多回答（1 个）

Jan
Jan 2023-2-7
移动：Jan 2023-2-7
θ'' - µ*θ'^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)));
% ^ ^

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by