- Email: [email protected]
A degree-increasing [N to N+1 ] homotopy for Chebyshev and Fourier spectral methods
Accepted Manuscript A degree-increasing [N to N + 1] homotopy for Chebyshev and Fourier spectral methods John P. Boyd PII: DOI: Reference:
S0893-9659(16)30001-5 http://dx.doi.org/10.1016/j.aml.2016.01.001 AML 4918
To appear in:
Applied Mathematics Letters
Received date: 13 November 2015 Revised date: 7 January 2016 Accepted date: 7 January 2016 Please cite this article as: J.P. Boyd, A degree-increasing [N to N + 1] homotopy for Chebyshev and Fourier spectral methods, Appl. Math. Lett. (2016), http://dx.doi.org/10.1016/j.aml.2016.01.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to view linked References
A Degree-increasing [N to N + 1] Homotopy for Chebyshev and Fourier Spectral Methods John P. Boyd Department of Climate and Space Sciences and Engineering University of Michigan 2455 Hayward Avenue, Ann Arbor MI 48109-2143 [email protected]
Abstract Hitherto, as a tool for tracing all branches of nonlinear differential equations, resolution-increasing homotopy methods have been applied only to finite difference discretizations. However, spectral Galerkin algorithms typically match the error of fourth order differences with one-half to one-fifth the number of degrees of freedom N in one dimension, and a factor of eight to a hundred and twenty~N five in three dimensions. Let u ~ N be the vector of spectral coefficients and R the vector of N Galerkin constraints. A common two-part procedure is to first ~ N (~uN ) = ~0 using resultants, Groebner basis methods or block find all roots of R matrix companion matrices. (These methods are slow and ill-conditioned, practical only for small N .) The second part is to then apply resolution-increasing continuation. Because the number of solutions is an exponential function of N , spectral methods are exponentially superior to finite differences in this context. Unfortunately, ~uN is all too often outside the domain of convergence of Newton’s iteration when N is increased to (N + 1). We show that a good option ~ u; τ ) ≡ R ~ N+1 (~u) − (1 − τ )R ~ N+1 (~uN ), is the artificial parameter homotopy H(~ τ ∈ [0, 1].) Narcgubg in small steps in τ , we proceed smoothly from the N -term to the N + 1-term approximations. Keywords: Chebyshev polynomials; continuation; homotopy; pseudospectral
1. Introduction When a nonlinear differential or integral equation is discretized with N degrees of freedom (grid point values or Fourier coefficients), the result is a set of N nonlinear equations. If the nonlinearity is a polynomial in the unknown u(x), the discretization is a system of N -variate polynomials. Unfortunately, such systems are often very challenging. Recently, however, reliable solvers for polynomial systems have become widely available. Most are very slow and therefore limited to small N . As noted by Allgower, Bates, Sommese and Wampler [1], though, small-N solutions are only the “opening
Preprint submitted to Elsevier
January 7, 2016
act”. The “continuation method” is the simple but powerful strategy of varying a parameter (either physical or numerical) in small steps, using the converged solution for the previous parameter value, or the extrapolation of several previous solutions, as the first guess for the next. In principle, one can track all the physically-interesting branches as accurately as one pleases by combining continuation in the physical parameters with resolution-increasing continuation in the number of grid points or spectral coefficients N . Unfortunately, N is always an integer and there is no compelling reason for the N -term solution to lie within the convergence domain of a rootfinding iteration in the (N + 1)-dimensional parameter space. The proposed remedy is an artificial parameter homotopy. If the system ~ N and its solution by u of nonlinear algebraic equations is denoted by R ~N , a vector of Fourier or Chebyshev coefficients, then the “Newton Degree-Increasing Spectral Homotopy” [DISH] is the one-parameter family of nonlinear polynomial or transcendental equations ~ u) ≡ R ~ N+M (~u) − (1 − τ )R ~ N+M (~uN ), H(~
τ ∈ [0, 1]
(1)
where M ≥ 1 is a positive integer. We mmostly choose M = 1, but there is no inherent restriction on its size. The first step is to make a vector of length (N + M ) whose first N elements are the spectral coefficients that solve the degree N system and the remaining M elements are zeros. When τ = 0, the ~ N+M (~u) − (1 − τ )R ~ N+M (~uN ) whose exact solution is u homotopy is R ~ N . When ~ u; τ = 1) = R ~ N+M (~u). τ = 1, the homotopy reduces to our target system H(~ Thus, the homotopy system does indeed smoothly morph from the solution we know to the solution we seek. We can march from N to N + M using steps in the parameter τ as small as we please. This “Newton homotopy” was “commonly used” by 1968 [6]. Newton-DISH is novel only as an instance of the family of Degree-Increasing Spectral Homotopies. Trajectories may collide if τ is real; the collision remedy [unneeded here] is to march from τ = 0 to τ = 1 through a semicircle or other user-chosen contour in the complex plane [2]. 2. Example The “Fifth-Degree Korteweg-deVries” [FKdV] equation is −ν uXXXXX + uXXX + (u − c)uX = 0
[FKdV Eq.]
(2)
subject to the periodic boundary condition that u(X) = u(X + 2π) [5, 4]. The Fourier approximation is uN (X) = A cos(X) +
N X
n=2
2
an cos(nX)
(3)
The solution branches are parameterized by “amplitude” A, the coefficient of cos(X). The phase speed c is an unknown. The coefficient-fixed parameterization excludes the trivial solution that all coefficients are zero. The residual function is the result of substituting the Fourier series into the differential equation R(x; c, a2, . . . aN ) = −ν uN,XXXXX + uN,XXX + (u − c)uN,X
(4)
Galerkin’s method minimizes the residual through the constraints that the first N spectral coefficients of the residual are zero, Z 1 π sin(X) R(x; c, a2 , . . . aN ) dX = 0, j = 1, 2, . . . N (5) Rj = π 0 The elements of the spectral homotopy from N = 2 to N = 3 are H1
=
2A + 2 ν A + 2 cA − Aa2 − a2 a3
(6)
H2
=
64 ν a2 + 16 a2 + 4 ca2 − A2 − 2 Aa3
(7)
H3
=
{54 + 486 ν + 6c} a3 − τ 3A a2
(8)
The terms in the boxes are those of the N = 2 system. We assume that by unspecified means, we have already calculated its solution u2 (x). The spectral homotopy for N = 3 differs from the standard Galerkin discretization only through the last component of the residual vector, evaluated with the zeropadded N = 2 solution so that H1 = R1 , H2 = R2 and H3 = R3 − (1 − τ )(−3Aa2 ). The third equation can be rewritten as a3 = τ
3 Aa2 54 + 486 ν + 6c
(9)
This shows that the spectral coefficient a3 varies smoothly from zero to its value in the N = 3 solution as the parameter τ varies from zero to one. Without the homotopy, the jump from N = 2 to N = 3 is the addition of a third spectral coefficient equal to 3Aa2 /(54 + 486 ν + 6c). Note that this can be arbitrarily large; the denominator is zero whenever ν = −1/9 − c/81. For infinitesimal amplitude (|A| 1), c = −1−ν and then the denominator is zero for ν = −1/10. Haupt and Boyd show that in the small amplitude limit, cos(kX) is resonant with cos(X) whenever ν = −1/(k 2 + 1) where k is an integer. These resonances persist at finite amplitude. Obviously, resonance in the third Fourier component can make the jump from N = 2 to N = 3 very large. Fig. 1 shows the consequences of resonance; Newton’s iteration without homotopy (initialized with the exact N = 2 solution) diverges everywhere within the wedge-shaped region [shaded] in the plane spanned by the wave amplitude, A, and the coefficient of the fifth derivative in the wave equation, ν. The wedge emanates from ν = −1/10, the point where infinitesimal traveling waves cos(X) and cos(3X) have identical phase speeds. 3
The divergent region spans a larger and larger interval in ν as the wave amplitude A increases. In contrast, the Newton Degree-Increasing Spectral Homotopy [NewtonDISH] never failed. Continuation in the artificial parameter was implemented by a predictor-corrector scheme. Newton’s iteration was the corrector; we chose as predictor the fourth order Runge-Kutta method applied to the Davidenko differential equation [2] d~u = − dτ
−1 ~ dH(~u) ~~ J(~uN ) dτ
(10)
~ where J~ is the usual Jacobian matrix [same as in Newton’s iteration] and all vectors and matrices are dimension (N + M ). Fig. 2 compares homotopy results for a resonant (ν = −1/10) and a non-resonant ν. The left panel shows that the variation of the spectral coefficient with τ is curved (though still very smooth) for resonance but straight in the absence of resonance. It is possible to compute roots entirely by Runge-Kutta integration of the Davidenko equation [2]; the curves on the right show that the Newton-free integration is indeed accurate, but error is a thousand times worse with resonance.
0 10
20
-1/20
Divergence
-1/10
ν
20 30
-3/20 10 -1/5 -1/4 5
10 A
15
20
Figure 1: The contours are labeled by the number of iterations required for Newton’s iteration (with line search underrelaxation [2]) to lower the maximum residual of the N = 3 Galerkin system, initialized with the N= 2 Fourier-Galerkin solution, below 10−8. The wedge-shaped shaded region is where the iteration — without amplitude continuation or the spectral homotopy —fails. The parameters are the coefficient of cos(X), A, and the coefficient ν of the highest derivative in the FKdV equation. Without line search, the [classic] Newton’s iteration diverges in the wider wedge bounded by the dashed lines.
4
RK4 predictor/Newton corrector, A=20 −4 14 10 12
ν=−0.1 ν=0
RK/Newt − RK
−5
10
ν=−0.1 ν=0
10 −6
a3
10 a3
8 6
−7
10
4 −8
10
2 0 0
−9
0.5
τ
10
1
0
0.5
τ
1
Figure 2: Left: Variation of Fourier coefficient a3 with τ for traveling waves of the FKdV equation when the two-term solution for A = 20 is continued to an N = 8 Fourier cosine Galerkin approximation for two different values of the ODE parameter ν. The right panel shows the difference between the predictor-corrector results [accurate to near-machine precision] and the Runge-Kutta integration of the Davidenko equation without iterative correction. The upper curve is for the resonance at ν = −1/10. The errors of pure Runge-Kutta continuation are small compared to the solution u(x), but large compared to the errors when ν is not close to resonance [lower curve].
3. Theory Our full-length follow-up will contain many additional examples [3]. It will also justify the following four useful assertions. 1. Newton-DISH is “regular” in the sense that it will succeed whenever the direct Newton’s iteration will succeed; there is thus no penalty for applying it even when it is not strictly necessary. 2. If the approximations to u(x) with N basis functions, uN (x), is sufficiently accurate, then a resolution-increasing homotopy is probably unnecessary; the raison d’etre of the spectral homotopy is to continue to higher resolution from an approximation whose accuracy is only moderate. 3. uN (x) is a poor approximation to uN+1 (x) when the (N +1)-th basis function is resonant; this usually wrecks the tactic of uN (x) as the initialization for the iteration for larger N . 4. The artificial parameter homotopy, Newton-DISH, converts one big change in the coefficient of the (N + 1)-th basis function into the end result of many small changes and is usually successful. 4. Summary Second order finite differences have been used in the “all-branches” attack on differential equations, but spectral methods are exponentially superior [3]. 5
This note makes it possible to pair an all-roots polynomial rootsolver with an equally advanced differential equation-solver by explicitly constructing a DISH: Degree-Increasing Spectral Homotopy. The simplest approach to resolution-increasing continuation is to use the N -term approximation, padded with M zero Fourier or Chebyshev polynomial coefficients, as the initial approximation for computing a higher degree N + M term approximation by Newton’s iteration. We have shown by example that when the degree is low, the jump from N basis functions to N +1 basis functions may be too much. In contrast, the artificial parameter homotopy introduced here splits the distance from one resolution to another into small steps. Resolution-increasing continuation is no longer restricted to the ingenious scheme for finite differences devised by Allgower, Bates, Sommese and Wampler [1]. Fourier and Chebyshev Galerkin discretizations, much more accurate than finite differences can be employed to convert a differential equation into a small system of polynomial equations which can be solved completely for all its solution branches. The homotopy described here extends these small-N solutions to larger degree to trace all the physically relevant solutions with spectral accuracy. Acknowledgments: NSF grants OCE 1059703 and DMS 1521158. I thank Charles Wampler, Jonathan Hauenstein and two reviewers for their help. References [1] E. L. Allgower, D. J. Bates, A. Sommese, and C. W. Wampler, Solution of polynomial systems derived from differential equations, Computing, 76 (2006), pp. 1–10. [2] J. P. Boyd, Solving Transcendental Equations: The Chebyshev Polynomial Proxy and Other Numerical Rootfinders, Perturbation Series and Oracles, SIAM, Philadelphia, 2014. 460 pp. [3]
, The advantages of spectral methods for all-branch root-tracing for differential and integral equations and an N -modes-to-(N +1)-modes homotopy, J. Comput. Phys., (2015). submitted.
[4] J. P. Boyd and S. E. Haupt, Polycnoidal waves: Spatially periodic generalizations of multiple solitary waves, in Nonlinear Topics of Ocean Physics: Fermi Summer School, Course LIX, A. R. Osborne, ed., North-Holland, Amsterdam, 1991, pp. 827–856. [5] S. E. Haupt and J. P. Boyd, Modeling nonlinear resonance: A modification to Stokes’ perturbation expansion, Wave Motion, 10 (1988), pp. 83–98. [6] G. H. Meyer, On solving nonlinear equations with a one-parameter operator imbedding, SIAM J. Numer. Anal., 5 (1968), pp. 739–740.
6
S0893-9659(16)30001-5 http://dx.doi.org/10.1016/j.aml.2016.01.001 AML 4918
To appear in:
Applied Mathematics Letters
Received date: 13 November 2015 Revised date: 7 January 2016 Accepted date: 7 January 2016 Please cite this article as: J.P. Boyd, A degree-increasing [N to N + 1] homotopy for Chebyshev and Fourier spectral methods, Appl. Math. Lett. (2016), http://dx.doi.org/10.1016/j.aml.2016.01.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to view linked References
A Degree-increasing [N to N + 1] Homotopy for Chebyshev and Fourier Spectral Methods John P. Boyd Department of Climate and Space Sciences and Engineering University of Michigan 2455 Hayward Avenue, Ann Arbor MI 48109-2143 [email protected]
Abstract Hitherto, as a tool for tracing all branches of nonlinear differential equations, resolution-increasing homotopy methods have been applied only to finite difference discretizations. However, spectral Galerkin algorithms typically match the error of fourth order differences with one-half to one-fifth the number of degrees of freedom N in one dimension, and a factor of eight to a hundred and twenty~N five in three dimensions. Let u ~ N be the vector of spectral coefficients and R the vector of N Galerkin constraints. A common two-part procedure is to first ~ N (~uN ) = ~0 using resultants, Groebner basis methods or block find all roots of R matrix companion matrices. (These methods are slow and ill-conditioned, practical only for small N .) The second part is to then apply resolution-increasing continuation. Because the number of solutions is an exponential function of N , spectral methods are exponentially superior to finite differences in this context. Unfortunately, ~uN is all too often outside the domain of convergence of Newton’s iteration when N is increased to (N + 1). We show that a good option ~ u; τ ) ≡ R ~ N+1 (~u) − (1 − τ )R ~ N+1 (~uN ), is the artificial parameter homotopy H(~ τ ∈ [0, 1].) Narcgubg in small steps in τ , we proceed smoothly from the N -term to the N + 1-term approximations. Keywords: Chebyshev polynomials; continuation; homotopy; pseudospectral
1. Introduction When a nonlinear differential or integral equation is discretized with N degrees of freedom (grid point values or Fourier coefficients), the result is a set of N nonlinear equations. If the nonlinearity is a polynomial in the unknown u(x), the discretization is a system of N -variate polynomials. Unfortunately, such systems are often very challenging. Recently, however, reliable solvers for polynomial systems have become widely available. Most are very slow and therefore limited to small N . As noted by Allgower, Bates, Sommese and Wampler [1], though, small-N solutions are only the “opening
Preprint submitted to Elsevier
January 7, 2016
act”. The “continuation method” is the simple but powerful strategy of varying a parameter (either physical or numerical) in small steps, using the converged solution for the previous parameter value, or the extrapolation of several previous solutions, as the first guess for the next. In principle, one can track all the physically-interesting branches as accurately as one pleases by combining continuation in the physical parameters with resolution-increasing continuation in the number of grid points or spectral coefficients N . Unfortunately, N is always an integer and there is no compelling reason for the N -term solution to lie within the convergence domain of a rootfinding iteration in the (N + 1)-dimensional parameter space. The proposed remedy is an artificial parameter homotopy. If the system ~ N and its solution by u of nonlinear algebraic equations is denoted by R ~N , a vector of Fourier or Chebyshev coefficients, then the “Newton Degree-Increasing Spectral Homotopy” [DISH] is the one-parameter family of nonlinear polynomial or transcendental equations ~ u) ≡ R ~ N+M (~u) − (1 − τ )R ~ N+M (~uN ), H(~
τ ∈ [0, 1]
(1)
where M ≥ 1 is a positive integer. We mmostly choose M = 1, but there is no inherent restriction on its size. The first step is to make a vector of length (N + M ) whose first N elements are the spectral coefficients that solve the degree N system and the remaining M elements are zeros. When τ = 0, the ~ N+M (~u) − (1 − τ )R ~ N+M (~uN ) whose exact solution is u homotopy is R ~ N . When ~ u; τ = 1) = R ~ N+M (~u). τ = 1, the homotopy reduces to our target system H(~ Thus, the homotopy system does indeed smoothly morph from the solution we know to the solution we seek. We can march from N to N + M using steps in the parameter τ as small as we please. This “Newton homotopy” was “commonly used” by 1968 [6]. Newton-DISH is novel only as an instance of the family of Degree-Increasing Spectral Homotopies. Trajectories may collide if τ is real; the collision remedy [unneeded here] is to march from τ = 0 to τ = 1 through a semicircle or other user-chosen contour in the complex plane [2]. 2. Example The “Fifth-Degree Korteweg-deVries” [FKdV] equation is −ν uXXXXX + uXXX + (u − c)uX = 0
[FKdV Eq.]
(2)
subject to the periodic boundary condition that u(X) = u(X + 2π) [5, 4]. The Fourier approximation is uN (X) = A cos(X) +
N X
n=2
2
an cos(nX)
(3)
The solution branches are parameterized by “amplitude” A, the coefficient of cos(X). The phase speed c is an unknown. The coefficient-fixed parameterization excludes the trivial solution that all coefficients are zero. The residual function is the result of substituting the Fourier series into the differential equation R(x; c, a2, . . . aN ) = −ν uN,XXXXX + uN,XXX + (u − c)uN,X
(4)
Galerkin’s method minimizes the residual through the constraints that the first N spectral coefficients of the residual are zero, Z 1 π sin(X) R(x; c, a2 , . . . aN ) dX = 0, j = 1, 2, . . . N (5) Rj = π 0 The elements of the spectral homotopy from N = 2 to N = 3 are H1
=
2A + 2 ν A + 2 cA − Aa2 − a2 a3
(6)
H2
=
64 ν a2 + 16 a2 + 4 ca2 − A2 − 2 Aa3
(7)
H3
=
{54 + 486 ν + 6c} a3 − τ 3A a2
(8)
The terms in the boxes are those of the N = 2 system. We assume that by unspecified means, we have already calculated its solution u2 (x). The spectral homotopy for N = 3 differs from the standard Galerkin discretization only through the last component of the residual vector, evaluated with the zeropadded N = 2 solution so that H1 = R1 , H2 = R2 and H3 = R3 − (1 − τ )(−3Aa2 ). The third equation can be rewritten as a3 = τ
3 Aa2 54 + 486 ν + 6c
(9)
This shows that the spectral coefficient a3 varies smoothly from zero to its value in the N = 3 solution as the parameter τ varies from zero to one. Without the homotopy, the jump from N = 2 to N = 3 is the addition of a third spectral coefficient equal to 3Aa2 /(54 + 486 ν + 6c). Note that this can be arbitrarily large; the denominator is zero whenever ν = −1/9 − c/81. For infinitesimal amplitude (|A| 1), c = −1−ν and then the denominator is zero for ν = −1/10. Haupt and Boyd show that in the small amplitude limit, cos(kX) is resonant with cos(X) whenever ν = −1/(k 2 + 1) where k is an integer. These resonances persist at finite amplitude. Obviously, resonance in the third Fourier component can make the jump from N = 2 to N = 3 very large. Fig. 1 shows the consequences of resonance; Newton’s iteration without homotopy (initialized with the exact N = 2 solution) diverges everywhere within the wedge-shaped region [shaded] in the plane spanned by the wave amplitude, A, and the coefficient of the fifth derivative in the wave equation, ν. The wedge emanates from ν = −1/10, the point where infinitesimal traveling waves cos(X) and cos(3X) have identical phase speeds. 3
The divergent region spans a larger and larger interval in ν as the wave amplitude A increases. In contrast, the Newton Degree-Increasing Spectral Homotopy [NewtonDISH] never failed. Continuation in the artificial parameter was implemented by a predictor-corrector scheme. Newton’s iteration was the corrector; we chose as predictor the fourth order Runge-Kutta method applied to the Davidenko differential equation [2] d~u = − dτ
−1 ~ dH(~u) ~~ J(~uN ) dτ
(10)
~ where J~ is the usual Jacobian matrix [same as in Newton’s iteration] and all vectors and matrices are dimension (N + M ). Fig. 2 compares homotopy results for a resonant (ν = −1/10) and a non-resonant ν. The left panel shows that the variation of the spectral coefficient with τ is curved (though still very smooth) for resonance but straight in the absence of resonance. It is possible to compute roots entirely by Runge-Kutta integration of the Davidenko equation [2]; the curves on the right show that the Newton-free integration is indeed accurate, but error is a thousand times worse with resonance.
0 10
20
-1/20
Divergence
-1/10
ν
20 30
-3/20 10 -1/5 -1/4 5
10 A
15
20
Figure 1: The contours are labeled by the number of iterations required for Newton’s iteration (with line search underrelaxation [2]) to lower the maximum residual of the N = 3 Galerkin system, initialized with the N= 2 Fourier-Galerkin solution, below 10−8. The wedge-shaped shaded region is where the iteration — without amplitude continuation or the spectral homotopy —fails. The parameters are the coefficient of cos(X), A, and the coefficient ν of the highest derivative in the FKdV equation. Without line search, the [classic] Newton’s iteration diverges in the wider wedge bounded by the dashed lines.
4
RK4 predictor/Newton corrector, A=20 −4 14 10 12
ν=−0.1 ν=0
RK/Newt − RK
−5
10
ν=−0.1 ν=0
10 −6
a3
10 a3
8 6
−7
10
4 −8
10
2 0 0
−9
0.5
τ
10
1
0
0.5
τ
1
Figure 2: Left: Variation of Fourier coefficient a3 with τ for traveling waves of the FKdV equation when the two-term solution for A = 20 is continued to an N = 8 Fourier cosine Galerkin approximation for two different values of the ODE parameter ν. The right panel shows the difference between the predictor-corrector results [accurate to near-machine precision] and the Runge-Kutta integration of the Davidenko equation without iterative correction. The upper curve is for the resonance at ν = −1/10. The errors of pure Runge-Kutta continuation are small compared to the solution u(x), but large compared to the errors when ν is not close to resonance [lower curve].
3. Theory Our full-length follow-up will contain many additional examples [3]. It will also justify the following four useful assertions. 1. Newton-DISH is “regular” in the sense that it will succeed whenever the direct Newton’s iteration will succeed; there is thus no penalty for applying it even when it is not strictly necessary. 2. If the approximations to u(x) with N basis functions, uN (x), is sufficiently accurate, then a resolution-increasing homotopy is probably unnecessary; the raison d’etre of the spectral homotopy is to continue to higher resolution from an approximation whose accuracy is only moderate. 3. uN (x) is a poor approximation to uN+1 (x) when the (N +1)-th basis function is resonant; this usually wrecks the tactic of uN (x) as the initialization for the iteration for larger N . 4. The artificial parameter homotopy, Newton-DISH, converts one big change in the coefficient of the (N + 1)-th basis function into the end result of many small changes and is usually successful. 4. Summary Second order finite differences have been used in the “all-branches” attack on differential equations, but spectral methods are exponentially superior [3]. 5
This note makes it possible to pair an all-roots polynomial rootsolver with an equally advanced differential equation-solver by explicitly constructing a DISH: Degree-Increasing Spectral Homotopy. The simplest approach to resolution-increasing continuation is to use the N -term approximation, padded with M zero Fourier or Chebyshev polynomial coefficients, as the initial approximation for computing a higher degree N + M term approximation by Newton’s iteration. We have shown by example that when the degree is low, the jump from N basis functions to N +1 basis functions may be too much. In contrast, the artificial parameter homotopy introduced here splits the distance from one resolution to another into small steps. Resolution-increasing continuation is no longer restricted to the ingenious scheme for finite differences devised by Allgower, Bates, Sommese and Wampler [1]. Fourier and Chebyshev Galerkin discretizations, much more accurate than finite differences can be employed to convert a differential equation into a small system of polynomial equations which can be solved completely for all its solution branches. The homotopy described here extends these small-N solutions to larger degree to trace all the physically relevant solutions with spectral accuracy. Acknowledgments: NSF grants OCE 1059703 and DMS 1521158. I thank Charles Wampler, Jonathan Hauenstein and two reviewers for their help. References [1] E. L. Allgower, D. J. Bates, A. Sommese, and C. W. Wampler, Solution of polynomial systems derived from differential equations, Computing, 76 (2006), pp. 1–10. [2] J. P. Boyd, Solving Transcendental Equations: The Chebyshev Polynomial Proxy and Other Numerical Rootfinders, Perturbation Series and Oracles, SIAM, Philadelphia, 2014. 460 pp. [3]
, The advantages of spectral methods for all-branch root-tracing for differential and integral equations and an N -modes-to-(N +1)-modes homotopy, J. Comput. Phys., (2015). submitted.
[4] J. P. Boyd and S. E. Haupt, Polycnoidal waves: Spatially periodic generalizations of multiple solitary waves, in Nonlinear Topics of Ocean Physics: Fermi Summer School, Course LIX, A. R. Osborne, ed., North-Holland, Amsterdam, 1991, pp. 827–856. [5] S. E. Haupt and J. P. Boyd, Modeling nonlinear resonance: A modification to Stokes’ perturbation expansion, Wave Motion, 10 (1988), pp. 83–98. [6] G. H. Meyer, On solving nonlinear equations with a one-parameter operator imbedding, SIAM J. Numer. Anal., 5 (1968), pp. 739–740.
6