Approximaths

Approximaths

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

  • Regularize V

    The previous post illustrates the general strategy: the admissible scalings are obtained from the intersections of the affine exponent functions. In this post, we explain why this construction works. Consider a perturbed polynomial of the form

    \displaystyle P_{\epsilon}(x)=\sum_{i=1}^{n}\epsilon^{\alpha_i}x^{m_i}

    After the scaling

    \displaystyle x=\epsilon^{p}w

    each monomial becomes

    \displaystyle \epsilon^{\alpha_i}(\epsilon^p w)^{m_i} = \epsilon^{\alpha_i+m_ip}w^{m_i}

    Thus every monomial defines an affine function

    \displaystyle L_i(p)=m_ip+\alpha_i

    The exponents of \epsilon are therefore completely described by the family of affine functions

    \displaystyle L_1(p),L_2(p),\ldots,L_n(p)

    The key observation is that, as \epsilon\to0, the smallest exponent dominates all the others. Indeed, if

    \displaystyle \alpha<\beta

    then

    \displaystyle \epsilon^\alpha\gg\epsilon^\beta

    Consequently, for a fixed value of p, only the monomials whose exponent is minimal contribute to the dominant part of the polynomial.

    Now suppose that one affine function lies strictly below all the others. Then a single monomial dominates, and the limit polynomial consists of only one term.

    The interesting situation occurs when two affine functions intersect and their common value is the minimum among all the exponents. At such a point,

    \displaystyle L_i(p)=L_j(p)=\min_k L_k(p)

    so two monomials have exactly the same asymptotic order. After normalization, both survive in the limit polynomial, producing a non-trivial dominant part.

    Not every intersection is relevant. Two affine functions may intersect while another affine function remains strictly below them. In that case, the common exponent is not minimal, and the corresponding scaling does not contribute to the dominant polynomial.

    Therefore, the admissible scalings are obtained precisely from the intersection points where the common exponent is minimal.

    Geometrically, this means that the admissible scalings correspond to the corners of the lower envelope of the affine functions (see previous post). This simple geometric picture explains why the graph of the exponents completely determines the regularization process.

  • Regularize IV (scaling)

    Let us consider the following perturbed polynomial:

    \displaystyle P(x)=x^3+x^2+x
    \displaystyle P_{\epsilon}(x)=\epsilon^3x^3+x^2+\epsilon x
    \displaystyle P_{\epsilon}(\epsilon^p w)=\epsilon^3(\epsilon^p w)^3+(\epsilon^p w)^2+\epsilon(\epsilon^p w)
    \displaystyle =\epsilon^{3p+3}w^3+\epsilon^{2p}w^2+\epsilon^{p+1}w

    In the figure below, we plot the graphs of the affine functions 3p+3, 2p, and p+1. Their intersection points determine the values of p for which two exponents coincide. Among these, we retain only the intersections where the common exponent is minimal. In this example, we obtain

    \displaystyle (p_j,\nu_j)=(-3,-6)

    We now compute T_{\epsilon}^{(j)}(w) defined as follows:

    \displaystyle T_{\epsilon}^{(j)}(w):=\epsilon^{-\nu_j}P_{\epsilon}(\epsilon^{p_j}w)
    \displaystyle T_{\epsilon}^{(j)}(w)=\epsilon^{6}P_{\epsilon}(\epsilon^{-3}w)
    \displaystyle =\epsilon^{6}\bigl(\epsilon^{3(-3)+3}w^3+\epsilon^{2(-3)}w^2+\epsilon^{-3+1}w\bigr)
    \displaystyle =\epsilon^{6}\bigl(\epsilon^{-6}w^3+\epsilon^{-6}w^2+\epsilon^{-2}w\bigr)
    \displaystyle =\underbrace{w^3+w^2}_{T_{0}^{(j)}(w)}+\underbrace{\epsilon^{4}w}_{E_{\epsilon}^{(j)}(w)}

    Observe that E_{0}^{(j)}(w)=0. Therefore, we may write

    \displaystyle T_{\epsilon}^{(j)}(w)=T_{0}^{(j)}(w)+E_{\epsilon}^{(j)}(w)

    The polynomial T_{0}^{(j)}(w) is the dominant part of T_{\epsilon}^{(j)}(w) as \epsilon\to0, since it is independent of \epsilon. In this example, the dominant polynomial is non-degenerate. Therefore, the scaling x=\epsilon^{p_j}w regularizes the perturbed polynomial P_{\epsilon}(x).

    Graphs of the affine functions 3p+3, 2p, and p+1. The point (-3,-6) is the intersection of the lines 3p+3 and 2p. At p=-3, the minimum exponent is attained simultaneously by two monomials.

  • Regularization III

    In this post, we justify the regularization procedure introduced in previous posts for singular polynomial problems. Consider:

    \displaystyle P(x)=a_0+a_1x+\dots+a_nx^n
    \displaystyle P(x)=1+A_1x+\dots+A_nx^n
    \displaystyle P_{\epsilon}(x)=f_0(\epsilon)+f_1(\epsilon)\epsilon^{\alpha_1}A_1x+\dots+f_n(\epsilon)\epsilon^{\alpha_n}A_nx^n
    \displaystyle P_{\epsilon}(x)=(1+b_0\epsilon+c_0\epsilon^2+\dots)+(1+b_1\epsilon+c_1\epsilon^2+\dots)\epsilon^{\alpha_1}A_1x+\dots+(1+b_n\epsilon+c_n\epsilon^2+\dots)\epsilon^{\alpha_n}A_nx^n

    where f_0(\epsilon),f_1(\epsilon),\dots,f_n(\epsilon) admit power-series expansions in \epsilon, \alpha_i‘s are rational numbers and b_i‘s, c_i‘s are constants.

    The last expression is therefore a generalized version of a perturbed polynomial

    . We can rearrange this expression as follows:

    \displaystyle P_{\epsilon}(x)=(1+b_0\epsilon+c_0\epsilon^2+\dots)
    \displaystyle +(1+b_1\epsilon+c_1\epsilon^2+\dots)\epsilon^{\alpha_1}A_1x
    \displaystyle +\dots
    \displaystyle +(1+b_n\epsilon+c_n\epsilon^2+\dots)\epsilon^{\alpha_n}A_nx^n

    Collecting terms according to their powers of \epsilon, we obtain:

    \displaystyle P_{\epsilon}(x)=1+\epsilon^{\alpha_1}A_1x+\dots+\epsilon^{\alpha_n}A_nx^n
    \displaystyle +b_0\epsilon+b_1\epsilon^{\alpha_1+1}A_1x+\dots+b_n\epsilon^{\alpha_n+1}A_nx^n
    \displaystyle +\dots

    Theorem (Newton–Puiseux). Each root of the polynomial above is of the form:

    \displaystyle x(\epsilon)=\epsilon^pw(\epsilon),\quad w(0)\neq0

    where w(\epsilon) is a continuous function near zero. For convenience, we simply write x=\epsilon^pw. If the theorem holds, the polynomial above becomes:

    \displaystyle P_{\epsilon}(\epsilon^pw)=1+\epsilon^{\alpha_1}A_1\epsilon^pw+\dots+\epsilon^{\alpha_n}A_n(\epsilon^pw)^n
    \displaystyle +b_0\epsilon+b_1\epsilon^{\alpha_1+1}A_1\epsilon^pw+\dots+b_n\epsilon^{\alpha_n+1}A_n(\epsilon^pw)^n
    \displaystyle +\dots

    We define Q_{\epsilon}(w):=1+\epsilon^{\alpha_1+p}A_1w+\dots+\epsilon^{\alpha_n+np}A_nw^n and have:

    \displaystyle P_{\epsilon}(\epsilon^pw)=Q_{\epsilon}(w)+\epsilon(b_0+b_1\epsilon^{\alpha_1+p}A_1w+\dots+b_n\epsilon^{\alpha_n+np}A_nw^n)+\dots

    If \epsilon^pw is a root of P_{\epsilon}(x), we must satisfy:

    \displaystyle \lim_{\epsilon\to0}P_{\epsilon}(\epsilon^pw)=\lim_{\epsilon\to0}Q_{\epsilon}(w)=0

    Since w(\epsilon) is a continuous function near zero, we inspect the structure of Q_{\epsilon}(w(0)) to determine the leading-order behavior:

    \displaystyle Q_{\epsilon}(w(0))=1+\epsilon^{\alpha_1+p}A_1w(0)+\dots+\epsilon^{\alpha_n+np}A_nw(0)^n

    Considering the set of its exponents of \epsilon:

    \displaystyle E=\{0,\alpha_1+p,\dots,\alpha_k+kp,\dots,\alpha_n+np\}

    In order for the condition \lim_{\epsilon\to0}Q_{\epsilon}(w)=0 to hold, the minimal value of the set E must be exactly 0, and this minimum must be shared by at least two distinct exponents. These identical minimal exponents define the dominant terms of Q_{\epsilon}(w(0)) which can then compensate each other to balance the equation at leading order.

  • Regularize II

    \displaystyle

    The underlying idea of regularization is to recover the roots that are lost in a singular limit. When some coefficients of a polynomial are multiplied by powers of a small parameter \epsilon, the limiting polynomial obtained by setting \epsilon=0 may have fewer roots than the original polynomial. In this situation, some roots have escaped to a different scale and are therefore invisible in the limit.

    To recover these missing roots, we introduce a rescaling

    \displaystyle x:=\epsilon^p w

    where the exponent p is chosen so that at least two dominant terms of the polynomial balance each other. Geometrically, this choice is determined from the Newton polygon (or Newton diagram) by examining the exponents of \epsilon. After an additional normalization by a factor \epsilon^{-m_j}, one obtains a new limiting polynomial in the variable w. This reduced polynomial captures the asymptotic behavior of the roots on the corresponding scale and typically restores the correct number of roots.

    In this sense, the singularity is not an intrinsic property of the problem but rather the consequence of observing the problem at an inappropriate scale.

    As a simple example, consider

    \displaystyle P(x) = x^2 -x+1
    \displaystyle P_\epsilon(x) =\epsilon x^2-x+1

    Setting \epsilon=0 gives

    \displaystyle P_0(x)=-x+1

    which has only one root, although the original polynomial is quadratic and therefore has two roots. To recover the missing root, let

    \displaystyle x:=\epsilon^{-1}w

    Substituting into the polynomial gives

    \displaystyle P_\epsilon(\epsilon^{-1}w) = \epsilon(\epsilon^{-2}w^2)-\epsilon^{-1}w+1 = \epsilon^{-1}(w^2-w+\epsilon)

    Multiplying by \epsilon yields

    \displaystyle w^2-w+\epsilon

    Taking the limit \epsilon\to0 gives the regular polynomial

    \displaystyle w^2-w=w(w-1)

    which has two roots. One corresponds to the finite root seen in the original limit, while the second reveals a root located on the larger scale x=O(\epsilon^{-1}). The rescaling has therefore recovered the root that was hidden in the singular limit.