On the Abel Equation of the Second Kind with Sinusoidal Forcing

Abstract. We consider the Abel equation of the second kind with sinusoidal forcing, which is a model equation for western boundary outflow in the Stommel mod el of the large scale ocean circulation. Series solutions of this equation indicate the pr es nce of resonances at certain discrete values of a parameter which measures the nonlinearity of the flow, but numerical solutions using a standard scheme show no evide nce of these resonances. We discuss and resolve this apparent contradiction.


Introduction
The large scale circulation of the oceans is both one of the most important and complex problems arising in geophysical fluid dynamics.One of the features of this circulation of particular interest over the years has been the flow near the western boundary [1,9] because of the obvious importance of the Gulf Stream in the North Atlantic and similar phenomena in other oceans.Since the Gulf Stream separates from the US east coast at Cape Hatteras, one particular aspect of the western boundary current that has been of interest to oceanographers is the issue of when and how it separates from the boundary.We should point out that in this paper we use the term separation to include a separation bubble or recirculation gyre as well as a fully detached boundary layer.Even some very simple models can mimic this separation, and in this study, we revisit one of the simplest of these models, that due to Stommel [10] which consists of two-dimensional flow in a rectangular basin on a beta plane with the frictional forces represented by an effective bottom drag −r∇ 2 ψ, so that the stream function ψ obeys the time-independent potential vorticity equation on a beta plane, where β is the Rossby parameter, W is the curl of the wind stress, and r is the friction parameter.The analysis presented here stems from the Stommel model, but we note that a slightly more complicated (and realistic) model is that due to Munk [7], in which the frictional forces are represented by lateral friction Re −1 ∇ 4 ψ.Although the Munk model is a more realistic model of the circulation of a real ocean, it is less tractable analytically.
Our experience [5] is that the boundary layers in the two models behave in a somewhat similar manner, and therefore it is likely that analysis for the Stommel model will provide additional insight to the Munk model.
In both the Stommel and Munk models, frictional forces are usually neglected in the interior of the flow, but become important near the walls, particularly near the western wall, where a thin viscoinertial boundary layer is necessary in order for the velocity to satisfy the no-normal flow condition at the wall.In this layer, there is a threefold balance between the Coriolis term, the viscous term and the nonlinear term.Usually, the flow in the interior is fairly simple, although one study [2] for the Munk model has indicated that this is not always the case.
In what follows, we shall consider flow in a rectangular basin, x W ≤ x ≤ x E and −L ≤ y ≤ L, where physically, the direction of increasing x corresponds to east and that of increasing y to north, and we shall assume for simplicity that the flow in the interior is zonally uniform, meaning independent of the east-west coordinate [11], with the stream function in the interior given by ψ I = 2 cos (πy/(2L)) and the curl of the wind stress by W = −W 0 cos (πy/(2L)), with W 0 = rπ 2 /(2L 2 ).This stream function/wind stress combination is a solution of (1), but does not satisfy the conditions at either the eastern or western walls, so that a boundary layer is required at each of these two walls.Traditionally, the boundary layer analysis for this problem has been concerned with the flow in the northwestern corner, as that is where separation first occurs.and therefore we expand about the point (x W , L), introducing the stretched coordinates, ξ = (x − x W )/δ, η = (y − L)/L and Ψ(ξ, η) = ψ(x, y), where δ = r/β is the boundary layer thickness, which is found by balancing the viscous and Coriolis terms inside the boundary layer.With these scalings, the potential vorticity equation inside the boundary layer becomes where λ = π/(βLδ 2 ) is a measure of the strength of the nonlinearity and ε = (δ/L) 2 measures the relative sizes of the derivatives with respect to ξ and η in the Laplacian.
If we make the usual boundary layer approximation [12][13][14][15] and let the boundary layer thickness δ → 0, we get The boundary conditions for Ψ, the stream function inside the boundary layer, are that Ψ is constant on both ξ = 0 and η = 0, and also that Ψ → ψ I as ξ → ∞, so that the flow inside the boundary layer matches onto the flow in the interior.If we use these boundary conditions, we can integrate (3) with respect to ξ, If we evaluate (4) at the western wall ξ = 0, where Ψ(0, η) = Ψ η (0, η) = 0, it reduces to where V (0, η) = V 0 (η) is the northward velocity V = Ψ ξ evaluated at the wall.The condition that Ψ is constant on η = 0 becomes V 0 (0) = 0.This is an Abel equation of the second kind, V 0 V 0ζ − V 0 = f (ζ), as discussed in §1.3 of [18].The better known solvable cases are listed in [18], and this does not appear to be one of them.
In [6], we used (5) as a simple model equation for western boundary outflow in the Stommel model.In that paper, we sought series solutions to (5) of the form V 0 (η) = ∞ j=1 v j η j and also to (4) of the form Ψ(ξ, η) = ∞ j=1 φ j (ξ)η j , with the two series related by v j = φ ′ j (0).We found that the model equation ( 5) captured some of the features of (4).For example, for both equations, solutions of the form sought could only exist if the parameter λ was less than a critical value λ c .The critical value λ c was reported in [5,6,8] for the Stommel model and in [3,4] for the Munk model.The failure of these equations to have a solution is coincident with the separation of the boundary layer [16], and in [17], it was shown that as λ is increased beyond λ c , the separation point moved further south.
A second feature that the series solutions to ( 5) and ( 4) shared was the presence of resonances at certain discrete values of λ, with the series solutions becoming singular as λ approached the resonant values, and those resonances are the topic of the present study, where we will pursue the asymptotics of (5) in more detail than in [6], using the power of symbolic computation, and also compare the series solutions to numerical solutions.

