网格数值计算中周期性边界条件的处理

       for i=1:N
         ip(i)=i+1;
         im(i)=i-1;
       end
       im(1)=N;
       ip(N)=1;
%
%      begin simulation loop
       for iter=1:nit
           for i=1:N
             for j=1:N
              dmx=(phi(i,j)-phi(im(i),j))/h;                     % x backward difference
              dpx=(phi(ip(i),j)-phi(i,j))/h;                     % x forward difference
              dmy=(phi(i,j)-phi(i,im(j)))/h;                     % y backward difference
              dpy=(phi(i,ip(j))-phi(i,j))/h;                     % y forward difference
              convx=max(u(i,j),0)*dmx+min(u(i,j),0)*dpx;        % if u(i,j)>0, use backward difference dmx, or else use dpx
              convy=max(v(i,j),0)*dmy+min(v(i,j),0)*dpy;
              phin(i,j)=phi(i,j)-(convx+convy)*dt;         % advance by dt
            end
           end 
           phi=phin;                                             % update
%
%         Plotting
          contour(X,Y,phi,[0,0],‘r‘);
          axis([-1 1 -1 1])
          axis(‘square‘)
          pause(.001)
       end