How to solve and plot this set of differential equations?

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
      
l=1.5;
a=180;

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:

Applied Sciences | Free Full-Text | Development of an Insect-like Flapping-Wing Micro Air Vehicle with Parallel Control Mechanism (mdpi.com) Appl. Sci. 2022, 12(7), 3509; https://doi.org/10.3390/app12073509

Answer 1

[Alan Stevens]

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.

[Ijaz Ahmed]

Well, for what value of rho would we be getting a graph similar to that one in the research article? I mean the rho is initially not zero, and setting it to a non-zero value is stil giving a graph with an L-shape.

[Alan Stevens]

"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).

Answer 2

[Sam Chak]

Hi @Ijaz Ahmed

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

[Ijaz Ahmed]

but still nothing like that of the pulley profile given by the red plot in the first figure?