ode45 --- Index exceeds matrix dimensions

THIS IS MY FUNCTION

function dzdt=sisedo(t,x,N,theta,mu,gamma,kappa,ka,ca)

      dzdt=[]

      for i=1:N
          j=(i-1)*6+1
          %Equations
          dzdt=[dzdt
                x(j+1)
                (-kappa*(2*x(j)-x(j+6))-(x(j)-x(j+2)-theta*x(j+4)))/mu  %-x(j-6)
                x(j+3)
                -(x(j+2)-x(j)-theta*x(j+4))-Faero(x(j+2),x(j+3),ca,ka)
                x(j+5)
                (-(2*theta*x(j+4)-x(j+2)+x(j))-Fric(x(j+4),x(j+5),i))/(theta*gamma)]; 
      end
  end

THIS IS THE MAIN SCRIPT

%%%%%%Parameters
N=10;
theta=0.001;
mu=1;
gamma=1;
kappa=1;
ka=0.001;
ca=-0.0001;

tspan=[0 3];

x0=zeros(6*N,1);

for j=1:6:6*N

      x0(j)=0;

      x0(j+1)=1;

      x0(j+2)=0;

      x0(j+3)=1;

      x0(j+4)=0;

      x0(j+5)=0;

end

opciones=odeset('abstol',1e-9,'reltol',1e-9)

[t,x]=ode45(@sisedo,tspan,x0,opciones,N,theta,mu,gamma,kappa,ka,ca);

[MANUEL PRECIADO VIDAL ARAGON]

Index exceeds matrix dimensions.

Error in sisedo (line 22)

 (-kappa*(2*x(j)-x(j+6))-(x(j)-x(j+2)-theta*x(j+4)))/mu  %-x(j-6)

Error in odearguments (line 87)

 f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

 odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in SIMULADOR_N (line 43) %THIS IS THE MAIN SCRIPT

 [t,x]=ode45(@sisedo,tspan,x0,opciones,N,theta,mu,gamma,kappa,ka,ca);

Answer 1

[Walter Roberson]

We cannot test that without Faero and Fric

Your error message is against a line number which exceeds the number of lines of code you show us.

You should be parameterizing the function rather than counting ode45 to pass extra parameters to the function. That behaviour of ode45 has not been documented for over a decade, having been functionally replaced as of MATLAB 5.1.

[t,x]=ode45(@(t,x) sisedo(t,x,N,theta,mu,gamma,kappa,ka,ca),tspan,x0,opciones);

For efficiency your sisedo should be pre-allocating the output instead of growing the output variable.

function dzdt = sisedo(~,x,N,theta,mu,gamma,kappa,ka,ca)

        dzdt = zeros(6*N,1);

        for i=1:N
            j=(i-1)*6+1;
            %Equations
            dzdt(j:j+5) = [
                  x(j+1)
                  (-kappa*(2*x(j)-x(j+6))-(x(j)-x(j+2)-theta*x(j+4)))/mu  %-x(j-6)
                  x(j+3)
                  -(x(j+2)-x(j)-theta*x(j+4))-Faero(x(j+2),x(j+3),ca,ka)
                  x(j+5)
                  (-(2*theta*x(j+4)-x(j+2)+x(j))-Fric(x(j+4),x(j+5),i))/(theta*gamma)]; 
        end
    end

As for the out of range:

Your x is defined as being 6*N long. Look at the behaviour when i reaches N in the for loop. j will become (N-1)*6+1 so with N = 10 that would be j = 55, allowing for x(j) = x(55), x(j+1) = x(56), x(j+2) = x(57), x(j+3) = x(58), x(j+4) = x(59), x(j+5) = x(60) . With x being 60 long, x(j+6) = x(61) will not exist. But you have

(-kappa*(2*x(j)-x(j+6))-(x(j)-x(j+2)-theta*x(j+4)))/mu

which refers to x(j+6) so that is always going to be out of bounds when i is maximum. You need to recheck your equations there to figure out why that one term is referring to the "next" group of x values.

[MANUEL PRECIADO VIDAL ARAGON]

I changed the code as you said, but with no result. I have the same errors.

But I understand what you mean. I will try to get a solution in a different way.

The Friction and Faero functions are the followings:

function F=Fric(y,yp,i)

      global F0 y0 yp0

      deltay=abs(y-y0(i));
      if deltay>2
          disp('deltay>2, explota')
      end
      deltaF=2*(1-(1-deltay/2)^(5/2));
      F=F0(i)+sign(yp0(i))*deltaF;

      if sign(yp)~=sign(yp0(i)) 
          F0(i)=F;
          y0(i)=y;
          yp0(i)=sign(yp);
      end
  end

function A=Faero(x,xp,ca,ka)

      A=-(ca*xp + ka*x)

   end

Thank you very much