Renormalized Coordinate Stretching: A Generalization of Shoot and Fit with Application to Stellar Structure

The standard shooting and fitting algorithm for non-linear two-point boundary value problems derives from conventional coordinate perturbation theory. We generalize the algorithm using the renormalized perturbation theory of strained coordinates. This allows for the introduction of an arbitrary function, which may be chosen to improve numerical convergence. An application to a problem in stellar structure exemplifies the algorithm and shows that, when used in conjunction with the standard procedure, it has superior convergence compared to the standard one alone.


Introduction
Consider the two-point boundary value problem governed by a non-linear system of k ordinary differential equations:  y f y (1) with k boundary conditions distributed at each end of the domain in x.Such cases may be solved numerically by varying the unknown conditions at each boundary and integrating each trial to a common fitting point, whence differences between the trial runs provide linearly independent functions by which to calculate corrections to the initial guesses [1][2][3].This method of shoot and fit is based on conventional perturbation theory, which expands the dependent variables: where   1 x y is a correction to   0 x y .The method of strained coordinates was first introduced by Lindstedt (1882) and other astronomers and now known as the Poincaré-Lighthill (PL) or Poincaré-Lighthill-Kuo method [4][5][6][7][8].It is well-known that the PL method was designed originally for systems   d d , ; x x   y f y that depend critically on a small parameter  , in order to render analytic series solutions uniformly convergent.However, we apply it here to problems like Equation (1) whose formulation is formally independent of a small parameter and for which standard coordinate perturbation methods are applicable.We use the adaptation in its renormalized form [9,10], which Dai [11] considers to be one of five important improvements in PL theory.Thus, we distinguish between standard and renormalized PL theory.
Section 2 shows the utility of the renormalized PL (rPL) theory in rendering variational equations homogeneous.Section 3 formulates the standard (std) fitting algorithm and Section 4 formulates its rPL generalization.Section 5 exemplifies the algorithm by means of a twopoint problem of stellar structure formulated in the Appendix, and compares the utility of each algorithm and the possibility of their combined usage.Section 6 discusses results.

Straining Coordinates
Apply Equation (2) to Equation (1) and separate orders to give to first order: Following the standard PL procedure, x may be expanded as well as y(x), each as a function of a new independent variable 0     Function   1 0 x x stretches the independent coordinate x.Equations ( 5) and ( 6), when applied to Equation (1), give to zeroth and first order: However, expansions ( 5) and ( 6) are not unique.Instead, apply Equation ( 6) to the equivalent standard expansion (2) [10], i.e. let:     so that Equation ( 9) becomes: On applying Equations ( 10) and (11) to Equation (1), and using the relation: we discover that all terms in   1 0 x x cancel and: These are equivalent to Equations ( 3) and ( 4), except that the independent coordinate is x 0 rather than x.The straining function   1 0 x x now appears only in expansions (10) and (11), and not in the differential Equation (8).On retaining terms to first order, these become the rPL perturbation expansions to first order:   Implementation of Equations ( 13) and ( 14), is straightforward as the conventional theory is recovered simply by setting   1 0 0 x x  , and x 0 reverts to its original meaning x, as in Equations ( 3) and (4).In other words, in rPL, the straining feature may be applied after, and not during, the integration of the variational equations.

Standard Shoot and Fit
Frequently, in two-point boundary value problems, it is impossible to integrate from one boundary to the other without encountering some catastrophic result.This might occur while integrating toward a singular boundary, or when solutions become imaginary or unphysical.The standard (std) shoot and fit algorithm is designed to counter these difficulties, and in anticipation of its rPL generalization it is useful to spell it out.
Assume disjoint explicit boundary conditions, a normalized range, and for clarity distinguish integrations from each boundary by lower case and upper case variables.Throughout, let For integrations from x = 0, Equation ( 13) is: and from X = 1, it is: x for 0 x    , and    , using initial condi- tions chosen as: where

and ,
ii jj     are Kronecker deltas.Assume that the differences: are linearly independent, where of course from Equations (17), ( 19), ( 21)-(24): Linear superposition of Equations ( 23) and (24) gives: where and i I j C C  are constants to be determined.Equations (2), ( 27) and (28) lead to: Take the computed mismatches at ξ to be: which, from the difference of Equations ( 27) and ( 28) at ξ and the desired outcome These are I J  linear algebraic equations in I J  unknowns, whose solution gives and
must satisfy two constraints: 1) In order that the independent variable be continuous at 0 X  , and vice versa; i.e.
2) With reference to Equations ( 39) and (40), 0 and 1 x X must be such as to nullify any singularities in     0 0 0 and 1 f F .

Stretching
Functions are otherwise arbitrary, but in the spirit of linearization let them vary linearly with their arguments.Equations ( 36) and (38) become: Let the value in Equation (42) be: where 1 I J C   is the amount of stretching/compression at the matching point.Thus:   and Equation (41) becomes: These are I + J equations in I + J + 1 unknowns, which require an extra equation for their solution.

