Physics Forums (http://www.physicsforums.com/index.php)
-   Differential Equations (http://www.physicsforums.com/forumdisplay.php?f=74)
-   -   Non linear ODE system matlab fourier analysis problem (http://www.physicsforums.com/showthread.php?t=580200)

 calum.g Feb22-12 09:45 AM

Non linear ODE system matlab fourier analysis problem

Hello all.

I am studying a system and want to investigate how the frequency of y(3) varies under different conditions. However, my the fft I perform on it tells me the frequency is zero, which must be incorrect. I have tried a stack of things but can't see what the problem is. I'm relatively new to the matlab environment.

The main code:

close all; clear all;
%--------------------------------------%
Time control and initial conditions
timeinmin=600; % Time in minutes
tfin=timeinmin*60; % Time in seconds
%--------------------------------------%
Initial Levels and system Saturation
totIKK=0.08; % Total IKKn at start
totNFkB=0.08; % Total NFkB cytoplasmic at start
T=1; % Saturation of the system
yint=[totNFkB,0,0,0,0,0,totIKK,0,0,0,0,totIKK,totNFkB,T]; % Initial conditions
%--------------------------------------%
ODEs
[t,y]=ode45(@nfkb,[0 tfin],yint); % calls Model
%--------------------------------------%
Extract Data from the ODE
fs = 1; % The sampling freqency in Hz
time = 0:1/fs:tfin; % Define time of interesrt
Datapoints = tfin*fs; % # of Data points being used
yi = interp1(t,y(:,3),time);
%plot(time/60,yi,'o',t/60,y(:,3)) % test to ensure it is working
%dbstop
%--------------------------------------%
Fourier Analysis from
NFFT = 2^nextpow2(tfin); % Next power of 2 from length of y
Y = fft(yi,NFFT)/Datapoints;
f = fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('Amplitude - |Y(f)|')
% Get value of frequency
[B,IX] = sort(2*abs(Y(1:NFFT/2+1)));
BFloor=0.1; %BFloor is the minimum amplitude value (ignore small values)
Amplitudes=B(B>=BFloor) %find all amplitudes above the BFloor
Frequencies=f(IX(1+end-numel(Amplitudes):end)) %frequency of the peaks

The ODEs:

function dydt=nfkb(t,y)
parameters;
dydt=[kd1*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-ka1*y(2)*y(1)-ki1*y(1)+...
kdegc*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))+ke1f*y(3)+kdegpin*y(11);
kd1*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-ka1*y(2)*y(1)-ki2*y(2)+...
kdegf*y(2)+ktria*y(6)-kc2*y(8)*y(2);
kd1*y(5)-ka1*y(4)*y(3)+kv*ki1*y(1)-kv*ke1f*y(3);
kd1*y(5)-ka1*y(4)*y(3)+kv*ki2*y(2)-kv*ke2*y(4)-kdegf*y(4);
ka1*y(4)*y(3)-kd1*y(5)-kv*ke1c*y(5);
kitria*(y(3).^h)/((y(3).^h)+(k.^h))-kdegtia*y(6);
kp*(y(12)-y(8)-y(7))*kbA20/(kbA20+y(10)*y(14))-y(14)*ka*y(7);
y(14)*ka*y(7)-ki*y(8);
kitra*(y(3).^h)/((y(3).^h)+(k.^h))-kdegta*y(9);
ktra*y(9)-kda*y(10);
kc3*y(8)*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-kdegpin*y(11);
0;
0;
0];

The parameters:

%rates
kv=3.3; % C:N ratio
kp=0.0006; % Conversion rate of IKKi
ka=0.004; %activation rate of IKK
ki=0.003; % Usage rate of IKK
kda=0.0045; % Deg rate of A20 protein
kbA20=0.0018; % A20 concentration
kc2=0.074; % phosopho of IkB
kc3=0.37; % phospho of IkB-NFkB
kdegpin=0.1; % splitting of Ikb-NFkB due to phospho
ka1=0.5; % nuclear IkB-NFkB formation
kd1=0.0005; % dissisociation of nuclear IkB-NFkB
kdegf=0.0005; % Deg rate of nuclear IkB
kdegc=0.000022; % dissisociation of cyto IkB-NFkB
ki1=0.0026; % nuclear entry of NF-kB
ke1f=52*10^(-6); % nuclear exit of NF-kB
ke1c=0.01; % Transport of IkB-NFkB out of nuc.
ki2=0.00067; % nuclear entry of IkB
ke2=3.35*10^(-4); % Nuclear exit of IkB
h=2; % hill coefficient of IkB transcription
k=0.065; % dissasociation constant for IkB tanscription with NFkB
kitria=1.4*10^(-7); % trans rate of IkB
ktria=0.5; % translation of IkB
kdegtia=0.0003; % Deg rate of IkB RNA
kitra=1.4*10^(-7); % Transcription rate of A20 RNA
ktra=0.5; % transcription rate of A20
kdegta=0.00048; % deg rate of A20 RNA

Thank you in advance for any insights.

 All times are GMT -5. The time now is 09:00 AM.