Past exam of the mathematics course of the University of Cambridge / 2014 / iii / Paper 74 (Perturbation and Stability Methods)

The small root and the root near $-2$ are regular perturbation roots: they stay finite as the parameter vanishes. Substitution of power series givesFor the third root, dominant balance for algebraic roots requires $\lambda =\mu /\epsilon$. The largest terms give ${\mu }^2(\mu +1)=0$, and the missing branch is $\mu =-1$. The sum of the three roots then yieldsThus the leading roots are $-\epsilon /2$, $-2$, and $-1/\epsilon$; only the last is a divergent perturbation root. The first regular branch happens to have a zero limit, so retaining its first nonzero term is essential on the long spatial scale.

First take $\epsilon \to 0^{+}$, the decaying half-line regime. If the three exact roots are $r_j$, the exact initial value problem solution isThis follows either by solving the three initial-value equations or from the Laplace transform $\hat{y}(s)=\epsilon /(\epsilon s^3+s^2+2s+\epsilon )$. The coefficient formula has sums $\sum A_j=\sum r_j{A}_j=0$ and $\sum r_j^2{A}_j=1$.

Retain the three leading roots but compute their coefficients without expanding their denominators. Set $a=\epsilon /2$, $b=2$, $d=1/\epsilon$. A convenient composite asymptotic expansion isIt satisfies all three initial conditions exactly and retains the fast transient, the ordinary decay, and the slow decay. Its first two coefficients are $\epsilon /2+O({\epsilon }^2)$ and $-\epsilon /2+O({\epsilon }^2)$; the fast coefficient is ${\epsilon }^2+O({\epsilon }^3)$. Consequently the simpler bulk expression is $\frac{\epsilon }{2}(e^{-\epsilon x/2}-e^{-2x})$, but that expression alone does not reproduce the initial derivative layer.

The absolute error estimate isTo see its uniformity, the slow-root error is $O({\epsilon }^2)$, while its magnitude is $O(\epsilon )$ and its coefficient is $O(\epsilon )$. The uniform error bound for nearby decaying exponentials therefore gives an $O({\epsilon }^2)$ contribution even for $x=O({\epsilon }^{-1})$. The middle-root error is $O(\epsilon )$ with coefficient $O(\epsilon )$, again giving $O({\epsilon }^2)$. The fast-root error is $O(1)$ with magnitude $O({\epsilon }^{-1})$ and coefficient $O({\epsilon }^2)$, giving only $O({\epsilon }^3)$. Coefficient errors are $O({\epsilon }^2)$ or smaller. This estimate concerns $y$ itself, rather than asserting the same uniform order for all derivatives or a relative error at its zero.

For the sketch, on $x=\epsilon X$ one obtainsThus $y\sim x^2/2$ at the origin and rises smoothly after the fast layer. For $x=O(1)$ it rises towards a plateau of height approximately $\epsilon /2$, then decays on the much longer scale $x=O({\epsilon }^{-1})$. The bulk maximum lies near $x=\log (4/\epsilon )/(2-\epsilon /2)$ and has height asymptotic to $\epsilon /2$.

The sign of the parameter matters on an infinite interval. If $\epsilon \to 0^{-}$ is allowed, the singular mode grows rather than decays. For each fixed $x>0$, its dominant contribution is $y(x)\sim {\epsilon }^2\exp [(-1/\epsilon +2)x]$. The extra $2$ in that exponent is needed for relative leading accuracy at fixed $x$. The positive-parameter uniform absolute bound and decaying sketch do not extend to that regime.

Apply the stationary phase method to the real part of the exponential integral for the Bessel function of the first kind. With $n$ fixed, use oscillatory phase $\sin \vartheta$ and oscillatory integral amplitude $e^{-in\vartheta }$. The unique stationary point is ${\vartheta }_0=\pi /2$, with second derivative $-1$. Its contribution isThe leading endpoint terms of the exponential integral are imaginary and do not contribute to its real part. Hence the large-argument asymptotic expansion of the Bessel function of the first kind isThis is an additive asymptotic formula: a relative ratio is inappropriate at zeros of the leading cosine. The envelope decreases like $x^{-1/2}$ while the oscillation frequency approaches one.

