Solving set of coupled non-linear ode using ode45

Hello everyone,

I am tryingt to solve a set of coupled non-linear differential equation using ode45 but i am not getting the desired results. By desired results I mean , setting all the initial conditions to be zero and setting torques for both joints , there should be no change in coordinate or change in velocity of the manipulator in other words if you plot the solution of the ode . It should be a horizontal line parallel to time axis. But this is not the case when i run the code. Given below are the set of equations that i am trying to solve numerically:

And this is the code that i am using to solve the above system :

function xdot = DynOde(t,y)
%% init constants;
m1 = 5;
m2 = 2;
a1 = 0.34;
a2 = 0.34;
g = 9.81;
T1 = 0;
T2 = 0;


x1dot  = y(2);
x1ddot = (T1*a2 - 2*a2 - 2*a1*cos(y(3)) - a1*a2*g*m1*cos(y(1)) - a1*a2*g*m2*cos(y(1)) + a1*a2^2*m2*sin(y(3))*y(2)^2 + a1*a2^2*m2*sin(y(3))*y(4)^2 + a1*a2*g*m2*cos(y(3))*cos(y(1) + y(3)) + a1^2*a2*m2*cos(y(3))*sin(y(3))*y(2)^2 + 2*a1*a2^2*m2*sin(y(3))*y(2)*y(4))/(a2*(a1^2*m1 + a1^2*m2 - a1^2*m2*cos(y(3))^2));
x2dot = y(4) ;
x2ddot = (T1*a2 - 2*a2 - 2*a1*cos(y(3)) - a1*a2*g*m1*cos(y(1)) - a1*a2*g*m2*cos(y(1)) + a1*a2^2*m2*sin(y(3))*y(2)^2 + a1*a2^2*m2*sin(y(3))*y(4)^2 + a1*a2*g*m2*cos(y(3))*cos(y(1) + y(3)) + a1^2*a2*m2*cos(y(3))*sin(y(3))*y(2)^2 + 2*a1*a2^2*m2*sin(y(3))*y(2)*y(4))/(a2*(a1^2*m1 + a1^2*m2 - a1^2*m2*cos(y(3))^2));

xdot = [x1dot;x1ddot;x2dot;x2ddot];
end

Any advice would be of great help.

Than you.

[David Goodmanson]

Hi Ajai,

what are your initial conditions?

[Ajai Singh]

To check if the code is working , give [ 0 0 0 0 ] . And if everything thing is fine then the plot of solution of the ode should be horizontal lines parallel to time axis.

Answer 1

[David Goodmanson]

Hi Ajai,

The angles appear to be measured from the horizontal, not the vertical. That's because the terms involving g use cos(theta) and not sin(theta). So if you put in theta = 0 for the initial conditions, the bars are horizontal. So of course they move due to the torque provided by gravity.

[Ajai Singh]

This is what I was thinking but shouldn't the plots come to a saturation if they are being pulled by gravity? Also I wanted to know if the function that I wrote correctly reflects the set of equations given in the picture.

[David Goodmanson]

By saturation do you mean that the system will eventually come to a halt? Unless there are frictional-type loss terms, this system should oscillate indefinitely, like a pendulum does.

[Ajai Singh]

Ah!!! That explains it, but if you try with setting torques to non-zero values. For example set T1=5,T2=7. Matlab throws a warning and the plots generated do not saturate. By saturation I mean : if we apply a torque at both joints then after some time the joints should have fixed coordinates in the XY plane since it's a planar Robot.

[David Goodmanson]

Hi Ajai,

are you saying that with fixed torques the system should settle down to a static solution? I don't think that follows. You can do it with time-dependent torques, but with constant torques there is still the opportunity for undamped oscillations.

[Ajai Singh]

Thank you so much for the explanation. I get it now.Just to confirm, all that I need to do is : instead of having constant torque, make it a function of time ?