Series solution
We need to find solutions to (5) which vanish at η = 0, and to this end, we will try a series solution of the form If we substitute this series into the Abel equation ( 5) and group powers of η, we get a hierarchy of equations.At leading order O (η), we find the coefficient v 1 must obey which has two solutions which are plotted in Fig. 1(a), where we take the "+" sign for v (+) 1 and the "−" sign for v (−) 1 .In order for v 1 to be real, which is required, the quantity under the square root in (8) cannot be negative, so we require λ ≤ 1  4 for a solution of the form (6) to exist.In what follows, we will write the coefficients in the series (6) in terms of the parameter µ which we used in [5,6,17], which is related to λ by λµ 2 − µ + 1 = 0, or λ = (µ − 1) /µ 2 , which has two roots which are plotted in Fig. 1(b), so that v 1 = −πµ. -20 λ At the next order O η 2 , we find that v 3 must obey This means that µ = 3 2 is a special case which must be considered separately.
If we continue in this fashion, at O(η n ), we find that the equation for v n is of the form If n = 2m is even, ( 14) will be of the form 14) will be of the form where c 2m+1 = 0, so that we get a contradiction if µ = 1 + 1/(2m + 1), meaning that additional terms must be included in the series solution for that value of µ.
The upshot of this is that if µ = 1 + 1/n for n = 1, 2, 3, . .., then there is a resonance in the series solution and additional terms must be included, while if µ = 1 + 1/n for n = 1, 2, 3, . .., then the series solution is of the form with These coefficients appear to be singular as we approach µ = 1 + 1/(2m + 1) for m = 1, 2, 3, . .., so we might expect that a plot of V 0 (η) as a function of µ for a given value of η would be singular at the values of µ.To examine this, we used a computer algebra software package to calculate the coefficients in the series (17) as far as v 55 , and we then evaluated the series at η = 0.1, keeping terms as far as v 55 η 55 , again using the computer algebra package and keeping 150 significant figures in the calculations.The results are shown in Fig. 2(a), where we plot V 0 (0.1) as a function of µ, and it can be seen that as expected, the series solution has very narrow singularities at µ = 1 + 1/(2m + 1) for n = 1, 2, 3, . .., and this figure would seem to confirm that V 0 (η) has singularities at those values of µ.In Fig. 2(b), we plot the series for values of µ close to 4  3 , which is one of the resonant values of µ, again keeping terms as far as v 55 η 55 , and it can be seen that the series blows up, with the blow-up happening earlier as µ approaches 4  3 .
- For the resonant values of µ, additional terms must be included in the series.If µ = 1 + 1/(2m), for example 3  2 , 5 4 , 7 6 , . .., then even powers of η are required, so that the series solution will be of the form with v 2j = 0 for j < m while v 2m is a free parameter.For example, for µ = 3 2 , we have where v 2 is a free parameter.If µ = 1 + 1/(2m + 1), for example 4 3 , 6 5 , 8 7 , . .., then logarithms are required, so that the series solution will be of the form where v 2m+1 is a free parameter.For example, for µ = 4 3 , we have where v 3 is a free parameter.Perhaps surprisingly, µ = 2, which corresponds to the critical value of λ = 1 4 in ( 9), is not a special case and is covered by the series solution ( 17), (18).

