Could anyone check my ode45() code please?

2 Ansichten (letzte 30 Tage)

Ältere Kommentare anzeigen

Avan Al-Saffar am 1 Sep. 2014

0
Verknüpfen

Direkter Link zu dieser Frage

https://de.mathworks.com/matlabcentral/answers/153051-could-anyone-check-my-ode45-code-please

Bearbeitet: Image Analyst am 1 Sep. 2014

In MATLAB Online öffnen

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 Kommentare
Keine anzeigenKeine ausblenden

dpb am 1 Sep. 2014

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 am 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.