请教，使用MATLAB编写了一下2D泊肃叶流动，结果有误

混沌

程序和结果如下，亦上传附件

主程序：
clear;clc;
%定义基本参数
global e w NX NY U Q tau_f u rho f
Q=9;                               %D2Q9模型
NX=81; NY=27;
U=0.1;                             %流动速度
e=[0 0;1 0;0 1;-1 0;0 -1;1 1;-1 1;-1 -1;1 -1];
w(1:9)=[4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36];

%模型初始化
dx=1; dy=1;
Lx=dx*NX; Ly=dy*NY;
dt=dx;
c=dx/dt;                           %格子速度
rho0=1.0;
Re=1000;
niu=U*Lx/Re;
tau_f=3*niu+5;                     %无量纲松弛时间
disp(['tau_f=',num2str(tau_f)]);

u=zeros(NY,NX);
u(:,:,1)=U*ones(NY,NX);
u(:,:,2)=zeros(NY,NX);
rho=ones(NY,NX)*rho0;
u(1,:,=0;
u(NY,:,=0;

tic
%i是横坐标，对应NX，在网格中表示列；j是纵坐标，对应NY，在网格中表示行，按照MATLAB的书写习惯，应为u(j,i,k)
for k=1: Q
    for i=1: NX
%         u(i,j,1)=0;
%         u(i,j,2)=0;
%         rho(i,j)=rho0;
%         u(NX+1,i,1)=U;
        for j=1: NY
            f(j,i,k)=feq(k,rho(j,i),u(j,i,);
        end
    end
end

kk=0;
for n=1:300                       %演化次数
    kk=kk+1                         %输出当前演化次数
    evolution();
%     wucha=error1()
    streamslice(u(:,:,1),u(:,:,2));  
    axis equal
    drawnow
    clf
%     if wucha<1e-6
%         break
%     end
%     if mod(n,10)==2
%         clf                         %清除当前窗口图形      
%         U0=sqrt((u(:,:,1))^2+(u(:,:,2))^2);
%         imagesc(U0');
%         set(gca,'YDir','normal');   %使y坐标变为从下向上依次递增
%         colorbar      
%       axis equal off;
%     end
    toc
end

toc
xlswrite('D:\u.xlsx',u(:,:,1));
xlswrite('D:\v.xlsx',u(:,:,1));
streamslice(u(:,:,1),u(:,:,2));
axis equal off
figure
quiver(u(:,:,1),u(:,:,2));    %2表示将箭头大小扩大两倍

% U0=sqrt(u(:,:,1).^2+u(:,:,2).^2);
% figure
% contour(U0,50);


evolution
function []=evolution()
global NX NY Q tau_f u rho f e
%先计算空间内的演化后的分布函数，不包括边界上的，计算过程包括碰撞和迁移
F=zeros(NY,NX,Q,'double');  %为F预分配内存(指定数据类型,...F=zeros(NX,NY,Q,'double'))，...
                            ...减少分配内存的时间，提高运算速度
for k=1: Q
    for i=2: NX-1
        for j=2: NY-1        %碰撞
            ip=i-e(k,1);    %矩阵平移
            jp=j-e(k,2);
            F(j,i,k)=f(jp,ip,k)+(feq(k,rho(jp,ip),u(jp,ip,:))-f(jp,ip,k))/tau_f;%注意公式要写对
        end
    end
end



%计算宏观量
%不计算误差,后面直接画图
e1=[0 1 0 -1 0 1 -1 -1 1];
e2=[0 0 1 0 -1 1 1 -1 -1];
e1_0=reshape(e1,1,1,9);
e2_0=reshape(e2,1,1,9);
e1_1=repmat(e1_0,NY-2,NX-2,1);
e2_1=repmat(e2_0,NY-2,NX-2,1);

f(2:NY-1,2:NX-1,:)=F(2:NY-1,2:NX-1,:);
rho(2:NY-1,2:NX-1)=sum(f(2:NY-1,2:NX-1,:),3);        %对三维数组的第三维所有元素求和并赋值给二维平面上的点
u(2:NY-1,2:NX-1,1)=sum(e1_1.*f(2:NY-1,2:NX-1,:),3);
u(2:NY-1,2:NX-1,2)=sum(e2_1.*f(2:NY-1,2:NX-1,:),3);
u(2:NY-1,2:NX-1,1)=u(2:NY-1,2:NX-1,1)./rho(2:NY-1,2:NX-1);
u(2:NY-1,2:NX-1,2)=u(2:NY-1,2:NX-1,2)./rho(2:NY-1,2:NX-1);


% for i=2: NX
%     for j=2: NY
%         u0(j,i,1)=u(j,i,1);
%         u0(j,i,2)=u(j,i,2);  %保证在误差计算的时候，u0始终是u的前一次演化的结果
%         rho(j,i)=0;
%         u(j,i,1)=0; u(j,i,2)=0;
%         for k=1
%             f(j,i,k)=F(j,i,k);               
%             rho(j,i)=rho(j,i)+f(j,i,k);      
%             u(j,i,1)=u(j,i,1)+e(k,1)*f(j,i,k);
%             u(j,i,2)=u(j,i,2)+e(k,2)*f(j,i,k);
%         end
%         u(j,i,1)=u(j,i,1)/rho(j,i);
%         u(j,i,2)=u(j,i,2)/rho(j,i);
%     end
% end

% 边界条件（反弹边界）
%上边界
f(NY,:,[5 8 9])= f(NY,:,[3 6 7]);
%下边界
f(1,:,[3 6 7])=f(1,:,[5 8 9]);

混沌

额。。。。怎么出现了那么多表情

混沌

不好意思忘记把平衡态分布函数发上来了
function f=feq(k,rho,u)
global e w
eu=u(1)*e(k,1)+u(2)*e(k,2);
u_2=u(1)*u(1)+u(2)*u(2);

f=w(k)*rho*(1+3*eu+4.5*eu^2-1.5*u_2);
end