Numerical solution
The analysis of Section 2.1 would seem to indicate that the Abel equation might be expected to exhibit resonant behavior when µ is equal to one of the resonant values discussed in Section 2.1, and to explore this more fully, we solved (23) numerically using a standard fourth order Runge-Kutta scheme, using the initial condition V 0 (0) = 0.These numerical solutions would seem to indicate that the resonances do not occur.In Fig. 2(a), we plotted V 0 (0.1) as a function of µ using the series solution.In Fig. 3(a), we repeat this plot but using the numerical solution rather than the series solution.For µ ≤ 1, to the left of all the resonances, and µ greater than about 1.35, just to the right of the resonance at µ = 4 3 , the plots from the series solution and the numerical solution are near identical.However, for µ between 1 and about 1.35, the numerical and series solutions are somewhat different.Between these values, V 0 (0.1) from the series solution appears to lie on a smooth curve between the branches outside this range, except at the resonances, at each of which there is a highly localized spike.By contrast, the numerical V 0 (0.1) tends to −∞ as µ → 0+, with nothing untoward happening at the resonances.
In Fig. 2(b), we plotted the series for values of µ close to the resonance at 4  3 .In Fig. 3(b), we repeat this plot but using the numerical solution rather than the series solution, and again the results are somewhat different.The numerical solution does not blow up, and all of the curves in Fig. 3(b) are very similar.
Taken together, Fig. 3(a) and 3(b) would seem to suggest that the resonances do not in fact occur.In Fig. 4, we compare the numerical and series solutions to (23) close to µ = 1.For µ slightly less than 1, in this case µ = 0.9, the series and numerical solutions are in excellent accord, but for µ slightly greater than 1, the series and numerical solutions are completely different.
There is a simple reason the two solutions are so different in the region 1 ≤ µ ≤ 1.35, which is that the series ( 17) is failing to converge in this region.The coefficients (18) in (17) involve negative powers of (µ − 4 3 ), (µ − 6 5 ), (µ − 8 7 ), . .., and in Table 1, we give the exponents of these terms for the first few coefficients in the series.The exponent of (µ − 4 3 ) decreases by 1 every coefficient, that of (µ − 6 5 ) every second coefficient, and that of (µ − 1 − 1 2m+1 ) every mth coefficient.Because (17) involves only odd powers of Table 1.Exponent of (µ − µ0) the coefficients (18) in the series (17) near the resonances at µ0 = 2m/(2m − 1) η, we would expect that the radius of converge of (17) behave like (µ − 4 3 ) 1/2 near the resonance at µ = 4  3 , like (µ − 6 5 ) 1/4 near µ = 6 5 , and like (µ − 1 − 1 2m+1 ) 1/(2m) near µ = 1 + 1 2m+1 , so that as µ approaches a resonance, the radius of convergence of (17) will approach zero.In this region then, the series can tell us the value of V 0 and each of its derivatives at η = 0, but is of little use as we move away from η = 0. Outside the region, there is good agreement between the series and numerical solutions.However, since the resonances at µ = 1 + 1 2m+1 lie inside that region, the region where the series is of little use is the region which is of most interest.
With regard to the solution at the resonances themselves, the numerical solutions discussed above showed no sign of them.However, this is a limitation of the numerical method used, and is an instance where the absence of evidence is not evidence of absence.The Runge-Kutta method used estimates V ′ 0 (η) at 4 points in each time step, with occur because the higher derivatives are undetermined, and to detect them, the numerical method must involve those higher derivatives.For example, the resonance at µ = 3 2 occurs because V ′′ 0 (0) is undetermined, and one way to investigate it is to differentiate the Abel equation ( 23) once with respect to η to get a second order equation, If we set µ = 3 2 and then introduce which we can again solve using a Runge-Kutta method by estimating both V ′ 0 (η) and V ′ 1 (η) at 4 points in each time step, with for η > 0, where v 2 is a parameter which must be specified.A similar approach, involving successively larger systems, can be used to investigate the resonances at µ = 5 4 , 7 6 , . ... In Fig. 5(a), we plot the numerical solutions using this approach.The solution using the original method, which is included for comparison, coincides with the new method when v 2 = 0.The series solution (19) is plotted in Fig. 5(b), and the solution with v 2 = 0 coincides with the general series (17).With the exception of the cases v 2 = ±10, for the cases shown, there is good agreement between the numerical results and the series up to about η = 0.7, when the series diverge.This can be taken as numerical confirmation of the existence of the resonances and the resultant multiple solutions.
The resonances at µ = 1 + 1 2m+1 occur because the higher derivatives are singular, and to detect them, we must remove the singularity at leading order.For example, for the resonance at µ = 4  3 , the series solution was so if we define V 0 = W 0 − π 3 6 η 3 ln η, then the series solution for W 0 will be so that the leading order logarithm has been removed.If we differentiate the Abel equation ( 23) twice with respect to η and set η = 4 3 , and then introduce W 1 (η) = W ′ 0 (η) and W 2 (η) = W ′ 1 (η), we have a first order system which we can again solve using a Runge-Kutta method by estimating each of , where v 3 is a parameter which must be specified.A similar approach, involving successively larger systems, can be used to investigate the resonances at µ = 6  5 , 8 7 , . ... In Fig. 6(a), we plot the numerical solutions using this approach.The solution using the original method is included for  comparison, but does not correspond to any value of v 3 .The series solution (21) is plotted in Fig. 6(b).The general series ( 17) was singular at µ = 4 3 .Once again, for the cases shown, there is good agreement between the numerical results and the series for small η, and again, this can be taken as numerical confirmation of the existence of the resonances and the resultant multiple solutions.

