Computing eigenfunctions of singular points in nonlinear parametrized two-point BVPs

Applied Numerical Mathematics 59 (2009) 671–676 www.elsevier.com/locate/apnum

Computing eigenfunctions of singular points in nonlinear parametrized two-point BVPs M. Hermann ∗ , Th. Milde Friedrich Schiller University, Institute of Applied Mathematics, Ernst-Abbe-Platz 2, D-07743 Jena, Germany Available online 20 March 2008

Abstract The iterative computation of singular points in parametrized nonlinear BVPs by so-called extended systems requires good starting values for the singular point itself and the corresponding eigenfunction. Using path-following techniques such starting values for the singular points should be generated automatically. However, path-following does not provide approximations for the eigenfunction if the singularity is a bifurcation point. We propose a new modification of this standard technique delivering such starting values. It is based on an extended system by which singular as well as nonsingular points can be determined. © 2008 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Boundary value problem; Bifurcation points; Turning points; Path-following; Augmented systems

1. Introduction Consider the parametrized two-point boundary value problem (BVP) for a system of n nonlinear ordinary differential equations   x  (t) − f t, x(t); λ = 0, a < t < b, (1) subject to n nonlinear boundary conditions   r x(a), x(b); λ = 0,

(2)

where x(t) ∈ Rn , f : Ω1 → Rn , Ω1 ⊂ (a, b) × Rn × R, r : Ω2 → Rn , Ω2 ⊂ Rn × Rn × R, and λ ∈ R is a control parameter. Assume that f and r are sufficiently smooth. In order to generate parts of the solution manifold of (1), (2) and to draw the associated bifurcation diagram, three problems have to be considered (see, e.g. [4,7,8]): (1) tracing the solution curves numerically, (2) recognizing and determining singular points (bifurcation and/or turning points), and (3) branch switching at the bifurcation points. * Corresponding author.

E-mail addresses: [email protected] (M. Hermann), [email protected] (Th. Milde). 0168-9274/$30.00 © 2008 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2008.03.010

672

M. Hermann, Th. Milde / Applied Numerical Mathematics 59 (2009) 671–676

In this paper we will deal with the second problem, namely the determination of singular points by so-called augmented or extended systems. Let (x, λ) be a solution of (1), (2). We define the matrices     A t, x(t), λ ≡ fx t, x(t); λ ,   Ba (x, λ) ≡ rx(a) x(a), x(b); λ ,   (3) Bb (x, λ) ≡ rx(b) x(a), x(b); λ . A solution q ≡ (x (1) , λ1 ) is called a singular point of (1), (2) if the corresponding variational problem ϕ  (t) − A(t, q)ϕ(t) = 0,

Ba (q)ϕ(a) + Bb (q)ϕ(b) = 0

(4)

has a nontrivial solution ϕ (1) (the eigenfunction of the BVP (1), (2)). Hence, a nonsingular point is a solution q which is not singular. Here, we consider only simple singular points, i.e. singular points for which (4) has a one-dimensional null space (simple turning and bifurcation points). Nonsimple singular points will be studied in a forthcoming paper. Let us now assume that a nonsingular point p ≡ (x (0) , λ0 ) of (1), (2) is given. This means that a point on a (a priori unknown) solution curve of the BVP is known. Starting from that point, pseudo-arc-length continuation (see, e.g. [4,5,7]), which is a kind of arc-length continuation, is usually used to compute further nonsingular points on this solution curve. The continuation method is also referred to as path-following technique. In our code RWPKV [4] the path-following is realized on the basis of Seydel’s predictor-corrector algorithm [6] and the multiple shooting algorithm RWPM (see, e.g. [2]). Observing a test-function (e.g. the determinant of the simple-shooting matrix) avoids jumping over a singular point, say q ≡ (x (1) , λ1 ). Let us here assume that there exists a so-called simple solution curve (consisting of nonsingular points only) between p and q and that q is a simple singular point. Usually so-called extended (augmented) systems are used to determine numerically the singular point q (see, e.g. [7–9]). These systems contain as an essential part the BVP (4). If standard solvers like the multiple shooting code RWPM [3] are applied to such an extended system, then good startingvalues for all components of the solution are required. Obviously, the path-following method provides approximations for x (1) and λ1 . In case of a simple turning point an approximation of ϕ (1) can be computed without extra costs by the ˆ − x(t; λ) ¯ (see, e.g. [7], p. 112). However, if the singularity is a bifurcation scaled difference of two curve points x(t; λ) point the situation is much more complicated. In this paper we propose a general numerical strategy to determine an approximation of ϕ (1) for all singular points with a one-dimensional null space of the variational problem (4). It is based on a special inhomogeneity v which is inserted into the differential equations (4). 2. A new extended system In the previous section we assumed, that the path between p = (x (0) , λ0 ) and q = (x (1) , λ1 ) is a simple solution curve. Moreover, we suppose that p belongs to the path and the simple singular point q lies outside. It follows from the Implicit Function Theorem, that the simple solution curve, consisting of solutions {(x(t), λ), λ ∈ [λ0 , λ1 )}, can be parametrized by the parameter λ. Thus we can write x = x(t, λ) and the matrices (3) are of the form   (5) A(t, λ) ≡ A t, x(t; λ), λ , Ba (λ) ≡ Ba (x, λ), Bb (λ) ≡ Bb (x, λ). Let us now consider the following modification of the linear BVP (4): ϕ  (t, λ) − A(t, λ)ϕ(t, λ) = v(t, λ),

