数学建模美赛集训的时候要用到一个海面模拟，分享一下海面模拟的MATLAB代码
先贴一下结果图：
下面是源代码~~~
1 function waterwave 2 3 4 n = 64; % grid size 5 g = 9.8; % gravitational constant 6 dt = 0.01; % hardwired timestep 7 dx = 1.0; 8 dy = 1.0; 9 nplotstep = 8; % plot interval 10 ndrops = 5; % maximum number of drops 11 dropstep = 500; % drop interval 12 D = droplet(1.5,21); % simulate a water drop 13 14 % Initialize graphics 15 16 [surfplot,top,start,stop] = initgraphics(n); 17 18 % Outer loop, restarts. 19 20 while get(stop,'value') == 0 21 set(start,'value',0) 22 23 H = ones(n+2,n+2); U = zeros(n+2,n+2); V = zeros(n+2,n+2); 24 Hx = zeros(n+1,n+1); Ux = zeros(n+1,n+1); Vx = zeros(n+1,n+1); 25 Hy = zeros(n+1,n+1); Uy = zeros(n+1,n+1); Vy = zeros(n+1,n+1); 26 ndrop = ceil(rand*ndrops); 27 nstep = 0; 28 29 % Inner loop, time steps. 30 31 while get(start,'value')==0 && get(stop,'value')==0 32 nstep = nstep + 1; 33 34 % Random water drops 35 if mod(nstep,dropstep) == 0 && nstep <= ndrop*dropstep 36 w = size(D,1); 37 i = ceil(rand*(n-w))+(1:w); 38 j = ceil(rand*(n-w))+(1:w); 39 H(i,j) = H(i,j) + rand*D; 40 end 41 42 % Reflective boundary conditions 43 H(:,1) = H(:,2); U(:,1) = U(:,2); V(:,1) = -V(:,2); 44 H(:,n+2) = H(:,n+1); U(:,n+2) = U(:,n+1); V(:,n+2) = -V(:,n+1); 45 H(1,:) = H(2,:); U(1,:) = -U(2,:); V(1,:) = V(2,:); 46 H(n+2,:) = H(n+1,:); U(n+2,:) = -U(n+1,:); V(n+2,:) = V(n+1,:); 47 48 % First half step 49 50 % x direction 51 i = 1:n+1; 52 j = 1:n; 53 54 % height 55 Hx(i,j) = (H(i+1,j+1)+H(i,j+1))/2 - dt/(2*dx)*(U(i+1,j+1)-U(i,j+1)); 56 57 % x momentum 58 Ux(i,j) = (U(i+1,j+1)+U(i,j+1))/2 - ... 59 dt/(2*dx)*((U(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ... 60 (U(i,j+1).^2./H(i,j+1) + g/2*H(i,j+1).^2)); 61 62 % y momentum 63 Vx(i,j) = (V(i+1,j+1)+V(i,j+1))/2 - ... 64 dt/(2*dx)*((U(i+1,j+1).*V(i+1,j+1)./H(i+1,j+1)) - ... 65 (U(i,j+1).*V(i,j+1)./H(i,j+1))); 66 67 % y direction 68 i = 1:n; 69 j = 1:n+1; 70 71 % height 72 Hy(i,j) = (H(i+1,j+1)+H(i+1,j))/2 - dt/(2*dy)*(V(i+1,j+1)-V(i+1,j)); 73 74 % x momentum 75 Uy(i,j) = (U(i+1,j+1)+U(i+1,j))/2 - ... 76 dt/(2*dy)*((V(i+1,j+1).*U(i+1,j+1)./H(i+1,j+1)) - ... 77 (V(i+1,j).*U(i+1,j)./H(i+1,j))); 78 % y momentum 79 Vy(i,j) = (V(i+1,j+1)+V(i+1,j))/2 - ... 80 dt/(2*dy)*((V(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ... 81 (V(i+1,j).^2./H(i+1,j) + g/2*H(i+1,j).^2)); 82 83 % Second half step 84 i = 2:n+1; 85 j = 2:n+1; 86 87 % height 88 H(i,j) = H(i,j) - (dt/dx)*(Ux(i,j-1)-Ux(i-1,j-1)) - ... 89 (dt/dy)*(Vy(i-1,j)-Vy(i-1,j-1)); 90 % x momentum 91 U(i,j) = U(i,j) - (dt/dx)*((Ux(i,j-1).^2./Hx(i,j-1) + g/2*Hx(i,j-1).^2) - ... 92 (Ux(i-1,j-1).^2./Hx(i-1,j-1) + g/2*Hx(i-1,j-1).^2)) ... 93 - (dt/dy)*((Vy(i-1,j).*Uy(i-1,j)./Hy(i-1,j)) - ... 94 (Vy(i-1,j-1).*Uy(i-1,j-1)./Hy(i-1,j-1))); 95 % y momentum 96 V(i,j) = V(i,j) - (dt/dx)*((Ux(i,j-1).*Vx(i,j-1)./Hx(i,j-1)) - ... 97 (Ux(i-1,j-1).*Vx(i-1,j-1)./Hx(i-1,j-1))) ... 98 - (dt/dy)*((Vy(i-1,j).^2./Hy(i-1,j) + g/2*Hy(i-1,j).^2) - ... 99 (Vy(i-1,j-1).^2./Hy(i-1,j-1) + g/2*Hy(i-1,j-1).^2)); 100 101 % Update plot 102 if mod(nstep,nplotstep) == 0 103 C = abs(U(i,j)) + abs(V(i,j)); % Color shows momemtum 104 t = nstep*dt; 105 tv = norm(C,'fro'); 106 set(surfplot,'zdata',H(i,j),'cdata',C); 107 set(top,'string',sprintf('t = %6.2f, tv = %6.2f',t,tv)) 108 drawnow 109 end 110 111 if all(all(isnan(H))), break, end % Unstable, restart 112 end 113 end 114 close(gcf) 115 116 % ------------------------------------ 117 118 function D = droplet(height,width) 119 % DROPLET 2D Gaussian 120 % D = droplet(height,width) 121 [x,y] = ndgrid(-1:(2/(width-1)):1); 122 D = height*exp(-5*(x.^2+y.^2)); 123 124 % ------------------------------------ 125 126 function [surfplot,top,start,stop] = initgraphics(n); 127 % INITGRAPHICS Initialize graphics for waterwave. 128 % [surfplot,top,start,stop] = initgraphics(n) 129 % returns handles to a surface plot, its title, and two uicontrol toggles. 130 131 clf 132 shg 133 set(gcf,'numbertitle','off','name','Shallow_water') 134 x = (0:n-1)/(n-1); 135 surfplot = surf(x,x,ones(n,n),zeros(n,n)); 136 grid off 137 axis([0 1 0 1 -1 3]) 138 caxis([-1 1]) 139 shading faceted 140 c = (1:64)'/64; 141 cyan = [0*c c c]; 142 colormap(cyan) 143 top = title('Click start'); 144 start = uicontrol('position',[20 20 80 20],'style','toggle','string','start'); 145 stop = uicontrol('position',[120 20 80 20],'style','toggle','string','stop');