Discussion
The results of this paper are a cautionary tale.In the previous section, we found both series and numerical solutions to the Abel equation of the second kind with sinusoidal forcing (23), which is a model equation for western boundary outflow in the Stommel model of the large scale ocean circulation [6].The series solution has resonances at certain discrete values of the parameter µ, at which values (23) has multiple series solutions and additional terms need to be included in the series.As µ approaches these resonant values, the series solution becomes singular.The numerical solutions tell a somewhat different story.When we solved (23) using a standard fourth order Runge-Kutta scheme, it was apparent that the numerical solutions were not singular near the resonant values, and in this conflict between the numerics and the asymptotics, the numerics came out on top because as µ approaches a resonant value, the radius of convergence of the series solution approaches zero.The standard numerical scheme also seemed to find no evidence of the resonances and multiple solutions predicted by the asymptotics, and on this issue, the asymptotics came out on top.The resonances occur because the higher derivatives are either undetermined or singular, and the standard numerical scheme was oblivious to these resonances because it was calculating only the first derivative.Once we converted the equation to a system of equations involving these higher derivatives, the numerics were able to verify the presence of the resonances and the existence of the multiple solutions, but without the information gleaned from the series solution about the location and nature of the resonances, we would not have been able to construct the appropriate numerical schemes and indeed we would have been as oblivious to the resonances as the original numerical scheme.The moral of this cautionary tale is that both series solutions and numerics can provide valuable information but both need to be used with caution: series which do not converge provide misleading information while numerics can fail to provide important information.