Case Study: Singular Perturbation Theory with SymPy

Mathematics | Programming | Python

Formal perturbation theory, and in particular singular perturbation theory, are standard topics in applied mathematics, but seem to be largely unknown outside that specific community. Which is a shame, because they are both useful, and intellectually fascinating.

Perturbation expansions, in particular for higher orders, have a reputation for being “cumbersome”, as the say: the algebra quickly becomes both tedious and error-prone. Which is true, but it turns out that the SymPy computer-algebra system can be put to good use in this context.

The Problem

Consider the equation (this is a strictly algebraic problem; extensions to differential equations are, of course, possible):

$$ \frac{1}{1 + x^2} = \epsilon x $$

Figure

What is desired is a simple, approximate expression for the solution $x^\star$ in terms of $\epsilon$. We will frequently, but not exclusively, assume $\epsilon$ to be small: $\epsilon \ll 1$.

Part 1: Theory

Ad-Hoc Approximations

Examining the graph, it is clear that when $\epsilon$ is small, $x^\star$ will be large. Hence we can approximate the left-hand side as follows:

$$ \frac{1}{1 + x^2} \approx \frac{1}{x^2} \qquad \qquad x \gg 1 $$

With this, equation becomes $1 = \epsilon x^3$ or:

$$ x^\star \sim \frac{1}{\epsilon^{1/3}} \qquad \qquad \epsilon \ll 1 $$

When $\epsilon$ is large, then the straight line in the figure is almost vertical, and hence $x^\star$ will be small. Hence we can expand on the left-hand side:

$$ \frac{1}{1 + x^2} \approx 1 \qquad \qquad x \ll 1 $$

and therefore:

$$ x^\star \sim \frac{1}{\epsilon} \qquad \qquad \epsilon \gg 1 $$

Exact Solution

The original problem is, in fact, a cubic equation. An equation of the form: $x^3 + px + q = 0$, with $(q/2)^2 + (p/3)^3 > 0$ has the explicit real solution:

$$ x = \sqrt[3]{-\frac{q}{2} + \sqrt{\left(\frac{q}{2}\right)^2 + \left(\frac{p}{3}\right)^3}} + \sqrt[3]{-\frac{q}{2} - \sqrt{\left(\frac{q}{2}\right)^2 + \left(\frac{p}{3}\right)^3}} $$

This formula, due to Tartaglia and Cardano, is (as Feynman once put it) “of very little use in itself”, but will nevertheless be interesting for comparison with the approximate expressions later on.

Applied to our problem, which is $x^3 + x - 1/\epsilon$, the condition $(q/2)^2 + (p/3)^3 > 0$ is always true, and the formula becomes:

$$ x^\star = \sqrt[3]{\frac{1}{2 \epsilon} + \sqrt{\left(\frac{1}{2 \epsilon}\right)^2 + \left(\frac{1}{3}\right)^3}} + \sqrt[3]{\frac{1}{2 \epsilon} - \sqrt{\left(\frac{1}{2 \epsilon}\right)^2 + \left(\frac{1}{3}\right)^3}} $$

Below is a figure of the exact solution, together with the two ad-hoc approximations. They are reasonable in their respective limits ($\epsilon \to 0$ or $\epsilon \to \infty$), but there does not seem to be a way to improve either of them systematically. Doing so will be the main purpose of the rest of this write-up.

Figure

Aside: Formal Perturbation Expansions

In formal perturbation theory, we are seeking an approximate solution to an equation, in terms of an expansion in some “small parameter” $\epsilon$:

$$ x \approx x_0 + \epsilon^{\alpha_1} x_1 + \epsilon^{\alpha_2} x_2 + \dots \qquad \qquad \alpha_1 < \alpha_2 < \dots $$

The powers $\alpha_1, \alpha_2, \dots$ are not restricted to be integers. Finding values for them will be a central part of the problem.

In order to make the problem well-defined, we do require that the exponents are ordered, with each successive power being strictly greater than any of the preceding ones.

In our case, we are seeking an approximate solution to:

$$ \epsilon x^3 + \epsilon x - 1 = 0 $$

in the limits of both small and large $\epsilon$.

