Approximaths

Approximaths

  • Stieltjes functions II

    In the previous post, we studied the direct problem of constructing a Stieltjes function from a given weight function and explored the associated Stieltjes series and Padé approximants. Now, let’s consider the opposite problem and first define the weight function. In this case, we choose:

    \displaystyle \rho(t)=\begin{cases}1 & \text{if } 0 < t < 1 \\ 0 & \text{otherwise}\end{cases}

    this implies that

    \displaystyle a_n = \int_{0}^{\infty} \rho(t) t^n \, dt
    \displaystyle = \int_{0}^{1} t^n \, dt = \frac{1}{n+1}

    This implies that the corresponding Stieltjes series is (see previous post):

    \displaystyle \sum_{n=0}^{\infty} \frac{(-1)^n x^n}{n+1}

    (which has radius of convergence 1 and coincides with f(x) for |x|<1) and the Stieltjes function is:

    \displaystyle f(x) = \frac{\log(1+x)}{x}

    (where \log denotes the principal branch, with cut along (-\infty,-1]). Thus f is analytic in \mathbb{C}\setminus(-\infty,-1].

    Properties of Stieltjes functions:

    • f(x) is analytic in the cut complex plane \mathbb{C}\setminus(-\infty,0] (in the present example the natural cut is actually shorter: \mathbb{C}\setminus(-\infty,-1]).

    • f(x)\ge 0 for all x>0.

    • f(x) \to 0 as |x| \to \infty in the cut complex plane, under appropriate assumptions on \rho(t).

    • f(x) \sim \sum_{n=0}^{\infty} (-1)^n a_n x^n as |x| \to 0 in the cut complex plane.

    • -f(x) is Herglotz (i.e., \mathrm{Im}(-f(x)) has the same sign as \mathrm{Im}(x)).

    Under appropriate additional hypotheses, a converse characterization also holds: a function satisfying the corresponding analyticity, asymptotic, positivity, and Herglotz properties is a Stieltjes function. It can be proved rigorously that appropriate Padé approximants of a Stieltjes series converge. This is a particularly useful result (which we have already illustrated with examples), as it implies that the Padé approximants constructed from the Stieltjes series converge to the Stieltjes function for x in the cut complex plane, under appropriate hypotheses.

  • Stieltjes function I

    The function f(x) is a Stieltjes function if it has the following form:

    \displaystyle f(x) = \int_{0}^{\infty} \frac{\rho(t)}{1+xt} \, dt

    where \rho(t) is a ‘weight function’ with \rho(t) \geq 0 for all 0 \leq t < \infty and the moments of \rho(t) exist. The moments are defined as:

    \displaystyle a_n = \int_{0}^{\infty} \rho(t) t^n \, dt

    Using the geometric series expansion of \frac{1}{1 + xt}, the Stieltjes function has the following asymptotic expansion:

    \displaystyle f(x) = \int_{0}^{\infty} \frac{\rho(t)}{1+xt} \, dt \sim \int_{0}^{\infty} \sum_{n=0}^{\infty} (-1)^n x^n t^n \rho(t) \, dt
    \displaystyle = \sum_{n=0}^{\infty} (-1)^n \left( \int_{0}^{\infty} \rho(t) t^n \, dt \right) x^n
    \displaystyle = \sum_{n=0}^{\infty} (-1)^n a_n x^n

    This series is called a ‘Stieltjes series’ and is often divergent since the moments a_n can grow sufficiently rapidly for the series to have zero radius of convergence. Let begin the exploration of Stieltjes functions with an example. Let set \rho(t) = e^{-t}. We have:

    \displaystyle f(x) = \int_{0}^{\infty} \frac{e^{-t}}{1+xt} \, dt
    \displaystyle \sim \sum_{n=0}^{\infty} (-1)^n \left( \int_{0}^{\infty} e^{-t} t^n \, dt \right) x^n
    \displaystyle = \sum_{n=0}^{\infty} (-1)^n n! x^n

    So \rho(t) = e^{-t} \implies a_n = n!.

    Stieltjes functions and their associated series offer a classic illustration of how divergent asymptotic expansions can still encode precise analytic information, a theme that continues to play a central role in asymptotic analysis and resurgence theory.

  • Padé to the Rescue

    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

  • Stieltjes integral II

    If we give \epsilon the value 1 for the series from the previous post and draw up a table of the values of the integral as a function of n, we can see that the values first converge towards the exact value and then move away from it (see table below). This reminds us of the “astronomers’ series” described by H. Poincaré. We also note that we cannot approach the exact value with arbitrary precision using this method.

    n Approximation Error
    0 1 0.4036
    1 0 -0.5963
    2 2 1.4037
    3 -4 -4.5963
    4 20 19.4037
    5 -100 -100.5963

    Given such a series, a natural question is: Up to what value of n should we sum in order to get as close as possible to the exact value ? In order to determine the ideal value of n, we must first note the following facts concerning the remainder of the series above:

    \displaystyle I_{n+1}(\epsilon) := \int_{0}^{\infty} \frac{e^{-t} t^{n+1}}{1 + \epsilon^2 t} dt
    \displaystyle |I_{n+1}(\epsilon)| > \int_{1}^{\infty} \frac{e^{-t} t^{n+1}}{1 + \epsilon^2 t} dt > \int_{1}^{\infty} \frac{e^{-t} t^{n+1}}{t + \epsilon^2 t} dt

    Defining t := 1+s (\implies dt = ds), we have:

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

    We observe that:

    \displaystyle \frac{e^{-1}}{1 + \epsilon^2} \int_{0}^{\infty} e^{-s} (1+s)^n ds > \frac{e^{-1}}{1 + \epsilon^2} \int_{0}^{\infty} e^{-s} s^n ds
    \displaystyle = \frac{e^{-1}}{1 + \epsilon^2}n!

    Therefore:

    \displaystyle |I_{n+1}(\epsilon)| > \frac{e^{-1}}{1 + \epsilon^2}n!
    \displaystyle \epsilon^{2n+2}|I_{n+1}(\epsilon)| > \epsilon^{2n+2}\frac{e^{-1}}{1 + \epsilon^2}n!

    Based on the above calculations, we can now define a lower bound (L_{n+1}(\epsilon)) and an upper bound (U_{n+1}(\epsilon)) for the error made when evaluating the Stieltjes integral using the sum derived above.

    \displaystyle L_{n+1}(\epsilon) := \epsilon^{2n+2}\frac{e^{-1}}{1 + \epsilon^2}n!
    \displaystyle < \left|I_{\epsilon} - \sum_{k=0}^{n} (-1)^k \epsilon^{2k} k!\right|
    \displaystyle < \epsilon^{2n+2} (n+1)! =: U_{n+1}(\epsilon)

    For \epsilon = 1, we have L_{n+1}(\epsilon) = \frac{e^{-1}}{2}n! and U_{n+1}(\epsilon) = (n+1)!. The approximation, error, and corresponding upper and lower bounds for the error are shown for \epsilon = 1 in the table below. In this table, we can clearly see that the approximations are getting worse and worse, while the lower and upper bounds are getting larger and larger. This makes the perturbative method, as implemented in this example, unattractive.

    n Approximation |Error| L_{n+1}(\epsilon) U_{n+1}(\epsilon)
    0 1 0.4036 0.18393 1
    1 0 0.5963 0.18393 2
    2 2 1.4037 0.36788 6
    3 -4 4.5963 1.10364 24
    4 20 19.4037 4.41455 120
    5 -100 100.5963 22.0728 720

  • Stieltjes integral I

    In this post we would like to illustrate the calculation of an integral using perturbation techniques. We would like to calculate the following Stieltjes integral:

    \displaystyle I = \int_{0}^{\infty} \frac{e^{-t}}{1+t}dt

    This integral can be evaluated exactly (I = 0.5963), but we’ll use perturbation techniques to compare the results obtained using these approximation techniques with the exactly calculated value. In order to do so, we will insert the parameter \epsilon as follows:

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

    For \epsilon = 0 we have I(\epsilon) = 1. This is the solution to the unperturbed problem. We observe that:

    \displaystyle \frac{1}{1-x} = 1 + x + x^2 + ...
    \displaystyle = 1 + x + x^2 + ... + x^n + x^{n+1}(1 + x + x^2 + ...)
    \displaystyle = 1 + x + x^2 + ... + x^n + \frac{x^{n+1}}{1-x}

    Defining x := -\epsilon^2 t and using \int_{0}^{\infty} e^{-t} t^n dt = n! in the formula above we obtain:

    \displaystyle \frac{1}{1 +\epsilon^2t} = 1 -\epsilon^2 t + \epsilon^4 t^2 + ... + (-1)^n \epsilon^{2n} t^{n} + \frac{(-1)^{n+1} \epsilon^{2n+2} t^{n+1}}{1+\epsilon^2 t}
    \displaystyle I(\epsilon) = \int_{0}^{\infty} e^{-t}(1 -\epsilon^2 t + \epsilon^4 t^2 + ... + (-1)^n \epsilon^{2n} t^{n} + \frac{(-1)^{n+1} \epsilon^{2n+2} t^{n+1}}{1 + \epsilon^2 t}) dt
    \displaystyle = \int_{0}^{\infty} (e^{-t} -\epsilon^2 e^{-t} t + \epsilon^4 e^{-t} t^2 + ... + (-1)^n \epsilon^{2n} e^{-t} t^{n} + \frac{(-1)^{n+1} \epsilon^{2n+2} e^{-t} t^{n+1}}{1 + \epsilon^2 t}) dt
    \displaystyle = 1 - \epsilon^2 + \epsilon^4 2! + ... + (-1)^{n}\epsilon^{2n} n! + (-1)^{n+1} \epsilon^{2n+2} I_{n+1}(\epsilon)

    where I_{n+1}(\epsilon) := \int_{0}^{\infty} \frac{e^{-t} t^{n+1}}{1 + \epsilon^2 t} dt. Finally we have:

    \displaystyle I(\epsilon) = \sum_{k=0}^{n} (-1)^{k} \epsilon^{2k} k! + (-1)^{n+1} \epsilon^{2n+2} I_{n+1}(\epsilon)

  • Quantum Mechanical Interpretation

    The differential equation presented in the previous posts

    \displaystyle y'' + (1 - x)y = 0

    appears in quantum mechanics when a particle is subjected to a linear potential. The time-independent Schrödinger equation is

    \displaystyle -\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + V(x)\psi(x) = E\psi(x)

    Consider a linear potential of the form

    \displaystyle V(x) = V_0 - Fx

    where F is a constant force. The Schrödinger equation then becomes

    \displaystyle -\frac{\hbar^2}{2m}\psi''(x) + \bigl(V_0 - Fx\bigr)\psi(x) = E\psi(x)

    Rearranging the terms gives

    \displaystyle \psi''(x) + \frac{2m}{\hbar^2}\bigl(E - V_0 + Fx\bigr)\psi(x) = 0

    We introduce the turning point x_0 (in quantum mechanics, it marks the transition between an oscillatory (classically allowed) region and an exponentially decaying (classically forbidden) region) defined by

    \displaystyle E - V_0 + Fx_0 = 0
    \displaystyle x_0 = \frac{V_0 - E}{F}

    The equation then becomes

    \displaystyle \psi''(x) + \frac{2mF}{\hbar^2}(x - x_0)\psi(x) = 0

    Introducing the dimensionless variable

    \displaystyle z = -\left(\frac{2mF}{\hbar^2}\right)^{1/3}(x - x_0)

    we finally obtain

    \displaystyle \frac{d^2\psi}{dz^2} - z\psi = 0

    Thus, the equation

    \displaystyle y'' + (1 - x)y = 0

    is, after a suitable rescaling and translation of the independent variable, an Airy equation. It describes, in particular, the behavior of a quantum wavefunction near a turning point, i.e. a position where

    \displaystyle E = V(x)

    In a region where E > V(x), the solution is oscillatory, whereas in a region where E < V(x), it exhibits exponential behavior. The Airy function provides the transition between these two regimes (see fig. below).

  • y” + (1-x) y= 0 (Part II)

    The first terms of the equation presented in the previous post are:

    \displaystyle a_0''(x) + a_0(x) - x a_0(x)\epsilon
    \displaystyle + a_1''(x)\epsilon + a_1(x)\epsilon - x a_1(x)\epsilon^2
    \displaystyle + a_2''(x)\epsilon^2 + a_2(x)\epsilon^2 - x a_2(x)\epsilon^3 + \cdots

    the differential equation can be written as:

    \displaystyle a_0''(x) + a_0(x) + (a_1''(x) + a_1(x) - x a_0(x))\epsilon + (a_2''(x) + a_2(x) - x a_1(x))\epsilon^2 + \cdots = 0
    \displaystyle a_0''(x) + a_0(x) = 0
    \displaystyle a_1''(x) + a_1(x) - x a_0(x) = 0
    \displaystyle a_2''(x) + a_2(x) - x a_1(x) = 0

    This implies

    \displaystyle a_0(x) = \cos(x)
    \displaystyle a_1(x) = \frac{1}{4}x^2 \sin(x) + \frac{1}{4} x \cos(x) - \frac{1}{4} \sin(x)

    We can check that a_1(x) is correct:

    \displaystyle a_1(x) = \frac{1}{4}x^2 \sin(x) + \frac{1}{4} x \cos(x) - \frac{1}{4} \sin(x)
    \displaystyle a_1'(x) = \frac{1}{4}x^2 \cos(x) + \frac{1}{2} x \sin(x) + \frac{1}{4} \cos(x) - \frac{1}{4} x \sin(x) - \frac{1}{4} \cos(x)
    \displaystyle a_1'(x) = \frac{1}{4}x^2 \cos(x) + \frac{1}{4} x \sin(x)
    \displaystyle a_1''(x) = \frac{1}{2}x \cos(x) - \frac{1}{4} x^2 \sin(x) + \frac{1}{4} \sin(x) + \frac{1}{4} x \cos(x)

    We can check that the condition:

    \displaystyle a_1''(x) + a_1(x) - x a_0(x) = 0

    is fulfilled. Therefore:

    \displaystyle y_{\epsilon}(x) = \cos(x) + \left(\frac{1}{4}x^2 \sin(x) + \frac{1}{4} x \cos(x) - \frac{1}{4} \sin(x)\right) \epsilon + \cdots

    Setting \epsilon = 1 gives the first-order perturbative approximation to the corresponding solution of the Airy equation.

    \displaystyle y(x) = \cos(x) + \frac{1}{4}x^2 \sin(x) + \frac{1}{4} x \cos(x) - \frac{1}{4} \sin(x) + \cdots

  • y” + (1-x) y= 0 (Part I)

    In this post we will solve the following basic second-order equation with corresponding initial condition:

    \displaystyle y'' + (1-x)y = 0, \quad y(0) = 1, \quad y'(0) = 0

    We rewrite the equation (inserting \epsilon):

    \displaystyle y'' + (1-\epsilon x)y = 0

    We expand y_{\epsilon}(x) as:

    \displaystyle y_{\epsilon}(x) = \sum_{n=0}^{\infty} a_n(x) \epsilon^n

    The differential equation becomes:

    \displaystyle y_{\epsilon}(x)'' + (1-\epsilon x)y_{\epsilon}(x) = 0
    \displaystyle \sum_{n=0}^{\infty} a_n''(x) \epsilon^n + (1-\epsilon x) \sum_{n=0}^{\infty} a_n(x) \epsilon^n = 0
    \displaystyle \sum_{n=0}^{\infty} a_n''(x) \epsilon^n + \sum_{n=0}^{\infty} a_n(x) \epsilon^n - x \sum_{n=0}^{\infty} a_n(x) \epsilon^{n+1} = 0

  • A nonlinear ODE (Part II)

    We could have directly solved the non-linear differential equation of the previous post exactly. This is now what we’re going to do, and we’ll be able to observe the equivalence between the exact solution and the corresponding perturbative series.

    \displaystyle y' + y + y^2 = 0
    \displaystyle y' = -(y + y^2)
    \displaystyle \frac{1}{y + y^2}\,dy = -\,dx
    \displaystyle \int \left(\frac{1}{y+y^2}\right) dy = - \int dx
    \displaystyle \int \frac{1}{y(1 + y)}\,dy = - \int dx
    \displaystyle \int \left(\frac{1}{y} - \frac{1}{1 + y}\right) dy = - \int dx
    \displaystyle \ln(y) - \ln(1 + y) = - x + C
    \displaystyle \frac{y}{y+1} = k e^{-x}
    \displaystyle \frac{y}{y + 1} - \frac{y+1}{y+1} = k e^{-x} -1
    \displaystyle - \frac{1}{y + 1} = k e^{-x} -1
    \displaystyle -\frac{1}{k e^{-x}-1} = y + 1
    \displaystyle y = \frac{-k e^{-x}}{k e^{-x}-1}

    Using the initial condition y(0) = 1 we have k=\frac{1}{2}. In a similar way, we can show that the exact solution to the perturbed differential equation :

    \displaystyle y' + y + \epsilon y^2 = 0

    is:

    \displaystyle y_{\epsilon} = \frac{-1}{\epsilon - \epsilon e^{x} - e^{x}}

    Which has the following Taylor series (in \epsilon) :

    \displaystyle \frac{-1}{\epsilon - \epsilon e^{x} - e^{x}} = e^{-x} - (e^{x}-1) e^{-2x} \epsilon + (e^{x}-1)^2 e^{-3x} \epsilon^2 + ...

    Which corresponds to the perturbative series derived previously

  • A nonlinear ODE (Part I)

    In this post we would like to solve the following non-linear first-order equation with corresponding initial condition:

    \displaystyle y' + y + y^2 = 0, y(0) = 1

    We rewrite the equation (inserting \epsilon):

    \displaystyle y' + y + \epsilon y^2 = 0

    and y(x)

    \displaystyle y_{\epsilon}(x) = \sum_{n=0}^{\infty} a_n(x) \epsilon^n

    The differential equation becomes:

    \displaystyle y' + y + \epsilon y^2 = 0
    \displaystyle \sum_{n=0}^{\infty} a_{n}'(x)\epsilon^{n} + \sum_{n=0}^{\infty} a_{n}(x)\epsilon^{n} + \epsilon\biggl(\sum_{k=0}^{\infty}\sum_{l=0}^{k}a_{l}(x) a_{k-l}(x)\epsilon^{k}\biggr) = 0

    Selecting the n-th term we have:

    \displaystyle a_n'(x) + a_n(x) + \sum_{k=1}^{n} a_{k-1} a_{n-k} = 0

    for n=0 we obtain:

    \displaystyle a_{0}'(x) + a_{0}(x) = 0
    \displaystyle a_{0}'(x) = - a_{0}(x)
    \displaystyle a_{0}(x) = e^{-x}

    for n=1 we obtain:

    \displaystyle a_{1}'(x) + a_{1}(x) + a_{0}^2(x) = 0
    \displaystyle a_{1}'(x) + a_{1}(x) + e^{-2x} = 0
    \displaystyle a_{1}'(x) = - a_{1}(x) - e^{-2x}

    which implies:

    \displaystyle a_{1}(x) = -e^{-x} + e^{-2x}

    Finally we have:

    \displaystyle y_{\epsilon}(x) = e^{-x} + (-e^{-x} + e^{-2x})\epsilon + ...

    This equation is related to the logistic equation, a classical model used in population dynamics to describe growth and decline. Here, the negative terms indicate that y(x) decreases over time and tends toward zero. The quadratic term adds a nonlinear effect to this decay. Similar equations can also be found in chemical kinetics and in other models involving nonlinear growth or decay.