ODE45, 3 degree of freedom functions, requires innoperable step size, or unrealistic starting conditions.

I have tried:

• Other ode solvers: the mostly give the same results except ode113, which never stops being busy.
• Removal of "options" changes nothing.
• simplifying more, not simplifying at all, no change.
• Unrealisticly fast accelerations works better.
• Replacing function with a tester that is much simpler, i.e. only 3 sin terms. This also fails to do anything much.

I have noticed

• Everything is moving faster than I would have expected, especially as acceperating from near rest.

Information I feel I am missing:

• Can ode45 actually deal with 3 degrees of freedom in a system?

I am trying to model them ovement of a trebuchet using the system described in my profile pic, using Lagrange. This has 3 degrees of freedom, and this is giving me problems with the ode45. The Lagrange expression is derived by Matlabs algebra and the rearanging depicted in the immage:

 is done:

Lagrange(1)=dL_ddtheta_dt-dL_dtheta;
Lagrange(2)=dL_ddphi_dt-dL_dphi;
Lagrange(3)=dL_ddeta_dt-dL_deta;

dotsol(1)=solve(Lagrange(1)==0,d2theta);
dotsol(2)=solve(Lagrange(2)==0,d2phi);
dotsol(3)=solve(Lagrange(3)==0,d2eta);

d2theta_wo_d2phi_eq0=subs(Lagrange(1),d2phi,dotsol(2));
d2theta_wo_d2eta_eq0=subs(Lagrange(1),d2eta,dotsol(3));

d2theta_wo_d2phi=solve(d2theta_wo_d2phi_eq0==0,d2theta);
d2theta_wo_d2eta=solve(d2theta_wo_d2eta_eq0==0,d2theta);

d2phi_wo_d2theta_eq0=subs(Lagrange(2),d2theta,d2theta_wo_d2eta);
d2eta_wo_d2theta_eq0=subs(Lagrange(3),d2theta,d2theta_wo_d2phi);

d2phi_wo_d2theta=solve(d2phi_wo_d2theta_eq0==0,d2phi);
d2eta_wo_d2theta=solve(d2eta_wo_d2theta_eq0==0,d2eta);

d2theta_wo_d2phi_wo_d2eta_eq0=subs(d2theta_wo_d2phi_eq0,d2eta,d2eta_wo_d2theta);
d2theta_wo_d2phi_wo_d2eta=solve(d2theta_wo_d2phi_wo_d2eta_eq0==0,d2theta);

Then these three expressions are put into this function as zdot(2), zdot(4) and zdot(6)

function [ zdot ] = sling_uncon_func( t,z )
%Using Lagrange expression for movemnt of a simple see-saw pendulum with fixed counterweight, for use with ode45
%and derived by matlab.
% replace zdot(2) with output from see_saw_fixedCW_LExpression   
   zdot=zeros(6,1);
  global EL el ELL ell pm CW
zdot(1)=z(2);
zdot(2)=(z(2)*((pm*(2*EL*sin(z(1))*(EL*z(2)*cos(z(1)) - ELL*z(6)*cos(z(5))) - 2*EL*cos(z(1))*(EL*z(2)*sin(z(1)) - ELL*z(6)*sin(z(5)))))/2 + (CW*(2*el*sin(z(1))*(el*z(2)*cos(z(1)) + ell*z(4)*cos(z(3))) - 2*el*cos(z(1))*(el*z(2)*sin(z(1)) + ell*z(4)*sin(z(3)))))/2) - (981*CW*el*cos(z(1)))/100 + (981*EL*pm*cos(z(1)))/100 + (CW*el*cos(z(1) - z(3))*(981*cos(z(3)) - 100*el*sin(z(1) - z(3))*z(2)^2))/100 + EL*ELL*pm*sin(z(1) - z(5))*z(6)^2 - CW*el*ell*sin(z(1) - z(3))*z(4)^2 + EL*ELL*pm*sin(z(1) - z(5))*z(2)*z(6) - CW*el*ell*sin(z(1) - z(3))*z(2)*z(4) - (EL*ELL*pm*sin(z(1) - z(5))*z(2)*((pm*z(6)*(2*ELL*sin(z(5))*(EL*z(2)*cos(z(1)) - ELL*z(6)*cos(z(5))) - 2*ELL*cos(z(5))*(EL*z(2)*sin(z(1)) - ELL*z(6)*sin(z(5)))))/2 + (981*ELL*pm*cos(z(5)))/100 + EL*ELL*pm*sin(z(1) - z(5))*z(2)^2 - (EL*ELL*pm*cos(z(1) - z(5))*((981*CW*el*cos(2*z(3) - z(1)))/2 - (981*CW*el*cos(z(1)))/2 + 981*EL*pm*cos(z(1)) + 50*CW*el^2*z(2)^2*sin(2*z(3) - 2*z(1)) + 100*EL*ELL*pm*sin(z(1) - z(5))*z(6)^2 - 100*CW*el*ell*sin(z(1) - z(3))*z(4)^2))/(50*CW*el^2 + 100*EL^2*pm - 50*CW*el^2*cos(2*z(3) - 2*z(1))) + EL*ELL*pm*sin(z(1) - z(5))*z(2)*z(6)))/(EL*ELL*pm*sin(z(1) - z(5))*z(2) + (50*EL^2*ELL^2*pm^2*z(2)*sin(2*z(5) - 2*z(1)))/(50*CW*el^2 + 100*EL^2*pm - 50*CW*el^2*cos(2*z(3) - 2*z(1)))))/(CW*el^2 + EL^2*pm - CW*el^2*cos(z(1) - z(3))^2);