Singular Perturbation Expansion for Small $\epsilon$

The case of small $\epsilon$ is the more interesting of the two, because in the limit of $\epsilon \to 0$, the solution vanishes — or rather, it escapes to infinity. An indication that there is something amiss is the fact that $\epsilon$ multiplies the highest power in the equation: when $\epsilon$ becomes zero, the nature of the equation changes.

This marks the problem as a singular perturbation problem, and the standard remedy is to re-write the problem in terms of a scaled variable:

$$ u = \frac{x}{\epsilon^\eta} \quad \Leftrightarrow \quad x = \epsilon^\eta u $$

with the exponent $\eta$ to be determined. Substituting $u$ into the original equation leads to:

$$ \epsilon^{3\eta + 1} u^3 + \epsilon^{\eta + 1} u - 1 = 0 $$

In this equation, $\epsilon$ is a free parameter and can, in principle, take on any value. On the other hand, the right side of the equation does not depend on $\epsilon$. For this equation to be valid, the terms in $\epsilon$ on the left side must, somehow, cancel.

The basic idea of formal perturbation theory is to adjust the exponents and coefficients in the perturbation expansion in such a way as to facilitate these cancellations, in strictly increasing order of the expansion parameter.

Because this is a singular perturbation problem, we first need to find the scaling exponent $\eta$, by balancing the term containing the highest power of the variable (that is, $u^3$) with one of the other terms. We have thus two choices:

  • Balance $\epsilon^{3\eta + 1} u^3$ with $\epsilon^{\eta + 1} u$: equating exponents of $\epsilon$ leads to $\eta = 0$.
  • Balance $\epsilon^{3\eta + 1} u^3$ with $-1$: equating exponents leads to $\eta = -1/3$.

Because we want an expansion in strictly increasing order of $\epsilon$, we choose the smallest of the possible exponents, which is $\eta = -1/3$. With this choice, the equation becomes:

$$ u^3 + \epsilon^{2/3} u - 1 = 0 $$

In other words, introducing the scaled variable $u$, and determining the exponent $\eta$ has transformed the original singular perturbation problem in $x$ into a regular problem in $u$. That is the central premise of singular pertubation theory.

Now taking the limit $\epsilon \to 0$ leaves us with $u^3 = 1$, so that we can conclude that the lowest order coefficient is $u_0 = 1$.

$\eta = -\frac{1}{3} \qquad \text{and} \qquad u_0 = 1$

First Order

To find the first order correction, we make the ansatz (remembering that we have already determined that $u_0 = 1$):

$$ u = u_0 + \epsilon^{\alpha_1} u_1 = 1 + \epsilon^{\alpha_1} u_1 $$

Using this expression in the scaled equation $u^3 + \epsilon^{2/3} u - 1 = 0$ leads to:

$$ (1 + \epsilon^{\alpha_1} u_1)^3 + \epsilon^{2/3} (1 + \epsilon^{\alpha_1} u_1) - 1 = 0 $$

or, multiplying out and canceling terms independent of $u_1$ and $\epsilon$:

$$ 3 \epsilon^{\alpha_1} u_1 + 3 \epsilon^{2 \alpha_1} u_1^2 + \epsilon^{3 \alpha_1} u_1^3 + \epsilon^{2/3} + \epsilon^{\alpha_1 + 2/3} u_1 = 0 $$

Now, again, we face the balancing problem: the right-hand side does not depend on $\epsilon$, but the left-hand side does. We therefore need to determine $\alpha_1$ in such a way that the lowest order terms on the left side cancel.

Equating the various exponents leads (in principle) to ten different combinations:

  • $\alpha_1 = 2 \alpha_1 = 3 \alpha_1 \Rightarrow \alpha_1 = 0$
  • $\alpha_1 = 2/3$
  • $\alpha_1 = \alpha_1 + 2/3 \Rightarrow \alpha_1 = 0$
  • $2 \alpha_1 = 2/3 \Rightarrow \alpha_1 = 1/3$
  • $2 \alpha_1 = \alpha_1 + 2/3 \Rightarrow \alpha_1 = 2/3$
  • $3 \alpha_1 = 2/3 \Rightarrow \alpha_1 = 2/9$
  • $3 \alpha_1 = \alpha_1 + 2/3 \Rightarrow \alpha_1 = 1/3$
  • $2/3 = \alpha_1 + 2/3 \Rightarrow \alpha_1 = 0$

