MATLAB: IFFT problem when approximating integrals

ifftnumerical integration

Hi all,

I have a problem when I use the ifft function to approximate an integral. For example, I want to integrate

exp(-1i*u*x)*exp(-x)

over x from 0 to infinity, for several values of u that are reals.

So I descritized and wrote

N = 2^12; %number of steps
eta = 0.01; %integration increment
lambda = 2*pi/(N*eta); %u increment
u_j = -(N-1)*lambda/2 + lambda.*(0:N-1); %list of u
x_k = eta.*(0:N-1); %grid of x for integration

Now given the ifft formula in the Matlab help, the corresponding X(k) I need to take is

f = exp(1i*u_j(end)*x_k).*exp(-x_k)*eta;

and the values of the integral for each u_j is then

sol = N*real(ifft(f));

The problem is that I do not obtain the correct approximated values (I checked with Mathematica and we know the exact solution of the integral: 1/(1+u^2) for the real part). For example, the approximation for u_j(2055)=0.9971 is sol(2055)=0.5892 instead of 0.5065 (closer to correct value 0.5). However, the correct approximation appears to be in sol(2056)=0.5065. So there is a lag of one. Is any one able to explain me why?

Best Answer

I finally found the problem. In version R2010b of Matlab, the help documentation for fft has a typo. The Nth root of unity should be defined as exp(-2i*pi/N) instead of exp(2i*pi/N).

Thus, I should use fft instead of ifft.

Related Solutions

MATLAB: Problem with Index (Dufort Frankel Scheme)

I strongly advice you to check and rethink on your code.

xl=0; xr=1;                 % x domain [xl,xr]
J = 10;                     % J: number of division for x
dx = (xr-xl) / J;           % dx: mesh size
tf = 0.1;                   % final simulation time
Nt = 50;                    % Nt: number of time steps
dt = tf/Nt;
k = 0;
t = 0;
mu = dt/(dx)^2;
if mu > 0.5         % make sure dt satisy stability condition
    error('mu should < 0.5!')
end
% Evaluate the initial conditions
x = xl : dx : xr;              % generate the grid point
% f(1:J+1) since array index starts from 1
f = sin(pi*x) + sin(2*pi*x);
% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);
u_ex = zeros(J+1,Nt) ;
% Find the approximate solution at each time step
for n = 1:Nt-1
    t = n*dt;         % current time
    
    % boundary condition at left side
    gl = sin(pi*xl)*exp(-pi*pi*t)+sin(2*pi*xl)*exp(-4*pi*pi*t);
    % boundary condition at right side
    gr = sin(pi*xr)*exp(-pi*pi*t)+sin(2*pi*xr)*exp(-4*pi*pi*t);
    if n==1    % first time step
        for j=2:J    % interior nodes

            u(j,n) = f(j) + mu*2*(f(j+1)-(f(j)+f(j))+f(j-1));
        end
        u(1,n) = gl;   % the left-end point

        u(J+1,n) = gr; % the right-end point

    else
        for j=2:J    % interior nodes
            u(j,n)=u(j,n-1)+mu*2*(u(j+1,n)-(u(j,n+1)+u(j,n-1))+u(j-1,n));
        end
        u(1,n) = gl;   % the left-end point
        u(J+1,n) = gr; % the right-end point
    end
    
    % calculate the analytic solution
    
    for j=1:J+1
        xj = xl + (j-1)*dx;
        u_ex(j,n)=sin(pi*xj)*exp(-pi*pi*t) ...
            +sin(2*pi*xj)*exp(-4*pi*pi*t);
    end
end
% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
colormap(gray);     % draw gray figure
surf(x,tt, u');     % 3-D surface plot

xlabel('x')
ylabel('t')
zlabel('u')
title('Numerical solution of 1-D parabolic equation')
figure(2)
surf(x,tt, u_ex');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Analytic solution of 1-D parabolic equation')
maxerr=max(max(abs(u-u_ex)));

MATLAB: Eulerian Position and Velocity Updates – Compressible Fluid Dynamics

Are u_j x_j acc vectors corresponding to 61 different points? Do you want to treat each point individually? In other words, do you want to find the new position,velocity, and ro of each index?

In that case you would do

for k = 1:10  
u_j = u_j +acc
x_j = x_j(k)+u_j_new
ro = m/(sqrt(pi)*sigma).*exp((-(x_i-x_j).^2)./sigma^2)
end

Related Question

Error Plotting a Multi-variable Function with Respect to a Single Variable.