Ba (λ)ϕ(a, λ) + Bb (λ)ϕ(b, λ) = 0.

(6)

Under our assumptions the BVP (6) has a unique solution for every given v. The problem is now to determine a function v such that v(t, λ) −−−−→ 0 λ↑λ1

(7)

and the corresponding solutions ϕ(t, λ) tend to an approximation of ϕ (1) . Let us assume that such a v is known. Then, (6) can be solved by the simple shooting technique for linear BVPs which is based on the following steps:

M. Hermann, Th. Milde / Applied Numerical Mathematics 59 (2009) 671–676

673

(1) Compute a fundamental matrix X(t, λ) which is defined by the homogeneous matrix-IVP X  (t, λ) = A(t, λ) X(t, λ),

X(a, λ) = I.

(2) Compute a particular solution g(t, λ) which is defined by the IVP g  (t, λ) − A(t, λ) g(t, λ) = v(t, λ),

g(a, λ) = 0.

(3) Determine the shooting-matrix M(λ) ≡ Ba (λ) + Bb (λ)X(b, λ) and compute ϕ(a, λ) ≡ c(λ) as the uniquely determined solution of the n-dimensional linear algebraic system M(λ) c(λ) = −Bb (λ) g(b, λ) ≡ d(λ).

(8)

Obviously, the vector d(λ) and the function v(t, λ) are interrelated. In a first step we will demonstrate how d(λ) has to be chosen that d(λ) −−−−→ 0

(9)

λ↑λ1

and the solutions c(λ) of (8) converge to a vector c = 0 which satisfies M(λ1 ) c = 0. Later on we will show how a function v satisfying (7) can be determined on the basis of d. The choice of d follows from the following theorem. Theorem 1. Let the continuous matrix function M(λ) ∈ Rn×n be invertible for λ ∈ [λ0 , λ1 ). Furthermore, assume that the null space N (M(λ1 )) is one-dimensional and that the vector s ∈ Rn does not belong to R(M(λ1 )). If c(λ) denotes the solution of the linear algebraic system   M(λ) c(λ) = det M(λ) s, (10) then the limit limλ↑λ1 c(λ) ≡ c = 0 exists and satisfies M(λ1 ) c = 0. Proof. Let M(λ) = [m1 (λ)|m2 (λ)| · · · |mn (λ)] be the column partition of the matrix M(λ) ∈ Rn×n . The j -th entry of the vector c(λ) ∈ Rn is determined by Cramer’s rule as det([m1 (λ), . . . , mj −1 (λ), det(M(λ))s, mj +1 (λ), . . . , mn (λ)]) det(M(λ))   = det m1 (λ), . . . , mj −1 (λ), s, mj +1 (λ), . . . , mn (λ) .

cj (λ) =

Since M and the determinant are continuous functions of λ the following limit exists: cj (λ) −−−−→ cj . λ↑λ1

The assumptions imply that cj = 0 for at least one j ∈ {1, . . . , n}.

2