Initially suppose $0<\alpha <\pi /2$, the usual oscillatory range. The oscillatory phase is $n\Phi (\vartheta )$ with $\Phi =\sec \alpha \sin \vartheta -\vartheta$. Its stationary point is $\vartheta =\alpha$, since $\sec \alpha \cos \vartheta =1$. At that pointThe stationary phase method therefore gives the Debye Bessel asymptoticThe error is additive for fixed positive $\alpha$ bounded away from $\pi /2$; this formula is not uniform as $\alpha \to 0$, when the stationary point joins the endpoint and its curvature vanishes.

The printed condition $\alpha >0$ alone includes other trigonometric branches. All defined fixed real cases can be covered as follows. If $0<|\cos \alpha |<1$, put $\beta =\arccos |\cos \alpha |\in (0,\pi /2)$. Use the displayed formula with $\beta$ instead of $\alpha$, and multiply it by $1$ if $\sec \alpha >0$, or by $(-1{)}^n$ if $\sec \alpha <0$. This follows from the integer-order parity $J_n(-x)=(-1{)}^n{J}_n(x)$, which is also obtained from the defining integral by $\vartheta \mapsto \pi -\vartheta$. If $|\cos \alpha |=1$, the argument is $\pm n$ and the turning-point answer below applies, with the same parity factor. If $\cos \alpha =0$, the stated argument is undefined.

The Bessel turning-point asymptotic comes from a cubic stationary endpoint, not an ordinary quadratic stationary point. Near $\vartheta =0$,Thus the contributing width is $\vartheta =O(n^{-1/3})$. On writing $\vartheta =(6/n{)}^{1/3}s$, the leading integral isThe oscillatory integral is understood with a vanishing damping factor. Substitution $u=s^3$ and the Gamma function Fourier integral giveThereforeThe equality uses the Gamma reflection formula. Contributions away from the degenerate endpoint are smaller. The $n^{-1/3}$ scale explains why the preceding $n^{-1/2}$ formula cannot be extended directly to zero angle.

Introduce the slow time $T=\epsilon t$ and treat $t,T$ independently in the method of multiple scales. WriteAt the next order,For fixed $R$ define the period average byHere $f$ inside the averages is evaluated at the leading position and velocity $R\cos \varphi ,-R\sin \varphi$. Orthogonality to both fundamental harmonics is the solvability condition in the method of multiple scales: it removes the resonant forcing that would otherwise generate a secular term. Since the squared sine and cosine averages are $1/2$, the amplitude-phase equations areEquivalently $\dot{R}=-\epsilon ⟨f\sin (t+\theta )⟩$ and $\dot{\theta }=-\epsilon ⟨f\cos (t+\theta )⟩/R$ at first order. These equations describe a bounded, weakly perturbed oscillation over $t=O({\epsilon }^{-1})$ while the amplitude remains in a range where the expansion is ordered. At zero amplitude the oscillation phase coordinate is singular; Cartesian harmonic coefficients or an equilibrium analysis should replace it.

For conservative forcing in averaged oscillator equations, $f=f(y)$, the average determining $R_T$ vanishes. If $F^{\prime }=f$, thenwhose period average is zero. Thus $R=R_0$ to this order, while a position-dependent force generally shifts the oscillation phase and frequency. The exact equation also conserves the energy $\frac{1}{2}{\dot{y}}^2+\frac{1}{2}{y}^2-\epsilon F(y)$, consistent with the absence of amplitude drift for bounded orbits.

For the cubic position force, $⟨f\cos \varphi ⟩=R^3⟨{\cos }^4\varphi ⟩=3R^3/8$. ThereforeThe leading oscillation is $R_0\cos [(1-3\epsilon R_0^2/8)t+{\theta }_0]$. For positive $\epsilon$ it is a softening cubic oscillator, obtained by weakly perturbing a harmonic oscillator. The weak-force condition is $\epsilon R_0^2\ll 1$; the constant $R$ is the leading harmonic amplitude, not a claim that the exact waveform remains sinusoidal.

For velocity-only forcing in averaged oscillator equations, $f=f(\dot{y})$, the oscillation phase average vanishes: $f(-R\sin \varphi )\cos \varphi$ is a periodic total derivative, using an antiderivative of $f$. Hence ${\theta }_T=0$ to first order, although the amplitude may change.

In the cubic case, $f=-R^3{\sin }^3\varphi$, soWith $R(0)=R_0>0$, integration givesFor positive $\epsilon$ this is cubic velocity anti-damping, not damping: the exact energy derivative is $d({\dot{y}}^2+y^2)/(2dt)=\epsilon {\dot{y}}^4\ge 0$. The slow solution grows and formally diverges at $T=4/(3R_0^2)$. The weakly nonlinear approximation loses validity before its amplitude becomes so large that $\epsilon R^2$ is not small; its divergence is not a controlled prediction of the exact late-time solution. If $R_0=0$, the exact zero solution should be used separately.