We rule out any solution of the form $\alpha_1 = 0$, because we already have found the zero-th order term: $u_0 = 1$ and we are therefore looking for the next higher order term. The solution $\alpha_1 = 1/3$ might seem like a promising candidate, but when we use it in the equation, we are left with a single term $3 \epsilon^{1/3} u_1$: no cancellation is taking place. The same holds true for $\alpha_1 = 2/9$. The consequence is that $\alpha_1$ must be $\alpha_1 = 2/3$. Using this value in the equation, the lowest order terms are:

$$ 3 \epsilon^{2/3} + \epsilon^{2/3} + \text{terms higher order in $\epsilon$} = 0 $$

Hence:

$\alpha_1 = \frac{2}{3} \qquad \text{and} \qquad u_1 = -\frac{1}{3}$

The expansion, to first order, is therefore:

$u = 1 - \frac{1}{3} \epsilon^{2/3}$

Second Order

To find the next order correction, we make the ansatz:

$$ u = 1 - \frac{1}{3} \epsilon^{2/3} + \epsilon^{\alpha_2} u_2 $$

and proceed as before. At this point, the algebra is starting to become tedious, and we may want to employ SymPy to help with it (see below). The result, after substituting the ansatz into the equation and canceling terms, is:

$$ -\frac{1}{3^3}\epsilon^2 + \epsilon^{3 \alpha_2} u_2^{3} + 3 \epsilon^{2 \alpha_2} u_2^{2} + 3 \epsilon^{\alpha_2} u_2 - \epsilon^{\alpha_2 + \frac{2}{3}} u_2 + \frac{1}{3} \epsilon^{\alpha_2 + \frac{4}{3}} u_2 - \epsilon^{2 \alpha_2 + \frac{2}{3}} u_2^{2} = 0 $$

At this point, we have again to adjust $\alpha_2$ to eliminate the lowest order terms on the left-hand side. The whole “balancing” procedure can seem somewhat arbitrary or mysterious, but the rules are not actually complicated. The idea is to find the smallest $\alpha$, which is strictly greater than the last order in the expansion so far, and which, when substituted into the equation, leaves at least two terms of lowest order, which can cancel each other. To put the last condition differently: find candidate values for $\alpha$ by equating exponents of different terms, then substitute every possible candidate value into the equation. In the resulting expression, examine the terms of lowest order in $\epsilon$: there should be enough to allow for cancellation; if there is only a single term, reject that value of $\alpha$ and try the next larger one.

If we apply the process to the expression above, we find:

$\alpha_2 = 2 \qquad \text{and} \qquad u_2 = \frac{1}{3^4}$

The expansion, up to and including second order, is therefore:

$u = 1 - \frac{1}{3} \epsilon^{2/3} + \frac{1}{3^4} \epsilon^2$

Higher Orders

The algebra becomes increasingly cumbersome, but the SymPy program discussed below can calculate higher-order terms without human intervention.

The 12-term expansion is:

$$ u = 1 - \frac{1}{3} \epsilon^{2/3} + \frac{1}{3^4} \epsilon^2 + \frac{1}{3^5} \epsilon^{8/3} - \frac{4}{3^8} \epsilon^4 - \frac{5}{3^9} \epsilon^{14/3} + \frac{77}{3^{13}} \epsilon^6 + \frac{104}{3^{14}} \epsilon^{20/3} - \frac{595}{3^{17}} \epsilon^8 - \frac{836}{3^{18}} \epsilon^{26/3} + \frac{5083}{3^{21}} \epsilon^{10} + \frac{7315}{3^{22}} \epsilon^{32/3} $$

Full Expansion

Finally, in order to get back to the quantities in the original problem, we have to replace the scaled variable $u$ with the original, un-scaled variable $x = \epsilon^{-1/3} u$. Substituting this in the full expansion yields the final result of this section:

$$ x = \frac{1}{\epsilon^{1/3}} - \frac{1}{3} \epsilon^{1/3} + \frac{1}{3^4} \epsilon^{5/3} + \frac{1}{3^5} \epsilon^{7/3} - \frac{4}{3^8} \epsilon^{11/3} - \frac{5}{3^9} \epsilon^{13/3} + \frac{77}{3^{13}} \epsilon^{17/3} + \frac{104}{3^{14}} \epsilon^{19/3} - \frac{595}{3^{17}} \epsilon^{23/3} - \frac{836}{3^{18}} \epsilon^{25/3} + \frac{5083}{3^{21}} \epsilon^{29/3} + \frac{7315}{3^{22}} \epsilon^{31/3} $$

Regular Perturbation Expansion for Large $\epsilon$

Things are a bit simpler in the opposite case of large $\epsilon$. If $\epsilon$ becomes large, then $\delta = 1/\epsilon$ becomes the new “small” perturbation parameter. The equation in this case is:

$$ x^3 + x - \delta = 0 $$

The problem is non-singular. As before, we begin with a simple ansatz:

$$ x = x_0 + \delta^\alpha_1 x_1 $$

Substituting this into the equation leads to:

$$ x_0^3 + 3 \delta^{\alpha_1} x_0^2 x_1 + 3 \delta^{2 \alpha_1} x_0 x_1^2 + \delta^{3 \alpha_1} x_1^3 + x_0 + \delta^{\alpha_1} x_1 - \delta = 0 $$

Taking the $\delta \to 0$ limit, we end up with the zero-order expression:

$$ x_0^3 + x_0 = 0 $$

Hence, $x_0 = 0$ and all terms containing $x_0$ vanish, so that we are left with:

$$ \delta^{3 \alpha_1} x_1^3 + \delta^{\alpha_1} x_1 - \delta = 0 $$

Again, we try to balance terms. The choice $\alpha_1 = 0$ is ruled out, because we already determined the zero-order term $x_0$. It may appear as if $\alpha_1 = 1/3$ might be an option (to balance $\delta^{3 \alpha_1} x_1^3$ with $\delta$), but this choice would lead to a single, unbalanced, lower-order term $\delta^{1/3} x_1$. Hence, the only viable choice is to balance $\delta^{\alpha_1} x_1$ with $- \delta$, leading to:

$\alpha_1 = 1 \qquad \text{and} \qquad x_1 = 1$

and:

$x = 0 + \delta$

Second Order

Because we already found the zero-th and first order terms, the next order is the second. We make the usual ansatz:

$$ x = \delta + \delta^\alpha_2 x_2 $$

which leads to:

$$ \delta^3 + 3 \delta^{2+\alpha_2} x_2 + 3 \delta^{1+2 \alpha_2} x_2^2 + \delta^{3 \alpha_2} x_2^3 + \delta^\alpha_2 x_2 = 0 $$

Trying to balance terms, we find that the exponent must be $\alpha_2 = 3$, which leads to:

$\alpha_2 = 3 \qquad \text{and} \qquad x_2 = -1$

and

$x = \delta - \delta^3$

Higher Orders

At this point, it makes sense to employ SymPy to take care of the algebra. The full expansion, to 12 non-zero terms (remember that the zero order term vanishes for this problem), is:

$$ x = \delta - \delta^{3} + 3 \delta^{5} - 12 \delta^{7} + 55 \delta^{9} - 273 \delta^{11} + 1428 \delta^{13} - 7752 \delta^{15} + 43263 \delta^{17} - 246675 \delta^{19} + 1430715 \delta^{21} - 8414640 \delta^{23} $$

Also keep in mind that this is an expansion in $\delta = 1/\epsilon$, which is valid for $\epsilon \gg 1$.

Examining the Results

So far, this has just been a lot of algebra. To believe it, we need to see some pictures.

Figure Figure

Two things stand out from these images immediately:

  • In their range of validity, the approximations are quite good.
  • As is generally the case with asymptotic expansions, more terms mean a closer approximation in the respective limit, but do not necessarily imply a more useful approximation over a larger parameter range!

With the benefit of having the exact solution, we can identify the three-term expansions as the most suitable. Here they are, together with the exact solution.

Figure