As we have seen, the conclusion above is based on the fact that the vector s ∈ Rn is not an element of R(M(λ1 )). It is well known that a vector s¯ belongs to R(M(λ1 )) if and only if s¯ is orthogonal to the null space N (M(λ1 )T ). The orthogonal complement of this null space is a one-dimensional subspace of Rn . Therefore, if the components of s are chosen as uniformly distributed numbers in the interval [−1, 1], the probability of s belonging to a specific one-dimensional subspace is zero. Thus, an arbitrary s will satisfy the assumptions of Theorem 1 with the utmost probability. Let us now return to the question how a function v(t, λ) can be found which satisfies the relation (7). From the theory of linear BVPs we know that a function v (1) (t, λ) exists for which the BVP (6), with λ = λ1 and v = v (1) as right hand side, has no solution. Such an inhomogeneity v (1) leads to a vector s = −Bb (λ1 )g(b, λ1 ), with M(λ1 ) c = s is not solvable, i.e. s ∈ / R(M(λ1 )). Obviously, the choice v(t, λ) ≡ det(M(λ))v (1) (t, λ) leads to     d(λ) ≡ − det M(λ) Bb (λ1 )g(b, λ1 ) = det M(λ) s. From Theorem (2) and Gronwall’s Lemma it now follows ϕ(a, λ) ≡ c(λ) −−−−→ c ≡ ϕ(a, λ1 ) = 0. λ↑λ1

Computing ϕ(t, λ1 ) from the ODE ϕ  (t, λ1 ) − A(t, λ1 ) ϕ(t, λ1 ) = 0 and the known initial vector ϕ(a, λ1 ) yields the required approximation for ϕ (1) .

674

M. Hermann, Th. Milde / Applied Numerical Mathematics 59 (2009) 671–676

3. The algorithm The approximation of ϕ (1) described in Section 2 can be realized by the following algorithm. Here, instead of the simple shooting method the more stable multiple shooting method is used. The corresponding multiple shooting (m) matrix Mk (λ) of the k-th iteration step has the same determinant as the matrix Mk (λ) of the simple shooting method (m) (see [2]). Therefore, the shooting matrix Mfinal (λ) computed with the code RWPM in the last iteration step is again denoted by M(λ). Algorithm. 1. Choose a fixed value λ0 of the control parameter λ and compute with the multiple shooting code RWPM a starting solution x (0) (t) of the nonlinear BVP (1), (2) at m + 1 shooting points a = τ0 < τ1 < · · · < τm = b. 2. Determine the matrices Ba (λ0 ) ≡ Ba (x (0) , λ0 ), Bb (λ0 ) ≡ Bb (x (0) , λ0 ) and det(M(λ0 )). 3. Choose an arbitrary nontrivial function u(t) satisfying the homogeneous equation Ba (λ0 )u(a) + Bb (λ0 )u(b) = 0. For instances, this can be realized by choosing a vector v = 0 and determining u(t) as a quadratic polynomial in each component, i.e. satisfying the conditions u(a) = 0, u(b) = 0 and u((b − a)/2) = v. 4. Determine the function xˆ (0) (t) as the piecewise linear polynomial interpolating the vectors x (0) (τ0 ), x (0) (τ1 ), . . . , x (0) (τm ) which have been computed in Step 1. Set   h(t) ≡ u (t) − fx t, xˆ (0) (t); λ0 u(t). 5. Compute the solutions (x, λ, ϕ) of the extended system   r x(a), x(b), λ = 0 x  (t) − f (t, x, λ) = 0,   ˆ h(t), Ba ϕ(a) + Bb ϕ(b) = 0 ϕ  (t) − fx (t, x, λ) ϕ(t) = det M(λ)

(11)