zdot(3)=z(4);
zdot(4)=(981*CW*el^2*cos(z(1) - z(3))*cos(z(1)) - 981*EL^2*pm*cos(z(3)) - 981*CW*el^2*cos(z(3)) + 100*CW*el^3*sin(z(1) - z(3))*z(2)^2 + 981*EL*ELL*pm*cos(z(1) - z(5))*cos(z(3)) - 981*EL*el*pm*cos(z(1) - z(3))*cos(z(1)) + 981*ELL*el*pm*cos(z(1) - z(3))*cos(z(5)) + 100*EL^2*el*pm*sin(z(1) - z(3))*z(2)^2 + 100*CW*el^2*ell*cos(z(1) - z(3))*sin(z(1) - z(3))*z(4)^2 + 100*EL*ELL*el*pm*cos(z(1) - z(3))*sin(z(1) - z(5))*z(2)^2 - 100*EL*ELL*el*pm*cos(z(1) - z(5))*sin(z(1) - z(3))*z(2)^2 - 100*EL*ELL*el*pm*cos(z(1) - z(3))*sin(z(1) - z(5))*z(6)^2)/(100*CW*el^2*ell + 100*EL^2*ell*pm - 100*CW*el^2*ell*cos(z(1) - z(3))^2 - 100*EL*ELL*ell*pm*cos(z(1) - z(5)));

zdot(5)=z(6);
zdot(6)=((pm*z(6)*(2*ELL*sin(z(5))*(EL*z(2)*cos(z(1)) - ELL*z(6)*cos(z(5))) - 2*ELL*cos(z(5))*(EL*z(2)*sin(z(1)) - ELL*z(6)*sin(z(5)))))/2 + (981*ELL*pm*cos(z(5)))/100 + EL*ELL*pm*sin(z(1) - z(5))*z(2)^2 - (EL*ELL*pm*cos(z(1) - z(5))*((981*CW*el*cos(2*z(3) - z(1)))/2 - (981*CW*el*cos(z(1)))/2 + 981*EL*pm*cos(z(1)) + 50*CW*el^2*z(2)^2*sin(2*z(3) - 2*z(1)) + 100*EL*ELL*pm*sin(z(1) - z(5))*z(6)^2 - 100*CW*el*ell*sin(z(1) - z(3))*z(4)^2))/(50*CW*el^2 + 100*EL^2*pm - 50*CW*el^2*cos(2*z(3) - 2*z(1))) + EL*ELL*pm*sin(z(1) - z(5))*z(2)*z(6))/(EL*ELL*pm*sin(z(1) - z(5))*z(2) - (100*EL^2*ELL^2*pm^2*cos(z(1) - z(5))*sin(z(1) - z(5))*z(2))/(50*CW*el^2 + 100*EL^2*pm - 50*CW*el^2*cos(2*z(3) - 2*z(1))));

end

(after using simplify())

When I try to solve this using ode45:

option=odeset('Reltol',1e-12,'MaxStep',0.001);

[T,Z]=ode45(@sling_uncon_func,[0,20],[5*pi/6,0,-pi/2,0,0,0],option);

This doesn't return an error, but I do get a matrix full of NaN for Z. If I change inputs tohave no 0 values I get between 0.2 and 1 seconds of the motion worked out. Such as:

[T,Z]=ode45(@sling_uncon_func,[0.1,20],[5*pi/6,-1e-8,-pi/2,1e-8,0,-1e-8],option);

error/warning:

Warning: Failure at t=1.000000e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t. > In ode45 at 309 In auto_maybe at 20 where line 20 is where I call the ode45.

I have tried:

• Other ode solvers: the mostly give the same results except ode113, which never stops being busy.
• simplifying more, not simplifying at all, no change.
• Unrealisticly fast accelerations works better.
• Replacing function with a tester that is much simpler, i.e. only 3 sin terms. This also fails to do anything much.

I have noticed

• Everything is moving faster than I would have expected, especially as acceperating from near rest.

Information I feel I am missing:

• Can ode45 actually deal with 3 degrees of freedom in a system?

I can supply more code, ubt I think this question is already clutterd enough, but I don't know how to explain it else. Thanks.