And the final expansion takes the form:

$$ x = \begin{cases}\dfrac{1}{\epsilon^{1/3}} - \dfrac{1}{3} \epsilon^{1/3} + \dfrac{1}{3^4} \epsilon^{5/3} & \epsilon < 5 \\ \delta - \delta^{3} + 3 \delta^{5} & \delta = \dfrac{1}{\epsilon} > 3 \end{cases} $$

There is a rather large overlap region, where the two different expansions are almost the same. Unfortunately (and in contrast to problems arising out of differential equations), there is no free parameter in the problem that we could adjust to “match” the two solutions up to form a smooth, uniform approximation.

Concluding Thoughts - Balancing Revisited

The critical step in the entire calculation is the “balancing” of terms, order by order, which determines the exponents and coefficients in the perturbation expansion. It might be worth to re-state the rules that apply to this step:

Find the smallest exponent $\alpha$, which is strictly greater than the last order in the expansion so far, and which, when substituted into the equation, leaves at least two terms of lowest order, which can then be made to cancel each other, so that the equation is correct, up to that order. To put the last condition differently: find candidate values for $\alpha$ by equating the exponents for all pairs of terms, then substitute every possible candidate value for the exponent in the equation. In the resulting expression, there should be at least two terms of lowest order, which can be made to cancel each other by adjusting their coefficients. If there is only a single term, the value for the exponent has to be rejected, and the next larger candidate must be tried.

This process is clear enough, but it certainly leaves some questions. For example, it is not clear whether this process will always “work”, meaning that there will always be terms that can be made to cancel each other explicitly. The prescription also has to be augmented when the problem involves transcendental, not just algebraic, functions. (See the books below for more information.)

Resources

There are many books on perturbation theory, but at times it can be difficult to see how to apply their examples to one’s own problems.

Two books that I have found valuable, because they emphasize the systematic method of perturbation theory, rather than spectacular examples, are the two volumes by Mark H. Holmes:

  • Introduction to Perturbation Methods, 2nd ed, Springer (2015)
  • Introduction to the Foundations of Applied Mathematics, 2nd ed, Springer (2019)

The first emphasizes the mathematical methods of formal perturbation theory (Boundary Layers and Matched Asymptotic Expansions, Two-Timing and Multiple Scales, WKB, and the possibly less familiar method of “Homogenization”). The second discusses the application of perturbation methods in various domains, including Diffusion, Traffic Flow, and Continuum Mechanics.

Part 2: SymPy

The SymPy script below calculates both the singular and the regular perturbation expansion for the problem $\epsilon x ( x^2 + 1 ) - 1 = 0$, and prints the results to screen. It is fully automated, and does not require manual, interactive intervention.

