Recently, much attention has been paid to the AAA algorithm's capability of returning poles readily available for subsequent approximation tasks. However, we must not forget that it also returns zeros, and this example will show it is competitive with other methods for rootfinding. Here we tackle some problems previously published in the Examples-Rootfinding category.

To begin with, let's consider the Bessel function $J_0$ on the interval $[0,100]$, with its AAA approximation on a set $X$ of just 400 equispaced points:

J0 = chebfun(@(x) besselj(0,x),[0 100]);
X = linspace(0,100,400);
tic, [j0,~,~,zer] = aaa(besselj(0,X),X); toc

Elapsed time is 0.058446 seconds.

We keep only the real zeros in the interval, and measure the maximum difference from those computed by Chebfun:

LW = 'linewidth'; MS = 'markersize'; FS = 'fontsize';
zer(imag(zer)~=0) = []; zer = sort(zer(zer>=0 & zer<=100))';
format long e
max_diff = norm(zer-roots(J0)','inf')
figure, plot(J0,'c'), hold on, plot(X,j0(X),'k:'), grid on
title('Roots of Bessel function J_0 in [0,100]')
plot(zer+1i*eps,'r.',MS,15)

max_diff =
     1.612932010175427e-12

Another interesting test is how the AAA algorithm behaves with random real polynomials, intended as a linear combination of Legendre polynomials (normalized by 2-norm on $[-1,1]$) with random independent coefficients drawn from the standard normal distribution [1]. We build ten random polynomials of degree 50 on the interval with 600 evenly distributed points.

X = linspace(-1,1,600); n = 50; diff = []; tic
rng(1)
for k = 1:10
  cleg = randn(n+1,1);                      % Legendre coeffs
  ccheb = leg2cheb(cleg,'norm');            % Chebyshev coeffs
  p = chebfun(ccheb,'coeffs');
  [pa,~,~,zer] = aaa(p(X),X,'tol',1e-9);    % AAA approximation
  zer(imag(zer)~=0) = []; zer = sort(zer(zer>=-1 & zer<=1));
  diff = [diff norm(roots(p)-zer,'inf')];
end, toc

Elapsed time is 0.545659 seconds.

Here is the maximum difference between roots for each iteration, and how the last random polynomial looks like together with its zeros:

diff = diff'
clf, plot(X,p(X),'c-',LW,2), hold on, plot(X,pa(X),'k:',LW,1)
plot(zer+1i*eps,'r.',MS,15), grid on
title('Roots in [-1,1] of random polynomial, degree = 50')

diff =
     1.343369859796439e-14
     5.107025913275720e-15
     1.872002552971708e-11
     5.224598531583524e-12
     1.776356839400250e-14
     4.329869796038111e-15
     3.566591466608315e-14
     6.550315845288424e-15
     3.144184912429182e-11
     9.117706589734098e-15

Let's now find the roots in the unit disk for some complex functions. Needless to say, the AAA algorithm catches them instantly.

Z = exp(1i*linspace(0,2*pi,1000));  % 1000 points on the disk boundary
ff = {}; zerff = {};                % cells for functions and zeros

tic
[ff{1},~,~,zer] = aaa((Z-0.5i).*exp(Z),Z);
zerff{1} = zer(abs(zer)<=1); zeros = cell2mat(zerff(1))

zeros =
      2.030217387992154e-16 + 4.999999999999997e-01i

[ff{2},~,~,zer] = aaa(cosh(pi*Z),Z);
zerff{2} = zer(abs(zer)<=1); zeros = cell2mat(zerff(2))

zeros =
     -4.737691658779396e-16 - 5.000000000000004e-01i
      1.503446612679992e-15 + 5.000000000000002e-01i

[ff{3},~,~,zer] = aaa(cosh(exp(Z)).*(Z-.3).*(1+4*Z.^2),Z);
zerff{3} = zer(abs(zer)<=1); zeros = cell2mat(zerff(3))

zeros =
     -8.511597100346472e-16 - 4.999999999999998e-01i
      2.045447541696151e-16 + 5.000000000000012e-01i
      2.999999999999996e-01 - 1.006559310013407e-15i

[ff{4},~,~,zer] = aaa((Z.^3-1/8).*exp((-1-2i)*Z),Z);
zerff{4} = zer(abs(zer)<=1); zeros = cell2mat(zerff(4))
toc

zeros =
     -2.500000000000026e-01 - 4.330127018922224e-01i
      5.000000000000022e-01 + 3.988939810727542e-15i
     -2.500000000000012e-01 + 4.330127018922205e-01i
Elapsed time is 0.108164 seconds.

Here's a beautiful phaseplot image, with the zeros highlighted in each subpicture.

clf, for k = 1:4
    subplot(2,2,k), phaseplot(ff{k}), hold on
    plot(Z,'k-',LW,1), plot(zerff{k},'r.',MS,20)
end