% GNU Octave/MATLAB script extKAmNat.m belonging to the article "Grazing away the resilience of patterned ecosystems": produces model runs as in Figure 2. % Scaled extended Klausmeier model on [0,L], with periodic boundary conditions, without slope: % w_t = ew_{xx}+a(t)-w-wn^2; % n_t = n_{xx}-(m0+g)n+wn^2. Gtype = 1; %% Set type of grazing (I=1,II=2,III=3) switch Gtype case 1 fprintf('local (type I) grazing') mloc=0.225; msus=0; %% no sustained (type II) grazing mnat=0; %% no natural (type III) grazing K=0; case 2 fprintf('sustained (type II) grazing') mloc=0; %% no local (type I) grazing msus=1.5; %% m=msus mnat=0; %% no natural (type III) grazing K=0.3; case 3 fprintf('natural (type III) grazing') mloc=0; %% no local (type I) grazing msus=0; %% no sustained (type II) grazing mnat=1.5; %% m=mnat K=0.3; otherwise error('Gtype has invalid value') end %% Decrease in rainfall a a = 3; %% Starting value for a TimeStep= 0.01; %% Time step used in Euler forward dadt = -0.0001; %% Rate of increase (decrease) of rainfall EndTime = 30000; %% Amount of time after which simulation stops, thus at value a=astart+dadt*EndTime (=0 if no values are changed) %% Spatial resolution and setup spectral (Fourier) method NumberMeshPointsSpace = 4096; %% Number of spatial discretization (grid) points, equidistant L=500; %% Length of spatial interval [0,L], rescaled to [0,1] below %% Klausmeier model parameters m0=0.225; %% non-grazing senescence ediff=500; %% diffusion coefficient e %% Possible wavenumbers k=[0:round(NumberMeshPointsSpace/2),-round(NumberMeshPointsSpace/2)+1:-1]*2*pi; k=k'; %% Calculation of operators needed below, only need to be computed once. LU = [-ediff*k.^2/L^2-1,-k.^2/L^2-m0]; %% Fourier transform of linear operator [e(d/dx)^2-1,(d/dx)^2-m0], after rescaling [0,L] to [0,1] EXP = exp(LU*TimeStep); ETD = (exp(LU*TimeStep)-1)./LU; %% Numerical computation of spatially homogeneous vegetated steady state (initial condition) vegetationfunc=@(n) a*n^2/(1+n^2)-m0*n-mloc*n-msus*n/(K+n)-mnat*n^2/(K^2+n^2); vegetationstart=fzero(vegetationfunc,100); waterstart=a/(1+vegetationstart^2); u=[waterstart*ones(1,NumberMeshPointsSpace);vegetationstart*ones(1,NumberMeshPointsSpace)]; u=u'; %% Setting of noise and save options N = round(EndTime/TimeStep); %% Number of time steps NoiseEvery = 100; %% Multiplicative noise every NoiseEvery time steps SaveEvery = 1000; %% Solution is saved every SaveEvery time steps NSaveEvery=floor(N/SaveEvery); %% Total number of solutions that is to be saved jSaveEvery=1; %% Counter of how many solutions have been saved PerturbAmp=10^-3; %% Amplitude of multiplicative perturbations %% Preparation before the first iteration fftu = fft(u); umat = zeros(NumberMeshPointsSpace,NSaveEvery+1,2); %% Matrix containing all saved (intermediate) solutions umat(:,1,:)=u; %% Save initial condition %% First order exponential time differencing (ETD) method for j=1:N a = a+dadt*TimeStep; %% Update to new rainfall value meann = ones(1,NumberMeshPointsSpace)*u(:,2)/NumberMeshPointsSpace; %% Mean vegetation density gloc = mloc; %% Local (type I) grazing gsus = msus/(K+meann); %% Sustained (type II) grazing gnat = mnat*meann/(K^2+meann^2); %% Natural (type III) grazing f = [a-u(:,1).*u(:,2).^2,-(gloc+gsus+gnat)*u(:,2)+u(:,1).*u(:,2).^2]; %% Nonlinear reaction term fftf = fft(f); %% Fourier transform fftu = fftu.*EXP + fftf.*ETD; %% Approximate solution of ODE in wavenumber space by first order ETD at new time level u = real(ifft(fftu)); %% Inverse Fourier transform yields solution at new time level if mod(j,SaveEvery) == 0 %% Displays and saves once every StoreEvery times plot([0:NumberMeshPointsSpace-1],u(:,1),'Color','blue'); hold on; plot([0:NumberMeshPointsSpace-1],u(:,2),'Color','green'); hold off; title(['a=',num2str(a)]); drawnow; jSaveEvery=jSaveEvery+1; umat(:,jSaveEvery,:)=u; if max(u(:,2))<=10^-6 break; %% Model run is stopped immediately if bare desert state is reached end end if mod(j,NoiseEvery) == 0 %% Multiplicative noise with amplitude PerturbAmp is applied every NoiseEvery time steps u(:,1)=u(:,1).*(1+(rand(NumberMeshPointsSpace,1)-1/2)*PerturbAmp); u(:,2)=u(:,2).*(1+(rand(NumberMeshPointsSpace,1)-1/2)*PerturbAmp); fftu=fft(u); end end imagesc(umat(:,:,2)) %% Plot of the vegetation component