Could anyone check my ode45() code please?

조회 수: 3 (최근 30일)
Avan Al-Saffar
Avan Al-Saffar 2014년 9월 1일
편집: Image Analyst 2014년 9월 1일
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