请教matlab程序错误

haorizi2012

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A LATTICE BOLTZMANN FORMULATION FOR THE ANALYSIS OF
% RADIATIVE HEAT TRANSFER PROBLEMS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in matlab
% Copyright (C) 2010-12-10 14:53:49 Zhao Lei .
% Address:
% E-mail: haorizi2008.cool@163.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clf

% RADIATION CONSTANS

ly     = 51;
lx     = 51;                        
delta_x= 1./(lx-1);
delta_y= 1./(ly-1);
delta_t= 1.0e-004;
sigma = 5.67e-008;          % constant s is the geometric distance
Thot  = 1;
Tcold = 0;
Ib    = sigma*Thot.^4/pi ;  % blackbody intensity
G     = 4*pi.*Ib ;          % volumetric absorption
IEq   = G/(4*pi);           % equilibrinm PDFs
beta  = 1.0;                % extinction coefficient
n   =8;
alpha = pi/4;      %
tPlot = 100;       % iterations between successive graphical outputs
tStatistics = 10;  % iterations between successive file accesses
Kn = delta_x*beta <= 0.05;  % Kn number is smaller than a threshold value that ensures stability and accuracy

% D2Q8 LATTICE CONSTANTS
cx  = [  1,   0,   -1,   0,   1,   -1,  -1,   1  ];
cy  = [  0,   1,    0,  -1,   1,   1,   -1,  -1  ];

t   = ones(1,8)*alpha/(2*pi);       % weight
az  = [0:2*pi/n:7/4*pi];      
[y,x] = meshgrid(1:ly,1:lx);


ux  = delta_x/delta_t;   % velocity
uy  = delta_y/delta_t;
omega = sqrt((cx*ux)*(cx*ux)'+(cx*uy)*(cx*uy)')*beta;    % relaxation frequency


% INITIAL CONDITION FOR RADIATION:  
IIn  = ones(8,lx,ly);   %8 is the direction number,2 is r
IOut = ones(8,lx,ly);
ts = 0;

% MAIN LOOP
while (ts<80000 && Kn <= 0.05 && abs((abs(IOut)-abs(IIn))./abs(IIn)))<1e-6 || ts<100
      
    %  BOUNDARY CONDITIONS FOR RADIATION:
    for i=1:8
        IIn (i , : ,1) = sigma*Thot^4/pi;  % sigma is Stefan-Boltzmann constant
        IOut(i , 1 ,:) = 0;
        IOut(i ,lx ,:) = 0;
        IOut(i ,: ,ly) = 0;
    end
      % COLLISION STEP

    for i=1:8
        IEq(i,:,:)  = zeros(i,lx,ly);  %让平衡分布函数IEq初始为同维的零矩阵
        IEq(i,:,:)  = IEq(i,:,:)+ reshape(t*reshape(IIn,8,lx*ly),i,lx,ly);       % equilibrinm PDFs
        IOut(i,:,:) = zeros(i,lx,ly);
        IOut(i,:,:) = IIn(i,:,:) - omega .* (IIn(i,:,:)-IEq(i,:,:));
    end
   
      % STREAMING STEP RADIATION
    for i=1:8
       IIn(i,:,:) = circshift(IOut(i,:,:), [0,cx(i),cy(i)]);
       tx  = pi*sin(alpha/2)*cos(az(i));     %weight x direction
       ty  = pi*sin(alpha/2)*sin(az(i));     %weight y direction
    end
    ts = ts+1;
end
    % radiative heat fluxes
    Qx = sum(tx.*IOut);
    Qy = sum(ty.*IOut);
  
    % VISUALIZATION
    if (mod(ts,tPlot)==0)
       Q = reshape(sqrt(Qx.^2+Qy.^2),lx,ly);
       imagesc(Q');colorbar
       title('heat flux');
       axis equal off; drawnow
    end


耽误一下大家的时间，希望能得到大家的关注

mhkmars

回复 1 # haorizi2012 的帖子

你不说你的问题，光贴个程序出来，你想让网友如何关注？

haorizi2012

回复 2 # mhkmars 的帖子

你运行一下就知道了啊，主要是  IEq(i,:,:)  = IEq(i,:,:)+ reshape(t*reshape(IIn,8,lx*ly),i,lx,ly);  这句，
??? Index exceeds matrix dimensions.  

meiyongyuandeze

haorizi2012 发表于 2011-3-30 19:03

                               
登录/注册后可看大图

回复 2 # mhkmars 的帖子

你运行一下就知道了啊，主要是  IEq(i,:,:)  = IEq(i,:,:)+ reshape(t*reshape(I ...

估计是你数组维数没有定义好，

雨人

你的错误原因是矩阵访问超出边界了。比如矩阵里面只要3个元素，你访问第四个就出现类似错误了。程序运行问题，你可以在M文件里面加断点，单步调试。

ChaChing

回复 3 # haorizi2012 的帖子

1. 请标题要明确!
2. Ref: 12F, 常见的程序出错问题整理 (eight)
http://forum.vibunion.com/forum/thread-46001-1-1.html