clf;nx=17; ny=17;MaxStep=120; Visc=0.1; dt=0.01; %resolution & parameters governing the run MaxIt=100; Beta=1.5; MaxErr=0.001; %parameters for SOR iteration sf=zeros(nx,ny); vt=zeros(nx,ny); w=zeros(nx,ny); h=1.0/(nx-1); t=0.0; for istep=1:MaxStep, %start the time integration for iter=1:MaxIt, %solve for the streamfunction w=sf; %by SOR iteration for i=2:nx-1; for j=2:ny-1 sf(i,j)=0.25*Beta*(sf(i+1,j)+sf(i-1,j)... +sf(i,j+1)+sf(i,j-1)+h*h*vt(i,j))+(1.0-Beta)*sf(i,j); end; end; Err=0.0; for i=1:nx,for j=1:ny, Err=Err+abs(w(i,j)-sf(i,j)); end; end; if Err <= MaxErr, break, end %stop if iteration has converged end; vt(2:nx-1,1)=-2.0*sf(2:nx-1,2)/(h*h); %set vorticity on bottom wall vt(2:nx-1,ny)=-2.0*sf(2:nx-1,ny-1)/(h*h)-2.0/h; %top wall vt(1,2:ny-1)=-2.0*sf(2,2:ny-1)/(h*h); %right wall vt(nx,2:ny-1)=-2.0*sf(nx-1,2:ny-1)/(h*h); %left wall for i=2:nx-1; for j=2:ny-1 % compute the w(i,j)=-0.25*((sf(i,j+1)-sf(i,j-1))*(vt(i+1,j)-vt(i-1,j))... % right hand side of the -(sf(i+1,j)-sf(i-1,j))*(vt(i,j+1)-vt(i,j-1)))/(h*h)... % vorticity equation +Visc*(vt(i+1,j)+vt(i-1,j)+vt(i,j+1)+vt(i,j-1)-4.0*vt(i,j))/(h*h); end; end; vt(2:nx-1,2:ny-1)=vt(2:nx-1,2:ny-1)+dt*w(2:nx-1,2:ny-1); %update the vorticity t=t+dt %print out t every step subplot(121), contour(rot90(fliplr(vt))), axis('square'); %plot vorticity subplot(122), contour(rot90(fliplr(sf))), axis('square');pause(0.01) %streamfunction end;