Evaluation-Interpolation using FFT algorithm

I'm trying to develop a FFT algorithm for evaluation-interpolation of polynomials.

I tried the simple function where the coefficients are expressed as but only the DFT seems to work. I've spent quite some time on this and I can't make it work. Any suggestions?

f = @(x) x^3;

Pf = [1 , 0 , 0 , 0];
yf = FFT(Pf,1);

y = FFT(yf,2)

function y = FFT(P,k)

% k = 1 -> DFT
% k = 2 -> IDFT

N = length(P);
omega = exp(2*pi*1i/N);

if k == 1
    l = 1;
    p = 1;
elseif k == 2
    l = 1/N;
    p = -1;
end

if N == 1
    y = P;
else
    n = N/2;
    P_e = P(2:2:end);
    P_o = P(1:2:end);
    y_e = FFT(P_e,k);
    y_o = FFT(P_o,k);
    y = zeros(N,1);
    for j = 1 : N/2
        y(j) = y_e(j) + (l*omega^(p*(j-1)))*y_o(j);
        y(j+n) = y_e(j) - (l*omega^(p*(j-1)))*y_o(j);
    end
end
end

[chicken vector]

For anyone having the same problem, below there's the fixed code for IFFT. I'm having some issues on dividing by N inside the recursive function, so it is done outside.

P = [%vector of the evaluations];
N = length(P);
y = IFFT(P)/N;

function y = IFFT(P)

% This works only if N = 2^k

N = length(P);
n = N/2;
omega = exp(-2*pi*1i/N);

if N == 1
    y = P;
else
    P_e = P(1:2:end);
    P_o = P(2:2:end);
    y_e = IFFT(P_e);
    y_o = IFFT(P_o);
    y = zeros(N,1);
    for j = 1 : n
        y(j) = y_e(j) + omega^(j-1)*y_o(j);
        y(j+n) = y_e(j) - omega^(j-1)*y_o(j);
    end
end
end

Answer 1

[Matt J]

A highly impractical thing to do. If you know the coefficients of the polynomial, you should just use polyval().

However, if you must use FFT interpolation, then interpft() will readily do it,

https://www.mathworks.com/help/matlab/ref/interpft.html

[chicken vector]

My problem is that I have an array of function handles and the determinant of this array is a 15th degree polynomial. I need to find the roots of this polynomial, but first I need to find the polynomial. Using a symbolic variable is computationally inefficient so I have to compute the coefficients of this polynomial, and FFT algorithm is the best option since it is optimized when the degree of the polynomial is of order .

[Matt J]

Finding the roots of a 15th order polynomial can be highly unstable numerically, e.g.,

rTrue=sort((rand(1,15))*5);

coeffsTrue=poly(rTrue), %true coefficients

coeffsTrue = 1×16

0.0000 -0.0000 0.0005 -0.0048 0.0297 -0.1319 0.4312 -1.0507 1.9172 -2.6054 2.5973 -1.8477 0.8961 -0.2745 0.0461 -0.0030

●

coeffs=coeffsTrue+[0,randn(1,15)]*1e-6*max(coeffsTrue), %add small errors to coefficients

coeffs = 1×16

0.0000 -0.0000 0.0005 -0.0048 0.0297 -0.1319 0.4312 -1.0507 1.9172 -2.6054 2.5973 -1.8477 0.8961 -0.2745 0.0461 -0.0030

●

rTrue, %true roots

rTrue = 1×15

0.1598 0.4384 0.6582 1.3390 1.5456 1.7830 2.1863 2.2286 2.2790 2.6051 2.9448 3.0386 3.6676 4.1255 4.5711

●

r=sort(real( roots(coeffs) )).' %calculated roots

r = 1×15

0.1596 0.4403 0.6541 1.0277 1.0277 1.1391 1.1391 1.2642 1.2642 1.4859 1.4859 2.3075 2.3075 8.5793 8.5793

[chicken vector]

I used 'roots' aswell and appears to have very good performances until now.

Thank you Matt for your help.