'Matrix dimensions must agree.' while invoking ode15i

I would like to solve a second order ode system consisting of two second order implicit differential equations (this system describes the motion of a point of mass and the centroid of a rigid body (a brick), if the mass is attached to the rigid body via a rope and swings as a pendulum; the motion is taken place in the xy plane and is given by two coordinates x and phi).

Here are the equations of motion:

(m1+m2)*x'' + m2*l*phi''*cos(phi) - m2*l*(phi')^2*sin(phi) = 0 x''*l*cos(phi) + l^2*phi'' = -g*l*sin(phi)

let X = [x x' phi phi'] and Y = X' = [x' x'' phi' phi''], then here is my function:

function Z = implicit(t, X, Y) m1 = 1; % mass of the brick m2 = 1; % mass of the l = 5; % length of the string g = 9.81; % gravitational constant Z = zeros(4,1); Z(1) = (m1+m2)*Y(2)+m2*l*Y(4).*cos(X(3))-m2*l*(Y(3)).^2.*sin(X(3)); Z(2) = Y(2).*cos(X(3))+l*Y(4)+g*sin(X(3)); Z(3) = Y(1)-X(2); Z(4) = Y(3)-X(4);

My first question: have I implemented the system well?

After, I called the solver this way: [t y] = ode15i(@implicit, [0 2], [0 0 pi/4 0], [0 0 0 0]);

Now the initial conditions are: initial coordinate, velocity, acceleration of the brick is 0, initial angle of the pendulum is pi/4, initial angular velocity and angular acceleration is 0.

And I get the following error:

??? Error using ==> max Matrix dimensions must agree.

Error in ==> odenumjac at 103 yscale = max(abs(y),thresh);

Error in ==> ode15ipdupdate at 30 [dfdyp,dfdyp_options.fac,NF] = odenumjac(odefun,{t,y,yp,extras{:}},f,dfdyp_options);

Error in ==> ode15ipdinit at 76 [dfdy,dfdyp,dfdy_options,dfdyp_options,nfcn] = ...

Error in ==> ode15i at 175 [Jac,dfdy,dfdyp,Jconstant,dfdy_options,dfdyp_options,nfcn] = ...

Please help me.

Thank you, András

Answer 1

[Zoltán Csáti]

I managed to solve those equations explicitly.

[Jan]

The explicit solution is not necessarily more efficient for numerical simulations or visualizations. Solving a small linear linear system in each step can be cheaper than a bunch of trigonometric functions.

However, currently the original question is hard to read due to the missing code formatting. Please edit the question - code formatting is trivial. Perhaps you want to post the solution?