gftool.hermpade

Hermite-Padé approximants from Taylor expansion.

See [fasondini2019] for practical applications and [baker1996] for the extensive theoretical basis.

We present the example from [fasondini2019] showing the approximations. We consider the cubic root f(z) = (1 + z)**(1/3), the radius of convergence of its series is 1.

Taylor series

Obviously the Taylor series fails for z<=1 as it cannot represent a pole, but also for larger z>=1 it fails:

>>> from scipy.special import binom
>>> an = binom(1/3, np.arange(17))  # Taylor of (1+x)**(1/3)
>>> def f(z):
...     return np.emath.power(1+z, 1/3)
>>> taylor = np.polynomial.Polynomial(an)

>>> import matplotlib.pyplot as plt
>>> x = np.linspace(0, 3, num=500)
>>> __ = plt.plot(x, f(x), color='black')
>>> __ = plt.plot(x, taylor(x), color='C1')
>>> __ = plt.ylim(0, 1.75)
>>> plt.show()

(png, pdf)

../_images/gftool-hermpade-1.png

Padé approximant

The Padé approximant can be used to improve the Taylor expansion and expands the applicability beyond the radius of convergence:

>>> x = np.linspace(-3, 3, num=501)
>>> pade = gt.hermpade.pade(an, num_deg=8, den_deg=8)
>>> __ = plt.plot(x, f(x).real, color='black')
>>> __ = plt.plot(x, pade.eval(x), color='C1')
>>> __ = plt.ylim(0, 1.75)
>>> plt.show()

(png, pdf)

../_images/gftool-hermpade-2.png

The Padé approximant provides a global approximation. For negative values, however, the Padé approximant still fails, as it cannot accurately represent a branch cut. The Padé approximant is suitable for simple poles and tries to approximate the branch-cut by a finite number of poles. It is instructive to plot the error in the complex plane:

>>> y = np.linspace(-3, 3, num=501)
>>> z = x[:, None] + 1j*y[None, :]
>>> error = abs(pade.eval(z) - f(z))

>>> import matplotlib as mpl
>>> fmt = mpl.ticker.LogFormatterMathtext()
>>> __ = fmt.create_dummy_axis()
>>> norm = mpl.colors.LogNorm(vmin=1e-16, vmax=1)
>>> __ = plt.pcolormesh(x, y, error.T, shading='nearest', norm=norm)
>>> cbar = plt.colorbar(extend='both')
>>> levels = np.logspace(-15, 0, 16)
>>> cont = plt.contour(x, y, error.T, colors='black', linewidths=0.25, levels=levels)
>>> __ = plt.clabel(cont, cont.levels, fmt=fmt, fontsize='x-small')
>>> for ll in levels:
...     __ = cbar.ax.axhline(ll, color='black', linewidth=0.25)
>>> __ = plt.xlabel(r"$\Re z$")
>>> __ = plt.ylabel(r"$\Im z$")
>>> plt.tight_layout()
>>> plt.gca().set_rasterization_zorder(1.5)  # avoid excessive files
>>> plt.show()

(png, pdf)

../_images/gftool-hermpade-3.png

Away from the branch-cut, the Padé approximant is a reasonable approximation.

Square Hermite-Padé approximant

A further improvement is obtained by using the square Hermite-Padé approximant, which can represent square-root branch cuts:

>>> herm2 = gt.hermpade.Hermite2.from_taylor(an, deg_p=5, deg_q=5, deg_r=5)

>>> __ = plt.plot(x, f(x).real, color='black')
>>> __ = plt.plot(x, herm2.eval(x).real, color='C1')
>>> __ = plt.ylim(0, 1.75)
>>> plt.show()

(png, pdf)

../_images/gftool-hermpade-4.png

It nicely approximates the function almost everywhere. Let’s compare the error to the Padé approximant:

>>> __ = plt.plot(x, abs(pade.eval(x) - f(x)), label="Padé")
>>> __ = plt.plot(x, abs(herm2.eval(x) - f(x)), label="Herm2")
>>> __ = plt.yscale('log')
>>> __ = plt.legend()
>>> plt.show()

(png, pdf)

../_images/gftool-hermpade-5.png

Let’s also compare the quality of the approximants in the complex plane:

>>> error2 = np.abs(herm2.eval(z) - f(z))

>>> __, axes = plt.subplots(ncols=2, sharex=True, sharey=True)
>>> __ = axes[0].set_title("Padé")
>>> __ = axes[1].set_title("Herm2")
>>> for ax, err in zip(axes, [error, error2]):
...     pcm = ax.pcolormesh(x, y, err.T, shading='nearest', norm=norm)
...     cont = ax.contour(x, y, err.T, colors='black', linewidths=0.25, levels=levels)
...     __ = ax.clabel(cont, cont.levels, fmt=fmt, fontsize='x-small')
...     __ = ax.set_xlabel(r"$\Re z$")
...     ax.set_rasterization_zorder(1.5)
>>> __ = axes[0].set_ylabel(r"$\Im z$")
>>> plt.tight_layout()
>>> cbar = plt.colorbar(pcm, extend='both', ax=axes, fraction=0.08, pad=0.02)
>>> cbar.ax.tick_params(labelsize='x-small')
>>> for ll in levels:
...     __ = cbar.ax.axhline(ll, color='black', linewidth=0.25)
>>> plt.show()

(png, pdf)

../_images/gftool-hermpade-6.png

Note, however, the square Hermite-Padé approximant contains the ambiguity which branch to choose. The heuristic can fail and should therefore be checked.

References

fasondini2019(1,2)

Fasondini, M., Hale, N., Spoerer, R. & Weideman, J. A. C. Quadratic Padé Approximation: Numerical Aspects and Applications. Computer research and modeling 11, 1017–1031 (2019). https://doi.org/10.20537/2076-7633-2019-11-6-1017-1031

baker1996

Baker Jr, G. A. & Graves-Morris, Pade Approximants. Second edition. (Cambridge University Press, 1996).

API

Functions

hermite2(an, p_deg, q_deg, r_deg)

Return the polynomials p, q, r for the quadratic Hermite-Padé approximant.

hermite2_lstsq(an, p_deg, q_deg, r_deg[, ...])

Return the polynomials p, q, r for the quadratic Hermite-Padé approximant.

pade(an, num_deg, den_deg[, fast])

Return the [num_deg/den_deg] Padé approximant to the polynomial an.

pade_lstsq(an, num_deg, den_deg[, rcond, fix_q])

Return the [num_deg/den_deg] Padé approximant to the polynomial an.

pader(an, num_deg, den_deg[, rcond])

Robust version of Padé approximant to polynomial an.

Classes

Hermite2(p, q, r, pade)

Square Hermite-Padé approximant with branch selection according to Padé.