MATLAB: Error using ==> mrdivide: Matrix dimensions must agree

matrixmrdivideodeode45

[EDIT: 20110609 14:53 CDT – reformat – WDR]
I've been working on an ODE modeling the slipping nature of sand in a cylindrical shell as it rolls down a hill. When I run my code, this error appears:
Error using ==> mrdivide
Matrix dimensions must agree.
Error in ==> sr>f2 at 135
II= (m*a*b*sin(y(1)+alpha)-Iw)/Q+a/J;
Can someone explain to me why this is happening? My initial conditions for this ODE are in a column vector (eg y0 = [PHI;0;0;0];) and I am using ode45 to solve this equation. y0 has been transformed so that it is not confused with a vector (yout = y0.';). Below is the code(note that the sticking system has no errors):
function sr
global T b My AREA a c h STEEL SAND M m alpha Ic Iw g
tstart = 0;
tfinal = 15;
T=0.05; %thickness of wedge
b=0.10; %inner radius
My=2/3*sqrt(T^3*(2*b-T)^3);
AREA= pi*b^2/2-b^2*asin(1-T/b)+(T-b)*sqrt(T*(2*b-T));
a=My/AREA; %centroid calculator (exact)
c=0.11; %outer radius
h=0.55; %height of cylinder
STEEL=7850;
SAND =1602;
M=pi*(c^2-b^2)*h*STEEL;
m=AREA*h*SAND; 
alpha=pi/180*5;
Ic=M/2*(c^2+b^2);
Iw=m*a^2;
g=9.81;
mus=0.5; mud=mus;
PHI=acos((1+M/m)*b/a*sin(alpha));
y0 = [pi/5;0;0;0]; %initial internal angle of wedge, internal angular
%velocity of wedge, position, and velocity of cylinder respectively.
tout = tstart;
yout = y0.'; % use .' at the end so it's not confused with a vector.
teout = [];
yeout = [];
ieout = [];
Q=m*a*sin(y0(1)+alpha)-Ic/b-b*(M+m);
W=Q+m*a*sin(y0(1)+alpha)-Iw/b;
A=1/W*(m*a*g*cos(y0(1)) -(m+M)*b*g*sin(alpha) +m*a*b*y0(2)^2);
PSIDOT=-A/b;
FN=(-A*sin(y0(1)+alpha) - a*PSIDOT + g*cos(y0(1)))/...
    (A*cos(y0(1)+alpha) + a*y0(2)^2 + g*sin(y0(1)));
Z= mus-FN;
for i = 1:30 
options = odeset('Events',@events2,'RelTol',1e-8,'AbsTol',1e-12);
    if tstart == tfinal
        break;
    end;
  [t,y,te,ye,ie] = ode45(@f2,[tstart tfinal],y0,options);
    nt = length(t);
    tout = [tout; t(2:nt)];
    yout = [yout; y(2:nt,:)];
    teout = [teout; te];
    yeout = [yeout; ye];
    ieout = [ieout; ie];
    y0(1) = y(end,1);
    y0(2) = y(end,2);
    y0(3) = y(end,3);
    y0(4) = y(end,4);      
    tstart = t(nt);
options = odeset('Events',@events1,'RelTol',1e-8,'AbsTol',1e-12);
    if tstart == tfinal
        break;
    end;
    [t,y,te,ye,ie] = ode45(@f1,[tstart tfinal],y0,options);
    nt = length(t);
    tout = [tout; t(2:nt)];
    yout = [yout; y(2:nt,:)];
    teout = [teout; te];     
    yeout = [yeout; ye];
    ieout = [ieout; ie];
    y0(1) = y(end,1);
    y0(2) = y(end,2);
    y0(3) = y(end,3);
    y0(4) = y(end,4);
    tstart = t(nt);
  end;
      figure(2)
      subplot(211)
      plot(tout,yout(:,1)) %,teout,yeout(:,1),'ro')
      xlabel('time');
      ylabel('\phi');
      subplot(212)
      plot(tout,yout(:,2)+yout(:,4)/b) %,teout,yeout(:,2),'ro')
      xlabel('time');
      ylabel('d\phi/dt');
      figure(1)
      subplot(211)
      plot(tout,yout(:,3)) %,teout,yeout(:,3),'ro')
      xlabel('time');
      ylabel('x');
      subplot(212)
      plot(tout,yout(:,4)) %,teout,yeout(:,4),'ro')
      xlabel('time');
      ylabel('v');
%for all systems, masses are in kg, lengths are in m, moments of
%intertia are in kg*m^2, and densities are in kg/m^3
function dydt = f1(t,y)
%sticking system
global b a M m alpha Ic Iw g
Q=m*a*sin(y(1)+alpha)-Ic/b-b*(M+m);
W=Q+m*a*sin(y(1)+alpha)-Iw/b;
A=1/W*(m*a*g*cos(y(1)) -(m+M)*b*g*sin(alpha) +m*a*b*y(2)^2*cos(y(1)+alpha));
PSIDOT=-A/b;
dydt = [-y(4)/b;PSIDOT;y(4);A];
function dydt = f2(t,y)
%slipping system
global b a M m alpha Ic Iw g mud
Q=m*a*sin(y(1)+alpha)-Ic/b-b*(M+m);
J=sin(y(1)+alpha)+mud*cos(y(1)+alpha)*sign(y(2)+y(4)/b);
II= (m*a*b*sin(y(1)+alpha)-Iw)/Q+a/J;
PSIDOT=1/II*((g*cos(y(1))-mud*sign(y(2)+y(4)/b)*(g*sin(y(1))+a*y(2)^2))/J-(m*a*g*cos(y(1))-(m+M)*b*g*sin(alpha)+m*a*b*y(2)^2*cos(y(1)+alpha))/Q);
A=1/Q*(m*a*g*cos(y(1))-(m+M)*b*g*sin(alpha)+m*a*b*y(2)^2*cos(y(1)+alpha)+PSIDOT*(m*a*b*sin(y(1)+alpha)-Iw));
dydt = [y(2);PSIDOT;y(4);A];
% --------------------------------------------------------------------------
function [value,isterminal,direction] = events1(t,y)
global b a M m alpha Ic Iw g mus
Q=m*a*sin(y(1)+alpha)-Ic/b-b*(M+m);
W=Q+m*a*sin(y(1)+alpha)-Iw/b;
A=1/W*(m*a*g*cos(y(1)) -(m+M)*b*g*sin(alpha) +m*a*b*y(2)^2);
PSIDOT=-A/b;
FN=(-A*sin(y(1)+alpha) - a*PSIDOT + g*cos(y(1)))/...
    (A*cos(y(1)+alpha) + a*y(2)^2 + g*sin(y(1)));
value = (mus-FN); %finds critical zero
isterminal = 1;   % stops the integration

direction = 0;  
function [value,isterminal,direction] = events2(t,y)
global b
value =y(2)+y(4)/b; 
isterminal = 1;   % stops the integration
direction = 0;

Best Answer

Change your line
II= (m*a*b*sin(y(1)+alpha)-Iw)/Q+a/J;
to
II= (m.*.a.*.b.*sin(y(1)+alpha)-Iw)./Q + a./J;
I did not trace far enough to see which of your variables is triggering the problem, but if that change gets rid of the division problem then you could put in a breakpoint and display the sizes of the variables.
Related Question