Matching Condition
Let the (I + J + 1) th equation follow from matching s(y,x) and its functional twin S(Y,X) at the fitting point ξ.Expansions using Equations (39) and (40) are: The desired outcome s(ξ) -S(ξ) = 0 and the computed difference: lead to: .
The I + J + 1 constants follow from the solution to Equations (48) and (52).With reference to Equations (39) and (40), improved guesses are: By Section 4.1, the last terms of Equations ( 53) and (54) must be finite or zero.If 0 and 1 x X nullify these terms, then the rPL Equations ( 53) and (54) resemble Equations ( 33) and (34) of std.In that case, the influence of stretching on the improved boundary conditions is nominally independent of C  through Equations ( 48) and (52).

Algorithm Option
For given trial values, relatively little extra computation time is involved in computing and solving the I + J + 1 algebraic Equations ( 48) and (52), compared to I + J Equations (32), whereupon the option exists to complete the next iteration using the rPL or the std algorithm.This requires devising a criterion by which to choose iterated values from one algorithm over values from the other.
From Equations ( 33) and (34), let the relative std corrections be: From Equations ( 39) and (40), let the relative rPL corrections be: As noted, quantities  are proportional to the stretching constant 1 I J C   and must be zero or finite.Define the sum of the squares of the relative corrections as: and let its value decide between the two options, i.e. between Equations ( 55) and (56), and Equations ( 57) and (58).Since relative errors obtain in Equation (59), it applies regardless of whether the converged solution is known (as it is in the present discussion), or not; i.e.Equation ( 59) is useful provided a way exists by which to decide between the available options.Section 5 contains empirical evidence concerning the better option.In practice, sequences of iterations generally have corrections derived from both algorithms; i.e. the two algorithms have the potential to widen the circle of convergence, since when one fails for any given iteration, the other may not.Thus, at the point of a failed iteration, an additional integration is the price to pay for the prospect of a wider circle of convergence.