Let $\delta =1-m\downarrow 0$. This is the complete elliptic integral of the first kind in parameter notation; the modulus used in some definitions is $\sqrt{m}$. Near the endpoint, $s=\pi /2-\vartheta$ givesUse matched asymptotic expansion with an intermediate cutoff $b$ satisfying $\sqrt{\delta }\ll b\ll 1$. Away from the endpoint the leading integral iswhile the endpoint integral isAdding them removes the arbitrary cutoff and gives the logarithmic endpoint asymptotic of the complete elliptic integralThe order of the next term is therefore $(1-m)\log (1/(1-m))$, rather than merely $1-m$. More explicitly, if $L=\log (4/\sqrt{\delta })$,The logarithmic correction arises from the next terms integrated through the overlap. Its coefficient can also be found by substituting $L+\delta (aL+b)$ into the Gauss hypergeometric equation satisfied here, $m(1-m)K_{mm}+(1-2m)K_m-K/4=0$, giving $a=1/4$, $b=-1/4$.

For fixed $x>0$, set $y=y_0+\epsilon y_1+{\epsilon }^2{y}_2+\cdots$ and impose the data at one at each order. The leading equation gives $y_0^{\prime }=0$, hence $y_0=1$. At the next two orders,The boundary conditions $y_1(1)=y_2(1)=0$ give the three-term outer expansionIts ordering fails when $\epsilon /x=O(1)$, locating a boundary layer of thickness $O(\epsilon )$ at zero. Put $x=\epsilon \xi$ and $y=Y_0(\xi )+\epsilon Y_1(\xi )+\cdots$. The rescaled equation isIts leading equation is ${\xi }^2{Y}_0^{\prime }=Y_0^3$. Matching to the positive outer value one fixes the sign and integration constant:At the next order the bracket $\xi (Y_0^2-1)+2Y_0^2$ vanishes identically, leavingThe general solution is $Y_1=(1+k\xi ){\xi }^{1/2}/(\xi +2{)}^{3/2}$. In the overlap its expansion isThe outer solution re-expressed on this scale has the corresponding terms $1-1/\xi +3/(2{\xi }^2)+\epsilon (1-2/\xi )+\cdots$. Thus matched asymptotic expansion determines $k=1$, consistently in both displayed orders, and the two-term inner expansion isThis positive branch tends continuously to zero at $x=0$, with a square-root cusp. An infinite endpoint derivative is compatible with the degeneracy of the original equation's $x^2{y}^{\prime }$ coefficient.

At fixed positive $x$ the leading outer expansion satisfies $z_0^{\prime }+z_0=e^{-x}/x$. Integrating and imposing the value at one yieldsIt diverges negatively at zero, and the neglected nonlinear shift in the derivative coefficient becomes important when $x\sim \epsilon |\log x|$. Thus the relevant logarithmically enhanced nonlinear boundary layer is larger than a plain $O(\epsilon )$ layer. Write $L=\log (1/\epsilon )$, $x=\epsilon LX$ andFor fixed $X$, the leading inner equation is $(X+1)Z_X=1$. Matching to the outer logarithm, for which $Z\sim \log X$, sets the integration constant and givesIn particular the exact equation at zero is $-\epsilon z(0)z^{\prime }(0)=1$, so

An implicit inner form also checks the matching constants. With $x=\epsilon \xi$, retain $z$ without assigning it a bounded size. The leading equation is $(\xi -z)z_{\xi }=1$. Inverting it gives $d\xi /dz-\xi =-z$, henceIts large-$\xi$ match to $z\sim 1+\log \epsilon +\log \xi$ fixes $B\sim e^{-1}/\epsilon$. Taking that matched leading value at zero giveswhere $W_0$ is the positive real branch of the Lambert W function. Its large-argument expansion reproduces the logarithm and log-log terms above. This refinement is a matched approximation, not an exact solution of the full equation; the leading slope conclusion follows directly from the first inner scaling and the exact endpoint identity.

Keep the prescribed outer slip velocity fixed when perturbing the unsteady Prandtl equation. Cancelling the base equation and retaining terms linear in the disturbance gives the linearized unsteady Prandtl equationThe no-slip boundary condition supplies $\tilde{u}=\tilde{v}=0$ at the wall; the fixed outer velocity supplies $\tilde{u}\to 0$ at infinity. There is no prescribed zero normal-velocity disturbance at infinity: a finite displacement-related value is permitted.

