2D-Fisher-Kolmogorov

Fisher-Kolmogorov.m

    %  Use of 2D FFT for Fourier-spectral solution of
    %
    %      ut = u_xx + u_yy 0 + u(1-u) , 0< x,y < 1
    %
    %  where
    % u0 = exp(-r^2/(2*sig^2)),  
    %      r^2 = (x-0.5)^2 + (y-0.5)^2, sig = 1
    %
    %  with periodic BCs on u in x,y, using N = 16 modes in each direction.
    %  Script makes a surface plot of u at the Fourier grid points.
    
    N = 32;      % No. of Fourier modes...should be a power of 2
    L = 10;       % Domain size (assumed square)
    sig = 1;   % Characteristic width of u0
    
    k = (2*pi/L)*[0:(N/2-1) (-N/2):(-1)]; % Vector of wavenumbers
    [KX KY]  = meshgrid(k,k); % Matrix of (x,y) wavenumbers corresponding
                                % to Fourier mode (m,n)
    delsq = -(KX.^2 + KY.^2); % Laplacian matrix acting on the wavenumbers
    delsq(1,1) = 1;           % Kluge to avoid division by zero of (0,0) waveno.
                                % of fhat (this waveno. should be zero anyway!)
    
    % Construct RHS f(x,y) at the Fourier gridpoints
    h = L/N;     % Grid spacing
    x = (0:(N-1))*h ;
    y = (0:(N-1))*h;
    [X Y] = meshgrid(x,y);
    rsq = (X-0.5*L).^2 + (Y-0.5*L).^2;
    sigsq = sig^2;
    u0 = exp(-rsq/(2*sigsq)).*1/(sigsq^2);
    
    %  Spectral inversion of Laplacian
    
    u = u0;
    dt = 0.01
    for it =0:100
        uhat = fft2(u);
        u = u+dt*real(ifft2(uhat./delsq))+dt*u.*(1-u);
        %u = u - u(1,1);   % Specify arbitrary constant by forcing corner u = 0.
    end
    
    %  Plot out solution in interior
    surf(X,Y,u)
    xlabel('x')
    ylabel('y')
    zlabel('u')
    title('Fourier spectral method for 2D Fisher-Kolmogorov Eqn')