by a path-following technique (e.g. with the code RWPKV [4]), starting with the known solution (x (0) (t), λ0 , ϕ (0) (t)) = (p, det(M(λ0 )) u(t)). Here, λˆ denotes the value of λ at the previous continuation step. However, in the first and second step λˆ should be set to λ0 . ¯ the corre6. If an appropriate test function identifies a simple singular point in the immediate vicinity of λ = λ, ¯ sponding functions ϕ(t), ¯ x(t) ¯ and λ can be used as starting approximations for the computation of the singular point q = (x (1) (t), λ1 ) by a suitable extended system. ˆ Remark 2. Since det(M(λ)) → 0 for λ → λ1 , the same is true for det(M(λ)). Remark 3. The algorithm can also be used in the case of simple turning points as it is shown in the next section. But the computational costs are more expensive than using the simple strategy mentioned in Section 1. 4. Numerical computations We report in this section on the results of three numerical experiments we have performed to illustrate some of the results given here. As a first example we solved the following parametrized BVP   x2 (t) = λ x1 (t) + x1 (t)2 , x1 (t) = x2 (t), x1 (0) = 0,

x1 (1) = 0.

(12)

Obviously, x1 (t) ≡ 0 and x2 (t) ≡ 0 is a solution of (12) for all λ ∈ R. It is well-known (see, e.g. [8]) that this parametrized BVP has simple primary bifurcation points at λ(k) = −(kπ)2 , k = 1, 2, . . . . The corresponding eigenfunctions are ϕ (k) (t) = C (sin(kπt), kπ cos(kπt))T , k = 1, 2, . . . . In Step 3 of our algorithm we used the vectorfunction u(t) = (t 2 − t, 2t)T . Then, the function h(t) in Step 4 is h(t) = (−1, t 2 − t + 2)T . We have set λ = λ0 = −1

M. Hermann, Th. Milde / Applied Numerical Mathematics 59 (2009) 671–676 (0)

675

(0)

and x1 (t) ≡ 0, x2 (t) ≡ 0. It follows det(M(λ0 )) = 1.17522. For the detection of singular points we have used the test function det(M(λ)). During the execution of Step 5 the parameter value λ = λ¯ = −9.86964 has been determined where the absolute value of the test function is 5.14494 × 10−5 . The small value indicated us the existence of a sin¯ gular point in the neighborhood of λ¯ (the known bifurcation point with λ = λ(1) ). The computed approximation ϕ(t) already agrees with the exact eigenfunction ϕ (1) (t) up to 6 significant digits, where C = −2.91813. To improve the result we have solved the extended system for primary simple bifurcation points [8] ϕ  (t) = fx (t, 0; λ)ϕ(t), 

ξ (t) = ϕ(t) ϕ(t), T

ϕ1 (0) = 0, ϕ1 (1) = 0, ξ(0) = 0, ξ(1) = 1,



λ (t) = 0

(13)