WARNING: This script has been developed and validated only for the equation given in the preceding problem statement. It makes various assumptions about the form of intermediate expressions. These assumptions are likely to be incorrect for arbitrary other problems.

  1from sympy import *
  2
  3# Turn on extra output
  4diagnostics = False
  5warnings = False
  6
  7# Print only if global var "warnings" is True
  8def warn( *args ):
  9    if warnings:
 10        print( "Warning:", *args )
 11
 12# Substitute sub for var in eq, expand, then collect powers of h
 13def makeAnsatz( eq, var, sub, h ):
 14    return collect( powsimp( expand( eq.subs( var, sub ) ) ), h )
 15
 16# Given a term like a*h**n, extract the exponent n
 17# This uses the identity n = h*diff(a*h**n, h)/(a*h**n)
 18def getPower( term, h ):
 19    return cancel(h*diff( term, h )/term)
 20
 21# Given eq == sum_i u_i h**a_i, return an array of exponents of h: [ a_i ]
 22# Returns can't be sorted, since not necessarily numeric.
 23def getExponents( eq, h ):
 24    return [ getPower(t, h) for t in eq.args ]
 25
 26# Given a set of exponents, form all distinct pairs, equate them, and
 27# solve for var. Return a sorted list of solutions (ascending), all of
 28# which are strictly greater than the supplied cutoff (such as the
 29# highest order in the expansion so far).
 30# Issue a warning if the number of solutions does not equal 1.
 31def equateExponents( exps, var, cutoff ):
 32    rs = {}
 33    for i in range(0, len(exps)):
 34        for j in range( i+1, len(exps)):
 35            r = solve( exps[i] - exps[j], var )
 36
 37            if len(r) != 1:
 38                warn( exps[i], "==", exps[j], "->", r )
 39                continue
 40
 41            r = r[0]
 42            if r > cutoff:
 43                rs[r] = True
 44
 45    return sorted(list(rs))
 46
 47# Given an eq == sum_i u_i h**a_i, with all a_i numeric, return the
 48# numerically smallest exponent encountered.
 49def minExpo( eq, h ):
 50    tmp = []
 51    for t in eq.args:
 52        tmp.append( getPower( t, h ) )
 53    return sorted(tmp)[0]
 54
 55# Perform the balancing operation.
 56# Given an expression eq and a SORTED set of candidate exponent values exps
 57# (as returned by equateExponents()), attempt to substitute each candidate
 58# for the quantity var1 in eq. Find the minimum resulting exponent of var2
 59# in eq. If greater zero, calculate the limit of eq/var2**mn. If at least
 60# two terms are left in the limit, return those terms and the current
 61# minimum exponent.
 62# Prints diagnostic output if enabled.
 63# This routine does NOT stop as soon as it finds a solution, but continues
 64# searching, in order to print complete diagnostic output. But the values
 65# returned eventually are the first matching values.
 66def balanceTerms( eq, var1, exps, var2 ):
 67    for ex in exps:
 68        tmp = powsimp(expand(eq.subs( var1, ex )))
 69        mn = minExpo(tmp, var2)
 70
 71        if diagnostics:
 72            print( ex, ":", mn )
 73            pprint(tmp);
 74            print()
 75
 76        terms = limit( tmp/var2**mn, var2, 0 )
 77        if mn > 0 and len(terms.args) > 1:
 78            r1, r2 = mn, terms
 79
 80    return r1, r2
 81
 82# Given the terms remaining from the limit process in balanceTerms(),
 83# equate them to zero, and solve for var.
 84# Issue a warning if the number of solutions is not 1.
 85def findCoeff( terms, var ):
 86    coeff = solve( terms, var )
 87
 88    if len(coeff) != 1:
 89        warn( terms, "--", coeff )
 90        return None
 91
 92    return coeff[0]
 93
 94
 95# ==================================================
 96
 97# These are placeholders in eqs only. Use "h" for "epsilon"
 98x, u, v, h, a = symbols( "x u v h a", Real=True )
 99
100# Abbreviations for two commonly occurring rationals
101r13, r23 = Rational(1,3), Rational(2,3)
102
103# The equations defining the problem
104original = h*x**3 + h*x - 1
105singular = u**3 + u*h**r23 - 1
106regular = x**3 + x - h
107
108
109# Solve the singular problem
110print( "=== Solving the singular problem ===\n" )
111
112s, expo = 1 + v*h**a, -r13   # expo is needed as cutoff in equateExponents()
113for i in range(11):
114    tmp = makeAnsatz( singular, u, s, h )
115    exps = getExponents( tmp, h )
116    exps = equateExponents( exps, a, expo )
117    expo, terms = balanceTerms( tmp, a, exps, h )
118    coeff = findCoeff( terms, v )
119
120    s = s + coeff*h**expo
121
122    # This is the solution for the SCALED variable u!
123    pprint( s.subs(h**a, 0) )
124    print("\n", "="*30, "\n", sep="" )
125
126
127# Solve the regular problem
128print( "=== Solving the reular problem ===\n" )
129
130s, expo = 0 + v*h**a, 0   # expo is needed as cutoff in equateExponents()
131for i in range(12):
132    tmp = makeAnsatz( regular, x, s, h )
133    exps = getExponents( tmp, h )
134    exps = equateExponents( exps, a, expo )
135    expo, terms = balanceTerms( tmp, a, exps, h )
136    coeff = findCoeff( terms, v )
137
138    s = s + coeff*h**expo
139
140    # This is the solution for the ORIGINAL variable x
141    pprint( s.subs(h**a, 0) )
142    print("\n", "="*30, "\n", sep="" )