haorizi2012 发表于 2011-3-30 14:32

请教matlab程序错误

本帖最后由 haorizi2012 于 2011-3-30 14:33 编辑

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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= ;
cy= ;

t   = ones(1,8)*alpha/(2*pi);       % weight
az= ;      
= 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,:,:), );
       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 发表于 2011-3-30 15:24

回复 1 # haorizi2012 的帖子

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

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

回复 2 # mhkmars 的帖子

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

meiyongyuandeze 发表于 2011-3-31 09:11

haorizi2012 发表于 2011-3-30 19:03 static/image/common/back.gif
回复 2 # mhkmars 的帖子

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

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

雨人 发表于 2011-3-31 23:24

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

ChaChing 发表于 2011-4-1 01:00

本帖最后由 ChaChing 于 2011-4-1 01:02 编辑

回复 3 # haorizi2012 的帖子

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

页: [1]
查看完整版本: 请教matlab程序错误