Evaluation-Interpolation using FFT algorithm

2020 12 21
1 个回答
更新时间：2020 12 22
6 次查看（30 天）

0 个投票

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

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

chicken vector
chicken vector 2020-12-22
编辑：chicken vector 2020-12-22
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

请先登录，再进行评论。

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

回答（1 个）

Matt J
Matt J 2020-12-22

0 个投票

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

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

chicken vector
chicken vector 2020-12-22
编辑：chicken vector 2020-12-22
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
Matt J 2020-12-22
Ran in:
Ran in:
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
chicken vector 2020-12-22
I used 'roots' aswell and appears to have very good performances until now.
Thank you Matt for your help.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Polynomials 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by