function call and getting array error

Dear Sir or Madam,

I recently learn matlab for my education purpose. I would like to know how to solve the Index error in my code as I call function in for loop and getting error. I have already check array boundary and can't find the answer for the case. I deeply request you how to solve this case.

clear
% nx and ny must be the same size with t, n_part of the DEM loop
nx=4;ny=9;
% f_ equation for collisions
f=zeros(nx,ny,9);
feq=zeros(nx,ny,9);feq_f=zeros(nx,ny,9);

% velocities for x and y direction
u=zeros(nx,ny);v=zeros(nx,ny);


% velocities on solid particle
% for x direction
U_s_x=zeros(nx,ny,9);U_px=[0,0];Omega_px=[0,0];Xpx=[0,0];
% for y direction
U_s_y=zeros(nx,ny,9);U_py=[0,0];Omega_py=[0,0];Xpy=[0,0];
% position of particle and omega_s
omega_s=zeros(nx,ny,9);

% force term Ff,Tf on solid particle
Ffx=zeros(nx,ny,9);Ffy=zeros(nx,ny,9);
Tfx=zeros(nx,ny,9);Tfy=zeros(nx,ny,9);

rho=ones(nx,ny);x=zeros(nx);y=zeros(ny);
w(9)=zeros;

% initialization of the system
w=[1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9];
cx=[1 0 -1 0 1 -1 -1 1 0];
cy=[0 1 0 -1 1 1 -1 -1 0];
c2=1./3.;

% derivative rate 
dx=1.0;dy=1.0;
x1=(nx-1)/(ny-1);y1=1.0;
dx=x1/(nx-1);
dy=y1/(ny-1);

x=(0:dx:x1);
y=(0:dy:y1);
u0=0.2;

% relaxiation time
alpha=0.02;
Re=u0*(ny-1)/alpha;
% omega = 1/tau
omega=1./(3.*alpha+0.5);


% toleriant and error
%count=0;tol=1.0e-4;error=10.;erso=0.0;

% setting velocity
for j=2:ny-1
    u(1,j)=u0;
end
t_lbm=50;
for t=1:t_lbm
    [f,Ffx,Ffy,Tfx,Tfy]=up_collition(nx,ny,u,v,cx,cy,U_s_x,U_px,Omega_px,Xpx,U_s_y,U_py,Omega_py,Xpy,feq,rho,w,feq_f,f,omega_s,Ffx,Ffy,Tfx,Tfy);
end

Index in position 1 exceeds array bounds. Index must not exceed 1.

Error in solution>up_collition (line 76)
U_s_x(i,j,k)=U_px(i,j)+Omega_px(i,j)*(x(i)+0.5*(cx(k)))-Xpx(i,j);

function[f,Ffx,Ffy,Tfx,Tfy]=up_collition(nx,ny,u,v,cx,cy,U_s_x,U_px,Omega_px,Xpx,U_s_y,U_py,Omega_py,Xpy,feq,rho,w,feq_f,f,omega_s,Ffx,Ffy,Tfx,Tfy)
% the weighting function of solid particles k Bk,
% the local solid ratio e_k divided 
e_k=[1/8 1/8 1/8 1/8 1/8 1/8 1/8 1/8 0];
e_total=1;

for j=1:ny
    for i=1:nx
        t1=u(i,j)*u(i,j)+v(i,j)*v(i,j);
        for k=1:9
            x=zeros(length(i));% later add with dx term, let dt =1;
            t2=u(i,j)*cx(k)+v(i,j)*cy(k);
            Bk(k)=e_k(k).*(0.56-0.5)/(1-e_total+(0.56-0.5));
            % velocity on solid particle
            % the equaiton is Us=UP+omega_px *(x +0.5 *c(k)*dt)-Xp
            U_s_x(i,j,k)=U_px(i,j)+Omega_px(i,j)*(x(i)+0.5*(cx(k)))-Xpx(i,j);
            U_s_y(i,j,k)=U_py(i,j)+Omega_py(i,j)*(x(i)+0.5*(cy(k)))-Xpy(i,j);
            t22=U_s_x(i,j,k)*cx(k)+U_s_y(i,j,k)*cy(k);
            t11=U_s_x(i,j,k)*U_s_x(i,j,k)+U_s_y(i,j,k)*U_s_y(i,j,k);
            feq(i,j,k)=rho(i,j)*w(k)*(1.0+3.0*t2+4.5*t2*t2-1.5*t1);
            feq_f(i,j,k)=rho(i,j)*w(k)*(1.0+3.0*t22+4.5*t22*t22-1.5*t11);
            
            f(i,j,k)=f(i,j,k)-(0.9091*(f(i,j,k)-feq(i,j,k)))+(Bk(k)*(feq_f(i,j,k)-feq(i,j,k)));
            % additional collision terms
            omega_s(i,j,k)=feq_f(i,j,k)-f(i,j,k)+((f(i,j,k)-feq(i,j,k)).*0.42);
            
            % force term and torque term for DEM
            Ffx(i,j,k)=sum(Bk(k)).*sum(omega_s(i,j,k))*cx(k);
            Ffy(i,j,k)=sum(Bk(k)).*sum(omega_s(i,j,k))*cy(k);
            Tfx(i,j,k)=sum((x(i)-Xpx)*Bk(k)).*sum(omega_s(i,j,k)).*cx(k);
            Tfy(i,j,k)=sum((x(i)-Xpy)*Bk(k)).*sum(omega_s(i,j,k)).*cy(k);
        end
    end
end
end

[Rik]

There are several indexing operations on that line. Have you tried using the debugger to stop execution so you can check the sizes of the variables?

[Thin Rupar Win]

Yes, I tried the debugger. The result shows at the same line and can't find the problem. I check only for loop in command line with n= input number For i= 1:n For j= 1:n For k= 1:n ..... ...... ..... End End End That shows no error. I think something wrong with calling function and can't find the reason.

[Rik]

You're indexing several variables. Did you check each one? One of them is resulting in an error. You best chance to solve the problem is to figure out which variable is causing the problem.

[Thin Rupar Win]

Sir, I tried checking the size and array in my file. I still cant find the reason for error.

Answer 1

[Sudharsana Iyengar]

U_s_x is a 4x9x9 matrix but U_px is 1x2 matrix. So there is no U_px(3,4) did you check on that.

            U_s_x(i,j,k)=U_px(i,j)+Omega_px(i,j)*(x(i)+0.5*(cx(k)))-Xpx(i,j);
            U_s_y(i,j,k)=U_py(i,j)+Omega_py(i,j)*(x(i)+0.5*(cy(k)))-Xpy(i,j);

[Thin Rupar Win]

Sir, I want to update that U-px with another matrix. So, I give initial condition to zero. Can you please suggest me how I should give that matrix with initial zero condition with same size of nx and ny?

[Thin Rupar Win]

Sir, thank you very much for your point. Now I can solve them with giving to initial zeros(x,y). Have a nice day.