Approximaths

Padé to the Rescue

Written by

approximaths

in

Recall that in the previous posts we attempted to calculate the following Stieltjes integral using perturbation techniques:

\displaystyle I(\epsilon) = \int_{0}^{\infty} \frac{e^{-t}}{1 + \epsilon^2 t}\,dt

For \epsilon = 1, we obtain the target integral I(1) = \int_{0}^{\infty} \frac{e^{-t}}{1+t}\,dt. Expanding (1 + \epsilon^2 t)^{-1} as a formal power series and integrating term by term leads to the series:

\displaystyle \sum_{k=0}^{\infty} (-1)^{k} k! \, \epsilon^{2k}

This series is divergent for \epsilon = 1 and does not allow us to directly approximate the exact value of the integral. We can, however, approximate the exact value using Padé approximants. By making the change of variable x = \epsilon^2, the series transforms into a dense power series in x, \sum_{k=0}^{\infty} (-1)^{k} k! \, x^k. We then construct its diagonal Padé approximants [N/N](x) evaluated at x = 1.

The following R code calculates and evaluates the [5/5] Padé approximant using the pracma package (note: coefficients must be supplied in decreasing order of powers for pracma::pade):

library(pracma)

# Series coefficients a_k = (-1)^k * k! in DECREASING order of powers
k <- 10:0
p1 <- (-1)^k * factorial(k)

# Compute [5/5] Pade approximant with respect to x = epsilon^2
Q  <- pade(p1, d1 = 5, d2 = 5)
r1 <- Q$r1
r2 <- Q$r2

# Rational function evaluation
f1 <- function(x) polyval(r1, x) / polyval(r2, x)

# Evaluate at epsilon = 1 (x = 1^2 = 1)
f1(1)

# Plot graph
xs <- seq(-1, 1, length.out=100)
ys1 <- f1(xs)
plot(xs, ys1, type = "l", col="blue", xlab="x", ylab="P[5,5](x)")
grid()

In the table below, the values of successive diagonal Padé approximants [N/N](x) evaluated at x = 1 are displayed. We can observe that these values converge rapidly towards the exact value of the integral (I \approx 0.596347).

Padé approximants have therefore made it possible to evaluate a Stieltjes integral represented by a divergent series. This result can be extended to the class of all Stieltjes functions, which will be the subject of the next posts

Padé [N/N] Padé approximant Approximation at x = 1 Error
[1/1] (1 + x)/(1 + 2x) 0.666667 7.03 × 10-2
[2/2] (1 + 5x + 2x2)/(1 + 6x + 6x2) 0.615385 1.90 × 10-2
[3/3] (1 + 11x + 26x2 + 6x3)/(1 + 12x + 36x2 + 24x3) 0.602740 6.39 × 10-3
[4/4] (1 + 19x + 102x2 + 154x3 + 24x4)/(1 + 20x + 120x2 + 240x3 + 120x4) 0.598802 2.46 × 10-3
[5/5] (1 + 29x + 272x2 + 954x3 + 1044x4 + 120x5)/(1 + 30x + 300x2 + 1200x3 + 1800x4 + 720x5) 0.597383 1.04 × 10-3
[10/10] (exact rational of degree 10) 0.596379 3.15 × 10-5

Comments

Leave a comment

More posts