MATLAB水波模拟源代码

functionwaterwave%WATERWAVE%2DShallowWaterModel%%Lax-Wendrofffinitedifferencemethod.%Reflectiveboundaryconditions.%Randomwaterdropsinitiategravitywaves.%Surfaceplotdisplaysheightcoloredbymomentum.%Plottitleshowst=simulatedtimeandtv=ameasureoftotalvariation.%Anexactsolutiontotheconservationlawwouldhaveconstanttv.%Lax-Wendroffproducesnonphysicaloscillationsandincreasingtv.%%Information:%Author:cuixing%qq:577737382%email:[e-mail]%Parametersn=64;%gridsizeg=9.8;%gravitationalconstantdt=0.01;%hardwiredtimestepdx=1.0;dy=1.0;nplotstep=8;%plotintervalndrops=5;%maximumnumberofdropsdropstep=500;%dropintervalD=droplet(1.5,21);%simulateawaterdrop%Initializegraphics[surfplot,top,start,stop]=initgraphics(n);%Outerloop,restarts.whileget(stop,'value')==0set(start,'value',0)H=ones(n+2,n+2);U=zeros(n+2,n+2);V=zeros(n+2,n+2);Hx=zeros(n+1,n+1);Ux=zeros(n+1,n+1);Vx=zeros(n+1,n+1);Hy=zeros(n+1,n+1);Uy=zeros(n+1,n+1);Vy=zeros(n+1,n+1);ndrop=ceil(rand*ndrops);nstep=0;%Innerloop,timesteps.whileget(start,'value')==0&&get(stop,'value')==0nstep=nstep+1;%Randomwaterdropsifmod(nstep,dropstep)==0&&nstep=ndrop*dropstepw=size(D,1);i=ceil(rand*(n-w))+(1:w);j=ceil(rand*(n-w))+(1:w);H(i,j)=H(i,j)+rand*D;end%ReflectiveboundaryconditionsH(:,1)=H(:,2);U(:,1)=U(:,2);V(:,1)=-V(:,2);H(:,n+2)=H(:,n+1);U(:,n+2)=U(:,n+1);V(:,n+2)=-V(:,n+1);H(1,:)=H(2,:);U(1,:)=-U(2,:);V(1,:)=V(2,:);H(n+2,:)=H(n+1,:);U(n+2,:)=-U(n+1,:);V(n+2,:)=V(n+1,:);%Firsthalfstep%xdirectioni=1:n+1;j=1:n;%heightHx(i,j)=(H(i+1,j+1)+H(i,j+1))/2-dt/(2*dx)*(U(i+1,j+1)-U(i,j+1));%xmomentumUx(i,j)=(U(i+1,j+1)+U(i,j+1))/2-...dt/(2*dx)*((U(i+1,j+1).^2./H(i+1,j+1)+g/2*H(i+1,j+1).^2)-...(U(i,j+1).^2./H(i,j+1)+g/2*H(i,j+1).^2));%ymomentumVx(i,j)=(V(i+1,j+1)+V(i,j+1))/2-...dt/(2*dx)*((U(i+1,j+1).*V(i+1,j+1)./H(i+1,j+1))-...(U(i,j+1).*V(i,j+1)./H(i,j+1)));%ydirectioni=1:n;j=1:n+1;%heightHy(i,j)=(H(i+1,j+1)+H(i+1,j))/2-dt/(2*dy)*(V(i+1,j+1)-V(i+1,j));%xmomentumUy(i,j)=(U(i+1,j+1)+U(i+1,j))/2-...dt/(2*dy)*((V(i+1,j+1).*U(i+1,j+1)./H(i+1,j+1))-...(V(i+1,j).*U(i+1,j)./H(i+1,j)));%ymomentumVy(i,j)=(V(i+1,j+1)+V(i+1,j))/2-...dt/(2*dy)*((V(i+1,j+1).^2./H(i+1,j+1)+g/2*H(i+1,j+1).^2)-...(V(i+1,j).^2./H(i+1,j)+g/2*H(i+1,j).^2));%Secondhalfstepi=2:n+1;j=2:n+1;%heightH(i,j)=H(i,j)-(dt/dx)*(Ux(i,j-1)-Ux(i-1,j-1))-...(dt/dy)*(Vy(i-1,j)-Vy(i-1,j-1));%xmomentumU(i,j)=U(i,j)-(dt/dx)*((Ux(i,j-1).^2./Hx(i,j-1)+g/2*Hx(i,j-1).^2)-...(Ux(i-1,j-1).^2./Hx(i-1,j-1)+g/2*Hx(i-1,j-1).^2))...-(dt/dy)*((Vy(i-1,j).*Uy(i-1,j)./Hy(i-1,j))-...(Vy(i-1,j-1).*Uy(i-1,j-1)./Hy(i-1,j-1)));%ymomentumV(i,j)=V(i,j)-(dt/dx)*((Ux(i,j-1).*Vx(i,j-1)./Hx(i,j-1))-...(Ux(i-1,j-1).*Vx(i-1,j-1)./Hx(i-1,j-1)))...-(dt/dy)*((Vy(i-1,j).^2./Hy(i-1,j)+g/2*Hy(i-1,j).^2)-...(Vy(i-1,j-1).^2./Hy(i-1,j-1)+g/2*Hy(i-1,j-1).^2));%Updateplotifmod(nstep,nplotstep)==0C=abs(U(i,j))+abs(V(i,j));%Colorshowsmomemtumt=nstep*dt;tv=norm(C,'fro');set(surfplot,'zdata',H(i,j),'cdata',C);set(top,'string',sprintf('t=%6.2f,tv=%6.2f',t,tv))drawnowendifall(all(isnan(H))),break,end%Unstable,restartendendclose(gcf)%------------------------------------functionD=droplet(height,width)%DROPLET2DGaussian%D=droplet(height,width)[x,y]=ndgrid(-1:(2/(width-1)):1);D=height*exp(-5*(x.^2+y.^2));%------------------------------------function[surfplot,top,start,stop]=initgraphics(n);%INITGRAPHICSInitializegraphicsforwaterwave.%[surfplot,top,start,stop]=initgraphics(n)%returnshandlestoasurfaceplot,itstitle,andtwouicontroltoggles.clfshgset(gcf,'numbertitle','off','name','Shallow_water')x=(0:n-1)/(n-1);surfplot=surf(x,x,ones(n,n),zeros(n,n));gridoffaxis([0101-13])caxis([-11])shadingfacetedc=(1:64)'/64;cyan=[0*ccc];colormap(cyan)top=title('Clickstart');start=uicontrol('position',[20208020],'style','toggle','string','start');stop=uicontrol('position',[120208020],'style','toggle','string','stop');