function call and getting array error

2 次查看（过去 30 天）
Thin Rupar Win
Thin Rupar Win 2021-11-9
评论： Steven Lord 2024-1-10
Ran in:
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
  4 个评论
显示 2更早的评论
Rik
Rik 2021-11-9
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
Thin Rupar Win 2021-11-9
Sir, I tried checking the size and array in my file. I still cant find the reason for error.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sudharsana Iyengar
Sudharsana Iyengar 2021-11-9
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);
  3 个评论
显示 1更早的评论
Thin Rupar Win
Thin Rupar Win 2021-11-9
Sir, thank you very much for your point. Now I can solve them with giving to initial zeros(x,y). Have a nice day.
Steven Lord
Steven Lord 2024-1-10
FYI, one way to help debug this type of issue is to set an error breakpoint, run the code to reproduce the error, and examine the variables used on the line where MATLAB enters debug mode to make sure you're not trying to retrieve an element that doesn't exist of one of those variables.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by