Could anyone check my ode45() code please?

3 次查看(过去 30 天)
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 个评论
dpb
dpb 2014-9-1
Please recast the question more clearly and provide an appropriate question title. I moved the long semi-question into body but can't interpret the issue well enough to try to take a stab at an answer, sorry...
Image Analyst
Image Analyst 2014-9-1
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.

请先登录,再进行评论。

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by