Derivatives using fft in Matlab

1.6k Views Asked by At

I have been used matlab for two weeks and I start to appreciate its power. At the moment I am studying Fourier series in order to solve PDEs.

The exercise below asks me to find the derivate of a function using the fast fourier transform function in Matlab.

I am not able to understand where I am wrong. Thank you!

l= 8;    
n = 100;    
x = (linspace(0,l,n))'; 

f = @(x) x/l.*sin(10*pi*x/l);    
d_f_ex = @(x) (1/l)*sin(10*pi*x/l)+ 10*pi/(l^2)*x.*cos(10*pi*x/l);    
fval = f(x);    
ft = (1/n)*fft(f(x));    
kk = [(0:floor(n/2))';(-ceil(n/2)+1:-1)'];    
ft_d = ft.* exp(1i*2*pi*kk/l);    
d_f_rec = ifft(n*ft_d);    
%%
plot (x, d_f_ex(x), 'or', x, d_f_rec, 'ob');
1

There are 1 best solutions below

9
On BEST ANSWER

In fact, I am ashamed because the first version of this answer was worst than unclear, it was not exact.

Having found meanwhile a very good answer (https://math.stackexchange.com/q/741238) I refer to it. I have (slightly) adapted the code that can be found there with some supplementary explanations:

clear all;close all;hold on;
a = 0;
b = 8*pi;
N = 4096;
dx =(b-a)/N;
t = a + dx*(0:N-1);
f = t.*sin(t);% your function
df1= sin(t)+t.*cos(t);% is derivative
plot(t,df1,'b')
fftx = fft(f);
k = (2*pi/(b-a))*[0:N/2-1, 0, -N/2+1:-1];% "swapping"
dffft = i*k.*fftx;% multiplication by ik
df2 = ifft(dffft);
plot(t,real(df2),'r'); % graphical "perfect coincidence"

Remark 1: It wasn't necessary, on the last line to take the real part. But it acts as a firewall in case imaginary parts come in for any reason. Besides, one can check that imag(df) is identically zero ; said otherwise, plot(t,imag(df),'r') gives ... the horizontal axis.

Remark 2 : It is interesting to plot the error df1-real(df2). It is rather small, but not negligeable, in the central part and explodes at the two ends. An interesting error analysis can be found in this question and its last answer : (https://math.stackexchange.com/q/869197).

Remark 3: As a general rule, it's better to express the mathematical expressions like $x.*\sin(x)$ without "implementation details".