Could anyone check my ode45() code please?

1 view (last 30 days)
Avan Al-Saffar
Avan Al-Saffar on 1 Sep 2014
Edited: Image Analyst on 1 Sep 2014
NB
Moved the following from title to body but have no idea what the real question is...dpb
Because nu represents the mean value so there is someone who expected to got my frequency at 1!!
ORIGINAL posting picks up here again with code...end dpb edits
My system is
function dgdt=stochasticnu(t,g,S0,Sr,nu)
% The system consists of two competing species,Grass density and Woody-vegetation density with cattle-stocking rate
dgdt(1)=1.5*g(1)*(1-(S0+Sr*cos(nu*t))-0.7*g(1)-g(2));
dgdt(2)=0.03 + g(2)*(1-2*g(1)-1.03*g(2));
dgdt=[dgdt(1) dgdt(2)]';
end
My code is:
function RunStocnnu
% Finding the DFT(Discrete Fourier Transform ) using FFT procedure to get
% the spectrum plot of the system of two competing species, Grass
% density and Woody-vegetation density after adding a small perturbation
% to the cattle-stocking rate.
tinit=0; % The initial time
tf=50; % The final time
h=0.1; % Step size
tspan=tinit:h:tf; % The integration time vector
S0=.3;
Sr=.1;
g0=[0.1,0.1]; % The initial conditions
nu=1;
[t,g]=ode45(@stochasticnu,tspan,g0,[],S0,Sr,nu); % ODE solver
% Begin a new figure
figure(01)
plot(g(:,1),g(:,2),'k*-') % plotting the phase plane of the two species
xlabel('Grass density') % The x axis
ylabel('Woody-vegetation density') % The y axis
title('A trajectory of stochastic system') % The title of the figure
n=length(g); % The length of the input data vector
dt=t(end)/(length(g)-1);
Fs=1/dt;
% T=1/Fs;
NFFT =n;
y=fft(g,NFFT)/n;
y=y(1:NFFT/2);
mx=abs(y);
f=(0:NFFT/2-1)*(Fs/NFFT); % The frequency range
figure(02)
plot(f(1:50),y(1:50)); % FFT is symmetric, throw away second half
xlabel('Frequency')
figure(03)
plot(f,mx)
  2 Comments
Image Analyst
Image Analyst on 1 Sep 2014
I have no idea what the question is. One good point though is I was able to copy and paste the code and run it, and it did successfully produce these plots and didn't give any error messages.

Sign in to comment.

Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!