Fisher-Kolmogorov1D

Fisher-Kolmogorov1D.m

%  Use of 2D FFT for Fourier-spectral solution of
%
%      ut = u_xx + u(1-u) , 0< x < 50
%
%  where
% u0 = exp(-r^2/(2*sig^2)),  
%      r^2 = (x-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*32;      % No. of Fourier modes...should be a power of 2
L = 50;       % Domain size (assumed square)
sig = 1;   % Characteristic width of u0

% Vector of wavenumbers
k = (2*pi/L)*[0:(N/2-1) (-N/2):(-1)]; 

% Laplacian matrix acting on the wavenumbers
% Kluge to avoid division by zero of (0,0) waveno.
% of fhat (this waveno. should be zero anyway!) 
delsq = -k.^2;
delsq(1,1) = 1;
            
% Construct RHS f(x,y) at the Fourier gridpoints
h = L/N;     % Grid spacing
x = (0:(N-1))*h ;
rsq = (x-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:250
    uhat = fft2(u);
    u = u+dt*real(ifft2(uhat./delsq))+dt*u.*(1-u);
    u = u - u(1);% Specify arbitrary constant by forcing corner u = 0.
end

%  Plot out solution in interior
plot(x,u)
xlabel('x')
ylabel('u')
title('Fourier spectral method for 1D Fisher-Kolmogorov Eqn')