With the frozen parallel profile $U=U(y)$, $V=0$, the coefficients are invariant under translations in $x$ and $t$. Normal modes therefore separate those variables. For the mode convention used here the equations areEliminate $u=-v^{\prime }/(i\alpha )$ to obtain the Prandtl normal-mode equationwith a bounded far-field $v$. Requiring $v(\infty )=0$ would incorrectly eliminate the proposed outer profile. For real $U$ and real $\alpha$, complex conjugation of the equation and its boundary conditions replaces $(\alpha ,c,v)$ by $(-\alpha ,c^{\ast },v^{\ast })$. This proves the stated eigenvalue symmetry. The paired modes have the same temporal growth rate $\alpha \mathrm{Im}c$, so choosing positive wavenumber loses no real physical disturbance.

Here the prescribed parallel profile is a local frozen-coefficient model. A nontrivial arbitrary $U(y)$ with $V=0$ and constant outer slip is not generally an exact steady solution of the unforced unsteady Prandtl equation; its base equation would require $U^{\prime \prime }=0$. The following stability calculation uses the simplifying local model stipulated for the mode analysis.

Set $c_0=U_c=U(y_c)$. Away from the critical point the dominant equation is $(U-c_0)v_0^{\prime }-U^{\prime }{v}_0=0$, whose solutions are multiples of $U-c_0$. Choosing the coefficient to be zero below $y_c$ and one above it satisfies the wall and tangential far-field conditions. The choice $c_0=U_c$ makes $v_0$ continuous across the joining point; $U^{\prime }(y_c)=0$ also makes its first derivative continuous. Its second derivative jumps. Thus it is a valid leading outer solution away from $y_c$, but viscosity must smooth the join in a critical layer in a shear flow.

Near the join, $U-c_0\sim (y-y_c{)}^2/2$. If the layer width is $\ell$, matching gives $v=O({\ell }^2)$. The advection side of the mode equation scales as $\alpha {\ell }^3$ and $v^{\prime \prime \prime }$ scales as ${\ell }^{-1}$. Balance gives ${\ell }^4\sim {\alpha }^{-1}=\epsilon /k$. Equivalently, the eigenvalue correction ${\epsilon }^{1/2}{c}_1$ balances $(y-y_c{)}^2$, so the quarter-power critical layer hasSubstitution, retaining the leading terms and cancelling their common factor, gives

At order ${\epsilon }^{1/2}$ outside the layer,A solution consistent with the chosen lower branch and the wall conditions is $v_1=0$ below the join. Above it a particular solution is $v_1=-c_1$; a multiple of $U-c_0$ represents an arbitrary amplitude renormalization and may be set to zero. Matching consequently requiresThe corresponding derivatives match as $w_Y\sim Y$, $w_{YY}\to 1$ on the positive side and as zero on the negative side. The matching statements require that growing homogeneous corrections be absent.

To symmetrize these different end conditions and remove $k$ and $i$, chooseThe polynomial $A(z)=z^2/2-C$ is itself a solution of the transformed homogeneous equation. Direct substitution, using $a^4=ik$, givesThe derivative is third order, as required by the original mode equation. Initially the two ends are along the rotated rays $z=aY$ with real $Y$. Continuing them to the real $z$ axis is compatible with the asymptotic end conditions: the decreasing homogeneous correction behaves exponentially as $\exp [-z^2/(2\sqrt{2})]$ and still decreases throughout the rotation from angle $\pi /8$ to zero. Thus the supplied real-axis spectral normalization uses the same recessive conditions; a complex coordinate change should not be mistaken for a real stretching alone.

For each supplied real eigenvalue $C_n=(4n+7)/\sqrt{2}$, the phase-speed correction isThe mode factor has temporal magnitude $e^{\alpha \mathrm{Im}ct}$. ThereforeAmong the listed indices, $n=-2$ gives positive growth $\gamma \sim \sqrt{\alpha }/2$; $n=-1,1,2,\dots$ give decay. Hence the frozen non-monotone flow has high-wavenumber instability of a non-monotone Prandtl layer, with arbitrarily rapid linear growth as the positive wavenumber increases. The normal-mode calculation demonstrates linear instability; it is not by itself a claim about the final nonlinear state.