Perturbation Theory part II

Now that we understand the basics of perturbation theory we can address some of the shortcomings in our original formulation. For example as it currently stands our theory of perturbation is unable to be used for so called singular problems, ones where we cannot simply set $\epsilon = 0$. Problems like this usually have the small parameter $\epsilon$ on the leading term. We are also unable to deal with secular terms in our perturbation series, terms that grow without bound. The original example I presented in my previous post was an example of such a problem and we shall attempt to address how to deal with it here.

Singular Perturbations

We start by considering a singular perturbation problem, first we shall attempt to use regular perturbation theory as we did previous and see how and why it fails. Consider the quadratic equation $\epsilon x^2 - x + 1 = 0, \ \epsilon \ll 1$ using the usual quadratic formula the solution to this problem can be trivially and is \[ x = \dfrac{1 \pm \sqrt{1 - 4 \epsilon}}{2 \epsilon} \] Again as we have done previously, let's assume we couldn't have just used the quadratic formula and attempt to find a solution using perturbation theory. Following the algorithm.

1. First let's identify the unperturbed problem, just set $\epsilon = 0$ and see what we're left with \[ -x + 1 = 0 \] The solution to this linear equation is \[ x = 1 \] Things already seem to be going wrong, by setting $\epsilon = 0$ we have lost a solution to the problem. Let's continue and see what happens.
2. Let's assume we can write $x$ in a perturbation series as \[x = x_0 + \epsilon x_1 + \epsilon^2 x_2 + \cdots \] Next we substituting this into our original equation to give us \[ \epsilon (x_0 + \epsilon x_1 + \epsilon^2 x_2 + \cdots)^2 - (x_0 + \epsilon x_1 + \epsilon^2 x_2 + \cdots) + 1 = 0 \] So we now have the full problem written in terms of a perturbation series, so let's proceed to the final step.
3. Now we solve the perturbation series by collecting terms of equal order in $\epsilon$ and solving the resulting system of equations. So collecting terms of equal orders \[ \left\{ \begin{array}{l l} \mathcal{O}(1) &: - x_0 + 1 &=& 0 \\ \mathcal{O}(\epsilon) &: {x_0}^2 - x_1 &=& 0 \\ \mathcal{O}(\epsilon^2) &: 2 x_0 x_1 - x_2 &=& 0 \\ \vdots \end{array} \right.\] We solve the system of equations up to second order in $\epsilon$, $\mathcal{O}(\epsilon^2)$, and find that the solution to the system is given by \[ x_0 = 1 , \quad x_1 = 1, \quad x_2 = 2 \] So putting this all together we find \[ x = 1 + \epsilon + \epsilon^2 + \cdots \] This is a catastrophic failure, by applying regular perturbation to this problem we've only managed to find one solution, when there must exist two.

So what went wrong? By setting $\epsilon=0$ we have reduced the order of our equation from a second order polynomial to a first order and thus lost a solution. To get around this issue we need to find a way to apply perturbation without lowering the order of the equation we want to solve. Fortunately there is a simple way to do this, we introduce a scaling factor.

Let's introduce the scaling factor $y = \delta(\epsilon) x$. If we apply this to our previous problem we find \[ \epsilon \delta^2 y^2 - \delta y + 1 = 0 \] Now we need to work out the scaling factor $\delta$. There are two methods to do this. We can either vary the size of $\delta$ from small to large until we find the appropriate size of $\delta$ that balances the equation and gives us both solutions. This method is something of an art as spotting the right intervals to choose can be difficult. This works best with larger problems. However here we shall use the simpler method which is better suited to smaller problems. We compare pairwise all combinations of terms until we find a meaningful balance.

1. When the first and second term are in balance, $\delta = \epsilon^{-1}$. The first and second terms are then of order $\mathcal{O}(\epsilon^{-1})$ which dominates the third term, which is of order $\mathcal{O}(1)$. This result is sensible and leads to two solutions.
2. When the first and third order terms are in balance, $\delta = \epsilon^{-1/2}$. The first and third terms are then of order $\mathcal{O}(1)$ and are dominated by the second order term which is of order $\mathcal{O}(\epsilon^{-1/2})$. A single term, the second term, dominates; it is of higher order than the balancing terms so this does not result in a dominant balance. So we obtain no new solutions.
3. When the second and third order terms are in balance $\delta = 1$. The second and third terms are then of order $\mathcal{O}(1)$ and dominate the first term which is of order $\mathcal{O}(\epsilon)$. Of course here we have just recovered the original singular solution that we obtained by setting $\epsilon=0$.

Our analysis here leads us to conclude that we should rescale our equation by using the change of variables $x = y/\epsilon$. This is the result of the dominant balance where we find two solutions. Rescaling our original equation and simplifying we find \[ y^2 - y + \epsilon = 0 \] Through our rescaling we've changed our singular equation into one that can be readily tackled using regular perturbation theory. For brevity I shall no go through solving this regular perturbation problem, following the method from my previous post we find that the solutions to this problem, to second order, are given by \[ \left\{ \begin{array}{l l} y_+ = 1 - \epsilon - \epsilon^2 \\ y_- = \epsilon + \epsilon^2 \end{array} \right.\] and if we transform this back into our original variable, $x$ we find \[ \left\{ \begin{array}{l l} x_+ = 1/\epsilon - 1 - \epsilon \\ x_- = 1 + \epsilon \end{array} \right.\] These are the perturbation series solutions of our original equation to second order in $\epsilon$, $mathcal{O}(\epsilon^2$. Of course now that we have two solutions, we need to check them against the actual answer to make sure our method gives us the correct solution. The analytic solution is \[ x = \dfrac{1 \pm \sqrt{1 - 4 \epsilon}}{2 \epsilon} \] Writing out the square root term as a power series gives us \[ \dfrac{1 \pm (1 - 2\epsilon - 2\epsilon^2 + \cdots}{2\epsilon} \] which is the same as the perturbation series solution. This is an excellent result as through a simple rescaling and balancing act we have found a way to solve singular problems and vastly expand the repertoire of problems we can solve.

Let's analyse more thoroughly why regular perturbation failed. Examining the solutions we can see that the standard ansatz of regular perturbation theory failed to capture the $1/ \epsilon$ term that appears in the solution. If instead we had chosen an ansatz of the form $x = \epsilon^{-1}x_{-1} + x_0 + \epsilon x_1 + \cdots$ we would have found both solutions, however without considerable insight or an exact solution it is very difficult to predict an expansion that will cover all of the solutions; so the more general method is to rescale and then apply regular perturbation theory. An interpretation of the root that we initially failed to capture is that as we let $\epsilon \to 0$ the root diverges and $x \to \infty$. Be careful as this interpretation is only true for algebraic equations.

Multiple scale analysis

As we are addressing singular perturbation theory now is an appropriate time to talk about dealing with secular terms. A secular term is one that grow without bound. These often appear when applying perturbation theory to differential equations with oscillatory solutions. This technique is very important for extending the range over which solutions are good, in fact in some cases it lets us lift the problem that perturbation series are only valid up until $x \sim 1/\epsilon$. A prime example of an equation we need this sort of method to fix is the non linear duffing equation. This is the differential equation I originally presented at the beginning of my previous post, and here we will develop the tools to solve it.

This non linear duffing equation with appropriate boundary values is \[ \left\{ \begin{array}{l l} y''(t) + y(t) + \epsilon y^3(t) = 0\\ y(0) = a \\ y'(0) = 0 \\ \end{array} \right.\] We initially approach this problem using regular perturbation theory, following the usual method

1. Let's set $\epsilon = 0$ and Identify the unperturbed problem \[ \left\{ \begin{array}{l l} y'' + y = 0\\ y(0) = a \\ y'(0) = 0 \\ \end{array} \right.\] This is standard second order ODE with constant coefficients which can be solved easily via the usual methods, it has the solution \[ y(t) = a \cos t \] Now we move on to identify the perturbed problem.
2. We assume that we can write $y(t)$ as a perturbation series \[ y(t) = y_0 + \epsilon y_1 + \epsilon^2 y_2 + \cdots \] We now substitute this into our original equation, including the boundary conditions and find \[ \left\{ \begin{array}{l l} (y_0 + \epsilon y_1 + \epsilon^2 y_2 + \cdots)'' + (y_0 + \epsilon y_1 + \epsilon^2 y_2 + \cdots) + \epsilon (y_0 + \epsilon y_1 + \epsilon^2 y_2 + \cdots)^3 = 0 \\ y_0(0) + \epsilon y_1(0) + \epsilon^2 y_2(0) + \cdots = a \\ {y_0}'(0) + \epsilon {y_1}'(0) + \epsilon^2 {y_2}'(0) + \cdots = 0 \end{array} \right.\] Remember that only the zeroth order boundary conditions are inhomogeneous in this problem. Now we've formulated the problem in terms of a perturbation series let's proceed to the final step.
3. Now we collect all the terms of equal order in $\epsilon$ and form a system of equations that we need to solve. The system is \[ \left\{ \begin{array}{l l} \mathcal{O}(1) &: {y_0}'' + {y_0} &=& 0 \\ \mathcal{O}(\epsilon) &: {y_1}'' + {y_1} + {y_0}^3 &=& 0 \\ \mathcal{O}(\epsilon^2) &: {y_2}'' + {y_2} + 3{y_0}{y_1}^2 &=& 0 \\ \vdots \end{array} \right.\] and has the corresponding system of boundary conditions \[ \left\{ \begin{array}{l l} \mathcal{O}(1) &: y_0(0) = a, \quad {y_0}'(0) = 0 \\ \mathcal{O}(\epsilon) &: y_1(0) = 0, \quad {y_1}'(0) = 0 \\ \mathcal{O}(\epsilon^2) &: y_2(0) = 0, \quad {y_2}'(0) = 0 \\ \vdots \end{array} \right.\] This system is very tedious to solve by hand, it can be done very easily using the method of undetermined coefficients and the usual techniques for solving homogeneous linear constant coefficient ODEs. I would recommend using a computer algebra system, such as Mathematica, to solve this system; in particular if you want to consider solving for higher order terms then a CAS really is a must. Regardless the solution to first order is \[ y(t) = a \cos t + \dfrac{\epsilon a^3}{32}(\cos 3t - \cos t - 12 t \cos t) + \mathcal{O}(\epsilon^2) \] We could do all of this before, but here we can identify the secular term, which is $t \cos t$; this term will grow without limit and causes the perturbation series to diverge from the actual solution as $x \sim 1/\epsilon$.

Now we introduce the idea of multiple scale analysis; here we assume the existence of many time scales in our problem, $( t, \tau, \tau_1, \cdots )$. We treat each of these time scales as its own independent variable. This is the general procedure, however for simplicity here let's consider only two timescales with the two variables $t$ and $\tau = \epsilon t$ and seek a perturbative solution of the form \[ y(t) = Y_0(t,\tau) + \epsilon Y_1(t,\tau) + \epsilon^2 Y_2(t,\tau) + \cdots \] Substituting this new two scale expansion into our original differential equation we find \[ \left\{ \begin{array}{l l} \mathcal{O}(1) &: \dfrac{\partial^2 Y_0}{\partial t^2} + Y_0 &=& 0 \\ \mathcal{O}(\epsilon) &: \dfrac{\partial^2 Y_1}{\partial t^2} + Y_1 + {Y_0}^3 + 2 \dfrac{\partial^2 Y_0}{\partial t \partial \tau} &=& 0 \\ \vdots \end{array} \right.\] Solving this system of equations we find, starting at zeroth order $\mathcal{O}(1)$ \[ Y_0 = A(\tau) e^{i t} + B(\tau) e^{-i t} \] Where $A(\tau)$ and $B(\tau)$ are functions we have yet to determine. We determine the two unknown functions by imposing the condition that no secular terms should appear in $Y_1$. When we solved this problem using regular perturbation there were of course secular terms in $Y_1$. Because of the functional freedom we have in $A$ and $B$ thanks to our new timescale we are free to do this. Before we move we should more thoroughly examine our solution, by invoking euler's formula we can actually reduce the freedom of the functions to \[ Y_0 = A(\tau) e^{i t} + A^*(\tau) e^{-i t} \] This lowers the number of unknowns we need to solve for which will simplify our calculations. Next we substitute this resulting into the first order equation, $\mathcal{O}(\epsilon)$ to find \[ \dfrac{\partial^2 Y_1}{\partial t^2} + Y_1 = - A^3 e^{3 i t} - (A^*)^3 e^{-3 i t} - \Big( 2i \dfrac{\partial A}{\partial \tau} + 3 A^2 A^* \Big) e^{i t} + \Big( 2i \dfrac{\partial A^*}{\partial \tau} - 3 A (A^*)^2 \Big) e^{-i t} = 0 \] We examine the above equation, the homogeneous equation has solutions that go like $Y_1 \sim e^{\pm i t}$, so we expect secular terms to arise in the solution due to the presence of inhomogeneous terms that go like $e^{\pm i t}$. However as we have the function freedom to choose $A$, we want to solve the follow system to remove the secular terms \[ \left\{ \begin{array}{l l} 2i \dfrac{\partial A}{\partial \tau} + 3 A^2 A^* &=& 0 \\ 2i \dfrac{\partial A^*}{\partial \tau} - 3 A (A^*)^2 &=& 0 \end{array} \right.\] However these equations are simply the complex conjugate of one another, so we need only find a solution to one of them. So solve this equations we express the complex values coefficient function $A$ as $A(\tau) = r(\tau) e^{i \theta(\tau)} $. Where $r$ and $\theta$ are real valued functions of $\tau$. Using the chain rule and the product rule we find \[ \Big( 2i \dfrac{dr}{d\tau} - 2 r(\tau) \dfrac{d \theta}{d \tau} + 3r^3(\tau) \Big) e^{i \theta(\tau)} = 0 \] Equating the real and imaginary parts we find \[ \left\{ \begin{array}{l l} \dfrac{d r}{d \tau} &= 0 \\ \dfrac{d\theta}{d\tau} &= \dfrac{3}{2} r^2(\tau) \end{array} \right.\] The solution to these can be found simply, $r(\tau) = r_0$ and $\theta(\tau) = 3/2 \ {r_0}^2 \tau + \theta_0 $. This gives us the general solution \[ A(\tau) = A(0) e^{i\theta_0 + i3/2 {r_0}^2 \tau} \] substituting this solution into our equation for $Y_0$ we get \[ Y_0 = 2 r_0 \cos \Big( \theta_0 + \dfrac{3}{2} {r_0}^2 \tau + t \Big) \] and now we just need to match the BC to complete our solution for the timescale $\tau = \epsilon t$. The BCs to match are $y(0) = a, \quad y'(0) = 0$, from which we find that \[ r_0 = a/2, \quad \theta_0 = 0 \] and we arrive at our two scale approximation to the non linear duffing equation \[ y(t) = a \cos \Big( \Big( 1 + \epsilon \dfrac{3a}{8} \Big) t \Big) + \mathcal{O}(\epsilon) \] We can repeat this process to determine the higher order terms by eliminating all the secular terms that appear via this method. The result of this procedure is that our perturbation series is now valid over a much larger range, as $\tau = \epsilon t$ our expansion remains good as long as $t \ll 1/\epsilon^2 $.

Final comment

When we've talked about perturbation theory it has always been for equations were there is some term that is small and we can set to zero. However perturbation theory also work when we have a very large term which dominates the equation. I mention this here as now we've talked about rescaling we can see why this is. Consider the quadratic equation we worked with earlier $ x^2 - x + \alpha = 0$, except this time $\alpha$ is a very large parameter $\alpha \gg 0$. We can transform this equation by making the variable change $x = \alpha y$ to get \[ y^2 - \dfrac{y}{\alpha} + \dfrac{1}{\alpha} = 0 \] If we now define $\epsilon = 1/\alpha$ we now have a regular perturbation theory problem \[ y^2 - \epsilon y + \epsilon = 0 \] Through a simple rescaling we can extend perturbation theory to deal with the very small and the very large.