wave_eqn_neumman_2d

function [udata] = wave_eqn_neumman_2d
% wave_eqn_neumman_2d
%
% Sovles the following 2D PDE with Neumann boundary conditions.
% PDE is numerically solved on the interval [0,L] X [0,M]
%
% utt = D * (uxx + uyy)
%
% Created By: David Gross
% Date Created: Apr 18, 2006
% Date Modified: Apr 21, 2006

% -------------------------------------------------------------------------
% Grid, variables, and initial data:
% -------------------------------------------------------------------------
Nx = 50;                                   % Number of mesh points
Ny = 50;
L = 1;                                     % Domain Size
M = 1;
dx = L/Nx;
dy = M/Ny;
x = (0:Nx)*dx;                              % Vector to solve PDE on [0,L]
y = (0:Ny)*dy;
D = 1;                                     % Diffusion constant
t = 0;                                     % Starting Time
dt = 0.005; tmax = 1;                    % Time step and final time
nux = D*dt/dx^2;
nuy = D*dt/dy^2;

[xx, yy] = meshgrid(x,y);
% u0 = exp(-40*((xx-L/2).^2+(yy - M/4).^2));

u0_x = '1*(-0.5*(cos(2*pi*(x-.25)/L)-1)).^4';
u0_y = '1*(-0.5*(cos(2*pi*(y-.25)/M)-1)).^4';  % Initial Data
% u0_x = 'sin(3*pi*x/(L))';
% u0_y = 'abs(sin(pi/1.5*y/M))';

% -------------------------------------------------------------------------
% Implicit Matrices
% -------------------------------------------------------------------------

vars = {'x' 'y'};
for i = 1:length(vars)
    eval(['row = zeros(1,N' vars{i} '+1);']);
    eval(['row(1) = 1+nu' vars{i} ';']);
    eval(['row(2) = -nu' vars{i} '/2;']);
    eval([vars{i} 'mat = toeplitz(row);']);
        % Boundary Conditions
    eval([vars{i} 'mat(1, 1:3) = [-3,4,-1];']);
    eval([vars{i} 'mat(N' vars{i} '+1, N' vars{i} '-1:N' vars{i} '+1) = [1, -4, 3];']);
end

% -------------------------------------------------------------------------
% Fix initial conditions and equations
% -------------------------------------------------------------------------
% u = eval(u0);
u = eval(vectorize(u0_x))'*eval(vectorize(u0_y));
% u = u0;
u_old = u;

% -------------------------------------------------------------------------
% Plot values
% -------------------------------------------------------------------------
nplots = round(tmax/dt);
frames = 50;
nplt = floor(nplots/frames);

% -------------------------------------------------------------------------
% Put initial values in memory
% -------------------------------------------------------------------------
udata(:,:,1) = u;
tdata = t;

        clf

        set(gcf, 'Renderer', 'zbuffer');
        [xxx, yyy] = meshgrid(0:1/32:L,0:1/32:M);
        uuu =  interp2(xx, yy, u, xxx, yyy, 'cubic');
        surf(xxx, yyy, uuu);
%         surf(x, y, u);
        % view(150,30)
        zlim([-.15 1]);
        colormap(jet);
        title(['t = ' num2str(t)]);
        xlabel('y');
        ylabel('x');
        drawnow;

        pause

% -------------------------------------------------------------------------
% Solve PDE and plot results
% -------------------------------------------------------------------------
tic
for n = 1:nplots

    t = t+dt;

    [uxx uyy] = laplacian(u);

    u_new = 2*u - u_old + dt^2*D*(uxx+uyy); % Leap-frogging for 2nd time deriv.

    u_old = u;
    u = u_new;

    % printing when we want to see it
    if mod(n,nplt) == 0
        udata(:,:,n) = u; tdata = [tdata t];

        clf

        set(gcf, 'Renderer', 'zbuffer');
        [xxx, yyy] = meshgrid(0:1/32:L,0:1/32:L);
        uuu =  interp2(xx, yy, u, xxx, yyy, 'cubic');
        surf(xxx, yyy, uuu);
%         surf(y, x, u);
%         view(150,30)
        zlim([-.15 1]);
        colormap(jet);
        title(['t = ' num2str(t)]);
        xlabel('y');
        ylabel('x');
        drawnow;

    end

    % Stop if we see NaNs or numbers > 2^120
    if any(any([isnan([u]) (abs([u]) > 2^120)]))
        disp('Calculation halted due to inaccuracies'); break;
    end

end
toc
% time_elapsed = toc

% -------------------------------------------------------------------------
% Plotting
% -------------------------------------------------------------------------



% -------------------------------------------------------------------------
% Laplacian for 2D derivative calculation:
% -------------------------------------------------------------------------
function [uxx uyy] = laplacian(u)
x_N = length(u(1,:));
y_N = length(u(:,1));

vars = {'x' 'y'};
for i = 1:length(vars)
    if (eval(['mod(' vars{i} '_N,2)==0']))
        eval(['k_' vars{i} ' = sqrt(-1)*[0:' vars{i} '_N/2-1 0 -' vars{i} '_N/2+1:-1];']);
    else
        eval(['k_' vars{i} ' = sqrt(-1)*[0:(' vars{i} '_N-1)/2 0 -(' vars{i} '_N-1)/2+1:-1];']);
    end
end

uxx = [];
uyy = [];

for i=1:y_N
    ux = real(ifft(k_x .* fft(u(i,:))));
    uxx = [uxx ; real(ifft(k_x .* fft([0 ux(2:x_N-1) 0 ])))];
end

for i=1:x_N
    uy = real(ifft(k_y .* fft(u(:,i)')));
    uyy = [uyy , real(ifft(k_y .* fft([0 uy(2:y_N-1) 0 ])))'];
end