An Illustrative Example
The shoot and fit method is applicable to an assortment of problems, including the two-point boundary problem of stellar structure [1,12].Coordinate stretching has also proved useful in the study of rotating stars [13].Stellar structure appears to be a suitable means for testing the rPL method against a known solution, for which it suffices to choose a relatively simple yet realistic model for a thirty-solar-mass (6 × 10 31 kg) star.This is identical to early models [14,15], which were constructed less precisely and by different algorithms.The underlying physics is stated in the Appendix.From transformations of Equations ( 76)-( 81), ( 84), (85) of the Appendix, consider the problem:     where  is an index (see Appendix).Equation ( 63 Physical arguments require The unit boundary values of Equation (66) across the normalized range of the independent variable may be used here because the intent of Section 6 of this paper is to investigate the relative efficacy of algorithms by comparing trial values against the known solution.For the fourth-order system considered, the solution is characterized by four eigenvalues that correspond to four normalized ranges for the set   , , , z z z z .An auxiliary eigenvalue q marks the transition from one of Equations ( 63) and (64) to the other.Integration proceeds from each boundary to an arbitrarily chosen fitting point ξ, chosen throughout as 0.5.Integrations used the Odesolve (vector, x, b, npoints) routine of MATHCAD 14 [16], where vector is the array of sought solutions, x is the independent variable, b its terminal value, and npoints is the integer number of points used to interpolate the solution functions.Odesolve dynamically detects whether a system is stiff or not, and uses BDF (backward differentiation formula) or AB (Adams-Bashforth) methods respectively.Both methods use non-uniform step sizes, adding more integration steps in domains of greater variation.Subroutine root with tolerance TOL = 10 −9 determined the location of the interface fraction q.The default precision is 17 significant figures.
Singular gradients at the boundaries require asymp-totic developments, here limited to one significant term owing to the high precision of the calculations.Following Sections 3 and 4, lower and upper case variables distinguish integrations from x = 0 and x = 1.For small δx, δX, starting values are: Figure 1 shows the normalized solution.Its derived boundary values and the auxiliary eigenvalue q agree to 1% or less with the aforementioned 30 solar mass model, with discrepancies explained by fitting by interpolation or by hand used previously.Section 6.1 applies the std and rPL formulations.Integration proceeds through interfaces brought about by changes between Equations ( 63) and (64), and terminate at ξ = 0.5.To test the sensitivity of each independent variable to guesses in its unknown boundary value, we assign in turn the correct value 1 to all but one variable per Equation ( 65), to which we give values progressively farther from 1. Trial values are changed at intervals of ± 0.001 in the log, and convergence is considered continuous until the iteration fails to converge.Sometimes, convergence resumes for further changes in the guessed value, but it suffices for our purposes to apply the criterion of first failure uniformly across all cases.Execution ceases when one (or two) of the corrected guesses turns negative or when the number of iterations exceeds 20.
For the extra matching condition s(x), we require the continuity of a combination of variables unaccounted for by the original system.Consider the 6 products and require: None of these conditions is explicitly required in the problem formulation.Table 1 compares results for the 6 cases that use these rPL constraints simultaneously with the standard std algorithm.Iterations ceased upon satisfaction of the convergence criteria: for i = 1,2 and j = 3,4.Results for entry 2 of Table 1 appear in Figures 2(a)-(d).
Devise two statistics by which to assess overall relative performance.For every i z , the extrema of the positive Choose the other statistic as the relative rms measure:

Discussion
The std algorithm sets the standard by which to judge the efficacy of the rPL algorithm.It appears as Case 1 in Table 1, for which by definition ψ = 1.When the rPL algorithm alone is applied, only in about half the cases there is there an improvement over the std; i.e. for the present experiment, there is not much to choose between the two.Little is gained by showing these cases.However, strength of the rPL lies in the similarity of Equations ( 4) and ( 14), which allows application of the rPL algorithm at little extra cost in computing time.We discover empirically (and counter-intuitively) that the better choice has the lesser value of Δ 2 in Equation ( 59).This appears to minimize the tendency for iterations to diverge by over-correcting.Under these conditions, when the std option is tested simultaneously with the rPL (Cases 2-7 of Table 1), there is an improvement in convergence characteristics in all but one case.In Case 2, convergence occurs for as little as ≈3 × 10 −2 of the correct value.This phenomenon results from the fact that the model investigated is only a factor of 8 less than the Eddington luminosity limit at which stars are dynamically unstable (see Appendix), so that underestimates are computationally beneficial.
The excellent statistics of Case 2 in Table 1 attract attention.Figures 2(a)-(d z z , where 1 z and 2 z happen to be the two variables with the least uniformly-changing traces, as seen in Figure 1.For the present experiment, the comparative pathologies turn out to be virtual double-mirror images of one another and the product is quite well-behaved.It would seem a priori that these two variables would benefit most from an independent condition of continuity that involves them.
The chief result based on the problem definition and methodology here articulated, is that coordinate stretching

Appendix
Consider a non-rotating, radially symmetric star of mass M = 6 × 10 31 kg (30 solar masses) in quasi-hydrostatic equilibrium.Suppose that it is chemically homogeneous and that thermonuclear conversion of hydrogen to helium has just begun and is in a steady state.Assume that the material is entirely ionized and that it obeys the ideal gas law throughout.Let the relative mass abundances of hydrogen, helium and heavier elements be X = 0.70, Y = 0.27, and Z = 0.03, and let the mean molecular weight of the fully ionized gas in units of the proton mass H be μ Constants are k (Boltzmann's constant), c (the speed of light), a = 8π 5 k 4 /15c 3 h 3 (the radiation density constant), h (Planck's constant) and G (the gravitational constant).Variables are: v (volume), r (radius), ρ (material density), ℓ (luminosity or power), p total pressure, equal to the sum of radiation pressure , where ρ is the material density, t temperature, and m mass fraction, i.e. the mass contained within volume v.In many applications, it is customary to choose m as the independent variable.Gradients are then: where Equations ( 79) and (80) pertain to radiative and adiabatic convective energy transport for which the lesser of |dt/dm| applies.Physical conditions set limits 0 < β < 1.The ratio of gas specific heats at constant pressure and volume is .Take the rate of energy generation to be , where ν = 15 and 0  = 10 -113.6 Jm 3 ·kg -2 ·s -1 .Nomenclature in this formulation is standard, and symbols ε, k, X, Y, and Z are not to be confused with occurrences out of context.
During computation, starting values δx and δX were diminished until there were no further changes in the converged solution to a relative accuracy <10 -6 .For trial integrations sufficiently close to the true solution, we discover through trial and error or otherwise that radiative equilibrium obtains in the envelope and convective equilibrium in the core.Near the center 0 m  and Equations ( 76) and (77) give , where it suffices to let p and t vanish.Equations ( 78) and (79) give:     Thus, boundary conditions are: The stretching functions in Equations (39) and ( 40) must be such as to cancel singularities at the boundaries.For the center, the terms of interest in Equation (35) both vary as x m M z v V z L z p P z t T For present purposes, it is convenient to proceed with foreknowledge of the converged solution so that the unknown boundary values are normalized.The factors are is the Eddington limit for the dynamical stability of a star whose opacity κ is constant, as it is in the present case, for which wellbehaved up to a common point ξ.From Equations (17)-(20), generate J solutions     and Equations (29) and (30) give improved values:


that correspond to the converged solution in order that unknown boundary values converge to:
) compare it to the std Case 1. Convergence in Case 2 is assisted by the extra rPL requirement of continuity in the product 1 2 the opacity throughout is due to Thompson scattering by free electrons, κ = 0.02(1 + X) (kg -1 ·m 2 ).Hydrogen fusion to helium is by the carbon-nitrogen-oxygen (CNO) cycle for which we let the combined C, N and O abundance be CNO 2 Z Z  Section 4.1 is met.For the surface, by Equations (81)-(83) in the limit m M L by a factor of only 8.