Approximaths

Approximaths

  • 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.

  • Building the Sine Function Term by Term

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

    \displaystyle y'' = -y, \qquad y(0) = 0, \qquad y'(0)=1

    This equation can be solved very easily exactly. The answer is y(x) = sin(x). Now, let’s pretend that we don’t know how to solve this equation exactly. We introduce a parameter \epsilon and consider the family of equations

    \displaystyle y'' = -\epsilon y
    \displaystyle y'' +\epsilon y = 0

    and the function y(x) in the form:

    \displaystyle y_{\epsilon}(x) = a_0(x) + a_1(x)\epsilon + a_2(x)\epsilon^2 + a_3(x)\epsilon^3 + ...
    \displaystyle y_{\epsilon}''(x) = a_0''(x) + a_1''(x)\epsilon + a_2''(x)\epsilon^2 + a_3''(x)\epsilon^3 + ...

    The differential equation becomes (up to three terms):

    \displaystyle y'' +\epsilon y = 0
    \displaystyle a_0''(x) + a_1''(x)\epsilon + a_2''(x)\epsilon^2 + a_3''(x)\epsilon^3 + \epsilon (a_0(x) + a_1(x)\epsilon + a_2(x)\epsilon^2 + a_3(x)\epsilon^3) = 0
    \displaystyle a_0''(x) + (a_1''(x) + a_0(x))\epsilon + (a_2''(x) + a_1(x))\epsilon^2 + (a_3''(x) + a_2(x)) \epsilon^3 = 0

    This implies that:

    \displaystyle a_0''(x) = 0 \implies a_0(x) = x
    \displaystyle a_1''(x) + x = 0 \implies a_1(x) = -\frac{1}{3!}x^3
    \displaystyle a_2''(x) -\frac{1}{3!}x^3 = 0 \implies a_2(x) = \frac{1}{5!}x^5
    \displaystyle a_3''(x) +\frac{1}{5!}x^5 = 0 \implies a_3(x) = -\frac{1}{7!}x^7

    and

    \displaystyle y_{\epsilon}(x) = x -\frac{1}{3!}x^3 \epsilon + \frac{1}{5!}x^5 \epsilon^2 - \frac{1}{7!}x^7 \epsilon^3

    Setting \epsilon =1 to recover the original problem leads to:

    \displaystyle y_{\epsilon}(x) = x -\frac{1}{3!}x^3 + \frac{1}{5!}x^5 - \frac{1}{7!}x^7

    which corresponds to the first terms of the Taylor series of sin(x) (the exact solution).

  • Perturbation Theory in Action

    We would like to use perturbation methods similar to those presented for solving polynomials in the previous posts, but this time to solve differential equations. We will solve the following basic first-order equation with the corresponding initial condition:

    \displaystyle y' = y,\; y(0)=1

    This equation can be solved exactly very easily. The answer is y(x)=e^x. Now, let us pretend that we don’t know how to solve this equation exactly and use perturbation theory techniques. We write the function y(x) in the form (in a similar way to how we calculated the roots of polynomials):

    \displaystyle y_{\epsilon}(x)=\sum_{n=0}^{\infty}a_n(x)\epsilon^n=a_0(x)+a_1(x)\epsilon+a_2(x)\epsilon^2+\ldots

    According to the definition above, we have:

    \displaystyle y_{\epsilon}(x)=a_0(x)+a_1(x)\epsilon+a_2(x)\epsilon^2+a_3(x)\epsilon^3+\ldots
    \displaystyle y'_{\epsilon}(x)=a_0'(x)+a_1'(x)\epsilon+a_2'(x)\epsilon^2+a_3'(x)\epsilon^3+\ldots

    Now we can perturb the differential equation above:

    \displaystyle y'=\epsilon y

    and write (up to terms of degree 3):

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

    Therefore, we have to solve the following system of equations:

    \displaystyle a_0'(x)=0
    \displaystyle a_1'(x)-a_0(x)=0
    \displaystyle a_2'(x)-a_1(x)=0
    \displaystyle a_3'(x)-a_2(x)=0

    Applying the initial conditions a_0(0)=1 and a_n(0)=0 for n\ge1, we obtain:

    \displaystyle a_1'(x)-1=0\implies a_1=x
    \displaystyle a_2'(x)-a_1(x)=0\implies a_2=\frac{1}{2}x^2
    \displaystyle a_3'(x)-a_2(x)=0\implies a_3=\frac{1}{3!}x^3

    Therefore:

    \displaystyle y_{\epsilon}(x)=1+x\epsilon+\frac{1}{2}x^2\epsilon^2+\frac{1}{3!}x^3\epsilon^3

    and y(x) is approximated as follows (setting \epsilon=1):

    \displaystyle y_{\epsilon}(x)=1+x+\frac{1}{2}x^2+\frac{1}{3!}x^3

    which corresponds to the first terms of the Taylor series of e^x (the exact solution).

  • The Newton Polygon

    The construction based on the affine exponent functions (see this post)

    \displaystyle L_i(p)=m_ip+\alpha_i

    has a simple geometric counterpart known as the Newton polygon. Instead of representing each monomial

    \displaystyle \epsilon^{\alpha_i}x^{m_i}

    by the affine function L_i(p), we associate it with the point

    \displaystyle (m_i,\alpha_i)

    in the (m,\alpha) plane. The Newton polygon is defined as the lower convex hull of these points. Its edges contain exactly the information needed to determine the admissible scalings. Indeed, if an edge joins the points

    \displaystyle (m_i,\alpha_i) \qquad\text{and}\qquad (m_j,\alpha_j)

    its slope is

    \displaystyle s=\frac{\alpha_j-\alpha_i}{m_j-m_i}

    and the corresponding scaling exponent is

    \displaystyle p=-s=-\frac{\alpha_j-\alpha_i}{m_j-m_i}

    Thus, each edge of the Newton polygon corresponds to a dominant balance between monomials. The successive edges recover exactly the same admissible scalings as the corners of the lower envelope of the affine functions L_i(p). The Newton polygon therefore provides a dual geometric interpretation of the regularization process. As an illustration, consider the perturbed polynomial (see this post)

    \displaystyle P_{\epsilon}(x)=\epsilon^3x^3+x^2+\epsilon x

    The corresponding points in the (m,\alpha) plane are

    \displaystyle (3,3),\qquad (2,0),\qquad (1,1)

    The Newton polygon is the lower convex hull (see figure below)

    \displaystyle (1,1)\longrightarrow(2,0)\longrightarrow(3,3)

    The first edge has slope

    \displaystyle s=\frac{0-1}{2-1}=-1

    and therefore gives the scaling

    \displaystyle p=1

    The second edge has slope

    \displaystyle s=\frac{3-0}{3-2}=3

    and gives the scaling

    \displaystyle p=-3

    These two values are exactly the admissible scalings obtained from the intersections of the affine exponent functions

    \displaystyle 3p+3,\qquad 2p,\qquad p+1

    The left panel of the figure below shows the affine exponent functions and their lower envelope, while the right panel shows the corresponding Newton polygon. The two pictures are dual representations of the same dominant balance structure.