with the shooting code RWPM. Here, we used 20 equidistributed shooting points τj ∈ [0, 1), j = 1(1)20. The starting ¯ j ), ξ(τj ) = τj and λ = λ¯ . Only 5 iteration steps were necessary to obtain an approximation trajectory was ϕ(τj ) = ϕ(τ of the primary simple bifurcation point q = (0, λ1 ) with an accuracy of 12 significant digits. As a second example we studied the BVP x1 (t) = x2 (t),      2 x2 (t) = λ 2 + x1 (t) + x1 (t)2 + λ2 t 1 − t + 2x1 (t)(1 − t) + λ3 t 2 − t , x1 (0) = 0,

x1 (1) = 0.

(14)

It is known that a simple solution curve begins at the point λ = 0, x1 (t) ≡ 0, x2 (t) ≡ 0 and reaches a secondary bifurcation point q near by λ = −9.8696. In Step 3 of our algorithm we have chosen the function u(t) = (t 2 − t, t 2 )T . As a test function we used det(M(λ)). During the execution of Step 5 of our algorithm the parameter value λ = λ¯ = −9.8681 has been determined where the absolute value of the test function is 4.53418 × 10−5 . The small value indicated us that a singular point must be in the neighborhood of λ¯ . Therefore, we used λ¯ and the corresponding ˆ which is the difference quotient of functions x(t), ¯ ϕ(t) ¯ computed by the algorithm, s¯ (t) = (x(t) ¯ − x(t))/( ˆ λ¯ − λ), the last two approximations for x(t) in our algorithm, as well as ξ¯1 (t) = t, ξ¯2 (t) = t − t 2 and μ¯ ≡ 0 as starting approximations for the iterative solution of the extended system for secondary simple bifurcation points [8]   x  (t) = f t, x(t); λ − μψ0 (t), x1 (0) = 0, x1 (1) = 0,    ϕ (t) = fx t, x(t); λ ϕ(t), ϕ1 (0) = 0, ϕ1 (1) = 0,      s (t) = fx t, x(t); λ s(t) + fλ t, x(t); λ , s1 (0) = 0, s1 (1) = 0, ξ1 (t) = ϕ(t)T ϕ(t),

ξ1 (0) = 0,

ξ1 (1) = 1,

ξ2 (t) = ϕ(t)T s(t), 

ξ2 (0) = 0,

ξ2 (1) = 0,

λ (t) = 0,

μ (t) = 0.

(15)

An appropriate choice for the function ψ0 (t) was the function ϕ(t) with interchanged components. With these values and 20 equidistributed shooting points the multiple shooting code RWPM required 6 iteration steps to compute a solution of (15) with an accuracy of 10 significant digits. The corresponding value of μ was μ = 1.68 × 10−11 , i.e. the computed solution is really the secondary bifurcation point. The third test problem was the Bratu–Gelfand equation for the stationary temperature distribution in a beam   x1 (t) = x2 (t), x2 (t) = −λ exp x1 (t) , x1 (0) = 0,

x1 (1) = 0.

(16)

It is well known (see, e.g. [1,8]) that a simple solution curve begins at the point λ = 0, x1 (t) ≡ 0, x2 (t) ≡ 0 and reaches a simple turning point q near by λ = 3.51382988. In Step 3 of our algorithm we have chosen the vector-function u(t) = (t 2 − t, 2t)T . As a test function we used again det(M(λ)). During the execution of Step 5 of our algorithm the value λ = λ¯ = 3.51383 has been determined where the absolute value of the test function is 3.75348 × 10−5 . The small ¯ Therefore, we used λ¯ and the corresponding value indicated us that a turning point must be in the neighborhood of λ.

676

M. Hermann, Th. Milde / Applied Numerical Mathematics 59 (2009) 671–676

functions x(t), ¯ ϕ(t) ¯ computed by the algorithm and ξ¯ (t) = t as starting approximations for the iterative solution of the extended system for simple turning points [8]   x1 (0) = 0, x1 (1) = 0, x  (t) = f t, x(t); λ ,    ϕ (t) = fx t, x(t); λ ϕ(t), ϕ1 (0) = 0, ϕ1 (1) = 0, ξ  (t) = ϕ(t)T ϕ(t), 

λ (t) = 0.

ξ(0) = 0,

ξ(1) = 1, (17)

With these values and 20 equidistributed shooting points the multiple shooting code RWPM required only 4 iteration steps to compute a solution of (17), i.e. the turning point, with an accuracy of 12 significant digits. References [1] [2] [3] [4] [5] [6] [7] [8] [9]

W.J.F. Govaerts, Numerical Methods for Bifurcations of Dynamical Equilibria, SIAM, Philadelphia, 2000. M. Hermann, Numerik gewöhnlicher Differentialgleichungen. Anfangs- und Randwertprobleme, Oldenbourg Verlag, München, 2004. M. Hermann, D. Kaiser, Shooting methods for two-point BVPs with partially separated endconditions, ZAMM 75 (1995) 651–668. M. Hermann, K. Ullrich, RWPKV: a software package for continuation and bifurcation problems in two-point boundary value problems, Appl. Math. Lett. 5 (1992) 57–62. H.B. Keller, Lectures on Numerical Methods in Bifurcation Problems, Springer-Verlag, Heidelberg, 1987. R. Seydel, A continuation algorithm with step control, in: T. Küpper, H.D. Mittelmann, H. Weber (Eds.), Numerical Methods for Bifurcation Problems, Birkhäuser Verlag, Basel, 1984, pp. 480–494. R. Seydel, From Equilibrium to Chaos (Practical Bifurcation and Stability Analysis), Elsevier, New York, 1988. W. Wallisch, M. Hermann, Numerische Behandlung von Fortsetzungs- und Bifurkationsproblemen bei Randwertaufgaben, Teubner-Texte zur Mathematik, Band 102, Teubner Verlag, Leipzig, 1987. H. Weber, Numerische Behandlung von Verzweigungsproblemen bei Randwertaufgaben gewöhnlicher Differentialgleichungen, PhD thesis, Johannes Gutenberg-Universität Mainz, 1978.