title: Equatorial Kelvin Wave using FDM in MATLAB
body: <p>The MATLAB code below is to simulate an equatorial Kelvin Wave on a beta plane using finite differencing. First and second timesteps are differenced using forward time centered space, the main loop has a centered time centered space scheme. The solution to the shallow water equations is periodic in x and has solid boundaries at the top and the bottom of the domain. </p>
  
  <pre><code>dt =  2000         % time step (s)
  dx = 40000           % x grid size (m)
  dy = 40000           % y grid size (m)
  ce =  10         % phase speed of the wave (m per second)
  g = 9.81;          % acceleration due to gravity (m s^-2)
  h = ce/g;         % equivalent depth (m) - not used
  
  X =  40000000        % distance in the x direction (m)
  Y = 6000000          % distance in the y direction (m)
  
  x = [0:dx:X];
  y = transpose([-Y/2:dy:Y/2]);
  y1 = y.^2;
  
  beta = 2.3e-11;     % value of the Coriolois parameter divided by the radius     of the Earth (m^-1 s^-1)
  
  N = X/(ce*dt)          % the number of timesteps to complete one latitudinal circle
  I = 1+X/dx;           % number of grid-steps in the x-direction for the u and v components
  J = 1+Y/dy;           % number of grid-steps in the y-direction for the u and v components
  
  A=ones(1,I);     
  a = beta*y; 
  B = a*A;
  
  na = [1:N];       % set up the grid for the time array  
  na2 = [2:N-1];    % set up the grid for the interiror of the time array
  
  ia = [1:I];       % set up the grid for the x array of the u and v components
  ia2 = [2:I-1];    % set up the grid for the interior for the x array of the u and v components
  ja = [1:J];       % set up the grid for the y array of the u and v components
  ja2 = [2:J-1];    % set up the grid for the interior points of the y array of the u and v components
  
  C=B;    
  
  for n=1:N
  C(:,:,N) = B;
  end
  
  b = (g*dt)/dx;  % simplification of some terms for the future equations 
  c = ((ce^2)*dt)/(g*dx);     %simplification of some terms for the future equations
  d = (g*dt)/dy; %simplification of some terms for the future equations for v component
  e = ((ce^2)*dt)/(g*dy);
  u0 = 1;   % amplitude of the perturbation zonal velocity at the Equator
  Z0 = (u0*ce)/g;
  lambda = 40000000; % typical wavelength of equatorial Kelvin waves
  k = (2*pi)/lambda; % wavenumber
  
  u(ja,ia,na) = NaN; % set up the array u of correct dimensions
  v(ja,ia,na) = NaN;
  Z(ja,ia,na) = NaN;
  
  % Initial Conditions
  
  u_i(ja,ia) = u0*exp(-(beta*y1/(2*ce)))*cos(k*x);
  v_i(ja,ia) = 0;
  Z_i(ja,ia) = Z0*exp(-(beta*y1/(2*ce)))*cos(k*x);
  
  %FTCS for timestep 1 of u
  
  u(1,1,1) = u_i(1,1)-0.5*b*(Z_i(1,2)-Z_i(1,I))+dt*B(1,1)*v_i(1,1);
  u(ja2,1,1) = u_i(ja2,1)-0.5*b*(Z_i(ja2,2)-Z_i(ja2,I))+dt*B(ja2,1).*v_i(ja2,1);
  u(J,1,1) = u_i(J,1)-0.5*b*(Z_i(J,2)-Z_i(J,I))+dt*B(J,1)*v_i(J,1);
  
  u(1,I,1) = u_i(1,I)-0.5*b*(Z_i(1,1)-Z_i(1,I-1))+dt*B(1,I)*v_i(1,I);
  u(ja2,I,1) = u_i(ja2,I)-0.5*b*(Z_i(ja2,1)-Z_i(ja2,I-    1))+dt*B(ja2,I).*v_i(ja2,I);
  u(J,I,1) = u_i(J,I)-0.5*b*(Z_i(J,1)-Z_i(J,I-1))+dt*B(J,I)*v_i(J,I);
  
  u(1,ia2,1) = u_i(1,ia2)-0.5*b*(Z_i(1,ia2+1)-Z_i(1,ia2-1))+dt*B(1,ia2).*v_i(1,ia2);
  u(ja2,ia2,1) = u_i(ja2,ia2)-0.5*b*(Z_i(ja2,ia2+1)-Z_i(ja2,ia2-1))+dt*B(ja2,ia2).*v_i(ja2,ia2);
  u(J,ia2,1) = u_i(J,ia2)-0.5*b*(Z_i(J,ia2+1)-Z_i(J,ia2-1))+dt*B(J,ia2).*v_i(J,ia2);
  
  %FTCS for timestep 1 of v (use backward difference for the top boundary and forward difference for the bottom boundary)
  
  v(J,ia,1) = 0;  % Boundary Conditions for the top of the domain
  v(1,ia,1) = 0;  % Boundary Conditions for the bottom of the domain
  
  v(1,1,1) = v_i(1,1)-d*(Z_i(2,1)-Z_i(1,1))-dt*B(1,1)*u_i(1,1); 
  v(ja2,1,1) = v_i(ja2,1)-0.5*d*(Z_i(ja2+1,1)-Z_i(ja2-1,1))-dt*B(ja2,1).*u_i(ja2,1);
  v(J,1,1) = v_i(J,1)-d*(Z_i(J,1)-Z_i(J-1,1))-dt*B(J,1)*u_i(J,1); 
  
  v(1,I,1) = v_i(1,I)-d*(Z_i(2,I)-Z_i(1,I))-dt*B(1,I)*u_i(1,I); 
  v(ja2,I,1) = v_i(ja2,I)-0.5*d*(Z_i(ja2+1,I)-Z_i(ja2-1,I))-dt*B(ja2,I).*u_i(ja2,I);
  v(J,I,1) = v_i(J,I)-d*(Z_i(J,I)-Z_i(J-1,I))-dt*B(J,I)*u_i(J,I); 
  
  v(1,ia2,1) = v_i(1,ia2)-d*(Z_i(2,ia2)-Z_i(1,ia2))-dt*B(1,ia2).*u_i(1,ia2); 
  v(ja2,ia2,1) = v_i(ja2,ia2)-0.5*d*(Z_i(ja2+1,ia2)-Z_i(ja2-1,ia2))-dt*B(ja2,ia2).*u_i(ja2,ia2);
  v(J,ia2,1) = v_i(J,ia2)-d*(Z_i(J,ia2)-Z_i(J-1,ia2))-dt*B(J,ia2).*u_i(J,ia2); 
  
  %FTCS for timestep 1 of Z (boundaries are discretised over a single spatial
  %step
  
  Z(1,1,1) = Z_i(1,1)-0.5*c*(u_i(1,2)-u_i(1,I))-e*(v_i(2,1)-v_i(1,1)); 
  Z(ja2,1,1) = Z_i(ja2,1)-0.5*c*(u_i(ja2,2)-u_i(ja2,I))-0.5*e*(v_i(ja2+1,1)-v_i(ja2-1,1)); 
  Z(J,1,1) = Z_i(J,1)-0.5*c*(u_i(J,2)-u_i(J,I))-e*(v_i(J,1)-v_i(J-1,1)); 
  
  Z(1,I,1) = Z_i(1,I)-0.5*c*(u_i(1,1)-u_i(1,I-1))-e*(v_i(2,I)-v_i(1,I)); 
  Z(ja2,I,1) = Z_i(ja2,I)-0.5*c*(u_i(ja2,1)-u_i(ja2,I-1))-0.5*e*(v_i(ja2+1,1)-v_i(ja2-1,1)); 
  Z(J,I,1) = Z_i(J,I)-0.5*c*(u_i(J,1)-u_i(J,I-1))-e*(v_i(J,I)-v_i(J-1,I)); 
  
  Z(1,ia2,1) = Z_i(1,ia2)-0.5*c*(u_i(1,ia2+1)-u_i(1,ia2-1))-e*(v_i(2,ia2)-v_i(1,ia2)); 
  Z(ja2,ia2,1) = Z_i(ja2,ia2)-0.5*c*(u_i(ja2,ia2+1)-u_i(ja2,ia2-1))-0.5*e*(v_i(ja2+1,ia2)-v_i(ja2-1,ia2)); 
  Z(J,ia2,1) = Z_i(J,ia2)-0.5*c*(u_i(J,ia2+1)-u_i(J,ia2-1))-e*(v_i(J,ia2)-v_i(J-1,ia2)); 
  
  % Define the second timestep for the u component using a FTCS scheme
  
  u(1,1,2) = u(1,1,1)-0.5*b*(Z(1,2,1)-Z(1,I,1))+dt*C(1,1,1)*v(1,1,1);
  u(ja2,1,2) = u(ja2,1,1)-0.5*b*(Z(ja2,2,1)-Z(ja2,I,1))+dt*C(ja2,1,1).*v(ja2,1,1);
  u(J,1,2) = u(J,1,1)-0.5*b*(Z(J,2,1)-Z(J,I,1))+dt*C(J,1,1)*v(J,1,1);
  
  u(1,I,2) = u(1,I,1)-0.5*b*(Z(1,1,1)-Z(1,I-1,1))+dt*C(1,I,1)*v(1,I,1);
  u(ja2,I,2) = u(ja2,I,1)-0.5*b*(Z(ja2,1,1)-Z(ja2,I-    1,1))+dt*C(ja2,I,1).*v(ja2,I,1);
  u(J,I,2) = u(J,I,1)-0.5*b*(Z(J,1,1)-Z(J,I-1,1))+dt*C(J,I,1)*v(J,I,1);
  
  u(1,ia2,2) = u(1,ia2,1)-0.5*b*(Z(1,ia2+1,1)-Z(1,ia2-1,1))+dt*C(1,ia2,1).*v(1,ia2,1);
  u(ja2,ia2,2) = u(ja2,ia2,1)-0.5*b*(Z(ja2,ia2+1,1)-Z(ja2,ia2-1,1))+dt*C(ja2,ia2,1).*v(ja2,ia2,1);
  u(J,ia2,2) = u(J,ia2,1)-0.5*b*(Z(J,ia2+1,1)-Z(J,ia2-1,1))+dt*C(J,ia2,1).*v(J,ia2,1);
  
  % Define the second timestep for the v component using a FTCS scheme
  
  v(J,ia,2) = 0;  % Boundary Conditions for the top of the domain
  v(1,ia,2) = 0;  % Boundary Cponditions for the bottom of the domain
  
  v(1,1,2) = v(1,1,1)-d*(Z(2,1,1)-Z(1,1,1))-dt*C(1,1,1)*u(1,1,1); 
  v(ja2,1,2) = v(ja2,1,1)-0.5*d*(Z(ja2+1,1,1)-Z(ja2-1,1,1))-dt*C(ja2,1,1).*u(ja2,1,1);
  v(J,1,2) = v(J,1,1)-d*(Z(J,1,1)-Z(J-1,1,1))-dt*C(J,1,1)*u(J,1,1); 
  
  v(1,I,2) = v(1,I,1)-d*(Z(2,I,1)-Z(1,I,1))-dt*C(1,I,1)*u(1,I,1); 
  v(ja2,I,2) = v(ja2,I,1)-0.5*d*(Z(ja2+1,I,1)-Z(ja2-1,I,1))-dt*C(ja2,I,1).*u(ja2,I,1);
  v(J,I,2) = v(J,I,1)-d*(Z(J,I,1)-Z(J-1,I,1))-dt*C(J,I,1)*u(J,I,1); 
  
  v(1,ia2,2) = v(1,ia2,1)-d*(Z(2,ia2,1)-Z(1,ia2,1))-dt*C(1,ia2,1).*u(1,ia2,1); 
  v(ja2,ia2,2) = v(ja2,ia2,1)-0.5*d*(Z(ja2+1,ia2,1)-Z(ja2-1,ia2,1))-dt*C(ja2,ia2,1).*u(ja2,ia2,1);
  v(J,ia2,2) = v(J,ia2,1)-d*(Z(J,ia2,1)-Z(J-1,ia2,1))-dt*C(J,ia2,1).*u(J,ia2,1); 
  
  % Define the second timestep for the Z component using a FTCS scheme
  
  Z(1,1,2) = Z(1,1,1)-0.5*c*(u(1,2,1)-u(1,I,1))-e*(v(2,1,1)-v(1,1,1)); 
  Z(ja2,1,2) = Z(ja2,1,1)-0.5*c*(u(ja2,2,1)-u(ja2,I,1))-0.5*e*(v(ja2+1,1,1)-v(ja2-1,1,1)); 
  Z(J,1,2) = Z(J,1,1)-0.5*c*(u(J,2,1)-u(J,I,1))-e*(v(J,1,1)-v(J-1,1,1)); 
  
  Z(1,I,2) = Z(1,I,1)-0.5*c*(u(1,1,1)-u(1,I-1,1))-e*(v(2,I,1)-v(1,I,1));
  Z(ja2,I,2) = Z(ja2,I,1)-0.5*c*(u(ja2,1,2)-u(ja2,I-1,1))-0.5*e*(v(ja2+1,I,1)-v(ja2-1,I,1));
  Z(J,I,2) = Z(J,I,1)-0.5*c*(u(J,1,1)-u(J,I-1,1))-e*(v(J,I,1)-v(J-1,I,1)); 
  
  Z(1,ia2,2) = Z(1,ia2,1)-0.5*c*(u(1,ia2+1,1)-u(1,ia2-1,1))-e*(v(2,ia2,1)-v(1,ia2,1)); 
  Z(ja2,ia2,2) = Z(ja2,ia2,1)-0.5*c*(u(ja2,ia2+1,1)-u(ja2,ia2-1,1))-0.5*e*(v(ja2+1,ia2,1)-v(ja2-1,ia2,1));
  Z(J,ia2,2) = Z(J,ia2,1)-0.5*c*(u(J,ia2+1,1)-u(J,ia2-1,1))-e*(v(J,ia2,1)-v(J-1,ia2,1));
  
  for na = 2:N
  % U component
  
  u(1,1,na+1) = u(1,1,na-1)-b*(Z(1,2,na)-Z(1,I,na))+dt*C(1,1,na)*v(1,1,na)*2;
  u(ja2,1,na+1) = u(ja2,1,na-1)-b*(Z(ja2,2,na)-Z(ja2,I,na))+dt*C(ja2,1,na).*v(ja2,1,na)*2;
  u(J,1,na+1) = u(J,1,na-1)-b*(Z(J,2,na)-Z(J,I,na))+dt*C(J,1,na)*v(J,1,na)*2;
  
  u(1,I,na+1) = u(1,I,na-1)-b*(Z(1,1,na)-Z(1,I-1,na))+dt*C(1,I,na)*v(1,I,na)*2;
  u(ja2,I,na+1) = u(ja2,I,na-1)-b*(Z(ja2,1,na)-Z(ja2,I-1,na))+dt*C(ja2,I,na).*v(ja2,I,na)*2;
  u(J,I,na+1) = u(J,1,na-1)-b*(Z(J,1,na)-Z(J,I-1,na))+dt*C(J,I,na)*v(J,I,na)*2;
  
  u(1,ia2,na+1) = u(1,ia2,na-1)-b*(Z(1,ia2+1,na)-Z(1,ia2-1,na))+dt*C(1,ia2,na).*v(1,ia2,na)*2;
  u(ja2,ia2,na+1) = u(ja2,ia2,na-1)-b*(Z(ja2,ia2+1,na)-Z(ja2,ia2-1,na))+dt*C(ja2,ia2,na).*v(ja2,ia2,na)*2;
  u(J,ia2,na+1) = u(J,ia2,na-1)-b*(Z(J,ia2+1,na)-Z(J,ia2-1,na))+dt*C(J,ia2,na).*v(J,ia2,na)*2;
  
  % V component
  
  v(J,ia,na) = 0;  % Boundary Conditions for the top of the domain
  v(1,ia,na) = 0;  % Boundary Cponditions for the bottom of the domain
  
  v(1,1,na+1) = v(1,1,na-1)-2*d*(Z(2,1,na)-Z(1,1,na))-dt*C(1,1,na)*u(1,1,na)*2;
  v(ja2,1,na+1) = v(ja2,1,na-1)-(d*(Z(ja2+1,1,na)-Z(ja2-1,1,na)))-dt*C(ja2,1,na).*u(ja2,1,na)*2; %*
  v(J,1,na+1) = v(J,1,na-1)-2*d*(Z(J,1,na)-Z(J-1,1,na))-dt*C(J,1,na)*u(J,1,na)*2;
  
  v(1,I,na+1) = v(1,I,na-1)-2*d*(Z(2,I,na)-Z(1,I,na))-dt*C(1,I,na)*u(1,I,na)*2;
  v(ja2,I,na+1) = v(ja2,I,na-1)-(d*(Z(ja2+1,I,na)-Z(ja2-1,I,na)))-dt*C(ja2,I,na).*u(ja2,I,na)*2; %*
  v(J,I,na+1) = v(J,I,na-1)-2*d*(Z(J,I,na)-Z(J-1,I,na))-dt*C(J,I,na)*u(J,I,na)*2;
  
  v(1,ia2,na+1) = v(1,ia2,na-1)-2*d*(Z(2,ia2,na)-Z(1,ia2,na))-dt*C(1,ia2,na).*u(1,ia2,na)*2;
  v(ja2,ia2,na+1) = v(ja2,ia2,na-1)-d*(Z(ja2-1,ia2,na)-Z(ja2-1,ia2,na))-dt*C(ja2,ia2,na).*u(ja2,ia2,na)*2;
  v(J,ia2,na+1) = v(J,ia2,na-1)-2*d*(Z(J,ia2,na)-Z(J-1,ia2,na))-dt*C(J,ia2,na).*u(J,ia2,na)*2;
  
  disp(v(60,1,na))
  
  % Z component
  
  Z(1,1,na+1) = Z(1,1,na-1)-c*(u(1,2,na)-u(1,I,na))-2*e*(v(2,1,na)-v(1,1,na));
  Z(ja2,1,na+1) = Z(ja2,1,na-1)-c*(u(ja2,2,na)-u(ja2,I,na))-e*(v(ja2+1,1,na)-v(ja2-1,1,na));
  Z(J,1,na+1) = Z(J,1,na-1)-c*(u(J,2,na)-u(J,I,na))-2*e*(v(J,1,na)-v(J-1,1,na));
  
  Z(1,I,na+1) = Z(1,I,na-1)-c*(u(1,1,na-1)-u(1,I-1,na))-2*e*(v(2,I,na)-v(1,I,na));
  Z(ja2,I,na+1) = Z(ja2,I,na-1)-c*(u(ja2,1,na)-u(ja2,I-1,na))-e*(v(ja2+1,I,na)-v(ja2+1,I,na));
  Z(J,I,na+1) = Z(J,I,na-1)-c*(u(J,1,na)-u(J,I-1,na))-2*e*(v(J,I,na)-v(J-1,I,na));
  
  Z(1,ia2,na+1) = Z(1,ia2,na-1)-c*(u(1,ia2+1,na)-u(1,ia2-1,na))-2*e*(v(2,ia2,na)-v(1,ia2,na));
  Z(ja2,ia2,na+1) = Z(ja2,ia2,na-1)-c*(u(ja2,ia2+1,na)-u(ja2,ia2-1,na))-e*(v(ja2+1,ia2,na)-v(ja2-1,ia2,na));
  Z(J,ia2,na+1) = Z(J,ia2,na-1)-c*(u(J,ia2+1,na)-u(J,ia2-1,na))-2*e*(v(J,ia2,na)-v(J-1,ia2,na));
  
  end
  </code></pre>
  
  <p>The problem I have with it is the following: once the model goes into the for loop section, the values for v start to grow at the East and West boundaries, when they should be staying in the O(10^-3 to 10^-4). And it is these values that start to impact the u and Z at around the 11th timestep and carry on further into the model. </p>
  
  <p>If I multiply the v(:,:,na-1) in the 2 relevant equations (commented in the code with a * symbol) by any constant between 0.1 and 0.5, the model produced results that I expected. </p>
  
  <p>I wish I could call that the solution to my problem, but I wonder if anyone from the community could help me with finding the actual source of my problem? </p>