How to solve and plot this set of differential equations?

2 views (last 30 days)
Hello everyone,
I am trying to plot this set of differential equations, but so far I am getting some weird or no graphs at all.
we have the following equations:
and we want the red plot as under which is for a pulley profile:
and for the reference, we could use the previous graph on that research article as well:
and the code which I used recently for this was:
% code
f = @(t,y) [(l*y(1)*cos(a-t)-(y(1))^2)/(l*sin(a-t));y(1)+t*((l*y(1)*cos(a-t)-(y(1))^2)/(l*sin(a-t)))];
tspan = [0, 360];
% p(0) m(0)
xinit = [0, 15];
ode45(f, tspan, xinit)
legend('p(t)', 'm(t)')
For reference, here is the link to that article:

Answers (2)

Alan Stevens
Alan Stevens on 5 Feb 2023
If rho starts at 0, then drho/dtheta will be zero and rho will not change from zero according to your first ode. If rho and drho/dtheta both stay at zero then m will not change either, according to your second ode.
Alan Stevens
Alan Stevens on 6 Feb 2023
Edited: Alan Stevens on 6 Feb 2023
"Well, for what value of rho would we be getting a graph similar to that one in the research article?."
From the graph it looks like something close to 2.
However, it seems you are also mixing radians (using sin and cos) with degrees (using 180 and 360).
Also, the graph is a polar plot, so it's probably plotting y vs x, where y is rho*sin(theta) and x is rho*cos(theta).

Sign in to comment.

Sam Chak
Sam Chak on 6 Feb 2023
I attempt to test the first of the first-order differential equations because it does not depend on m.
The relationship between m and ρ is also given by the Law of Cosines:
Note that singularities (division by zero) occur at , where .
l = 1.5;
odefcn = @(t, x) cot(pi - t)*x - 1/(l*sin(pi - t))*x.^2;
tspan = [0.1 3.122]; % theta from 0.1 rad to 3.122 rad
x0 = 1; % initial rho
[t, x] = ode45(odefcn, tspan, x0);
plot(t, x), grid on, xlabel('t')
  1 Comment
Ijaz Ahmed
Ijaz Ahmed on 25 Feb 2023
but still nothing like that of the pulley profile given by the red plot in the first figure?

Sign in to comment.


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