Fourier transformation of a function solved by ODEs?

Question

Yuhong Jin 2019-11-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/493347-fourier-transformation-of-a-function-solved-by-odes

评论： Yuhong Jin 2019-11-28

Hi,

I was trying to fourier transform a function solved by ode15s using the code below:

function main2()
   % need function since in a script, one cannot define function vdp1
   
  [t,X]=ode15s(@rate,[0:0.01:10],[1.86;3.2;0.4]);
  N = size(X(:,3));
  T = 0.01;
  P = abs(fft(X(:,3)))./(N/2);
  P = P(2:N/2-1).^2;
  F = [2:N/2-1]'/T;
  loglog(F,P)
  
  
 
end
function dXdt=rate(t,X)
 
      if t>0 && t<0.0225
          c=0.1.*t;
      elseif t>0.0225 && t<0.045
          c=-0.1.*t+0.0045;
      else 
          c=0;
      end
      
 
  F = 10.*(1-0.5.*X(1)+log((1+c./0.02)./(1+c./0.5)));
  A = 1./(1+exp(F));
  
  dmdt= 0.1.*(1-A)-0.2.*A;
  dYdt= 0.182*17.9.*A.*(9.7-X(2))-2.2.*X(2);
  dPdt=(0.75+0.1.*(X(2)-3.2)).*(1-X(3))-(1.13-10.*(X(2)-3.2)).*X(3);
  dXdt=[dmdt;dYdt;dPdt];
end

Here the second function is not important as I run and plotted the ode solved function and this worked successfully.

Could anyone please helped to whether the fourier transform part below the ode solver in the main2 function is correct (shown above and pasted below)?

N = size(X(:,3));
  T = 0.01;
  P = abs(fft(X(:,3)))./(N/2);
  P = P(2:N/2-1).^2;
  F = [2:N/2-1]'/T;
  loglog(F,P)
  

What I am doing here is to fourier transform the third function calculated in the [t , X] matrix and plot a log-log graph.

I also omitted the first wave number, which just corresponds to the sum of all the time-course data.

I am not sure if the T value I used here, equal to the time step for ode solver is correct?

I got this graph below but it does not look right

I would be really grateful if anyone could help me.

Thanks in advance

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Ridwan Alam 2019-11-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/493347-fourier-transformation-of-a-function-solved-by-odes#answer_403501

编辑：Ridwan Alam 2019-11-26

在 MATLAB Online 中打开

Not sure why you squared the fft, but the frequency range needs to be as follows:

N = length(X(:,3));
T = 0.01;
Fs = 1/T;
P1 = abs(fft(X(:,3)))/N;
P = P1(1:N/2+1);
P(2:end-1) = 2*P(2:end-1);
F = Fs*(0:N/2)/N;
loglog(F,P)

Hope this helps!

3 个评论
显示 1更早的评论隐藏 1更早的评论

Ridwan Alam 2019-11-28

编辑：Ridwan Alam 2019-11-28

The reason I multiply the P(2:end-1) by 2: FFT generates a two sided (symmetric around f=0) spectrum (here P1) including negative frequencies. P is the one-sided part, where P(1) refers to f=0. Since P is halved in spectrum, the components for f>0 needs to be doubled to reflect the symmetric components in the negative spectrum.

The power spectrum part was not in your question, so I got confused.

I am also not sure what you mean by wave number in an FFT. Please explain and I will try to help.

Yuhong Jin 2019-11-28

B6F1AF2F-74B9-481F-9499-58E8E11118BE.jpeg

</matlabcentral/answers/uploaded_files/251072/B27786C0-289E-41E7-A641-728387941404.jpeg> This is the full question I want to solve, and the wavenumber here means frequency (I suppose), so I started from the second one.

请先登录，再进行评论。