Schwerpunktprogramm der Deutschen Forschungsgemeinschaft

Ergodentheorie, Analysis und effiziente Simulation dynamischer Systeme

Numerical Taylor expansions of invariant manifolds in large dynamical systems

W.-J. Beyn and W. Kleß

Preprint 6/96

Fakultät für Mathematik

Universität Bielefeld Postfach 100131 D-33501 Bielefeld Germany

version from 6 Aug. 1997 submitted to ”Numerische Mathematik”

Numerical Taylor expansions of invariant manifolds in large dynamical systems

Wolf–Jürgen Beyn^{∗}

email: beyn@mathematik.uni-bielefeld.de

Winfried Kleß^{∗}

email: kless@mathematik.uni-bielefeld.de

Fakultät für Mathematik, Universität Bielefeld Postfach 100131, D–33501 Bielefeld, Germany

Abstract

In this paper we develop a numerical method for computing higher order local approximations of invariant manifolds, such as stable, unstable or center manifolds near steady states of a dynamical system. The underlying system is assumed to be large in the sense that a large sparse Jacobian at the equilibrium occurs, for which only a linear (black box) solver and a low dimensional invariant subspace is available, but for which methods like the QR–Algorithm are considered to be too expensive. Our method is based on an analysis of the multilinear Sylvester equations for the higher derivatives which can be solved under certain nonresonance conditions. These conditions are weaker than the standard gap conditions on the spectrum which guarantee the existence of the invariant manifold. The final algorithm requires the solution of several large linear systems with a bordered Jacobian. To these systems we apply a block elimination method recently developed by Govaerts and Pryce (1991, 1993).

Keywords: Invariant manifolds, numerical methods, higher order approximations, large dynamical systems.

Mathematics Subject Classification (1991): 34C35, 58F14, 65F50, 65H17.

1 Introduction

Our primary goal is the local approximation of invariant manifolds near a steady state of a smooth dynamical system

u̇ = F (u), F : ℝ^{m} ⟶ ℝ^{m}.  (1.1)

We assume that the spectrum of the linearization DF (u_{0}) at the equilibrium u_{0} can be decomposed as

σ(DF (u_{0})) = σ_{l} ∪ σ_{h}, σ_{l} ∩ σ_{h} = ∅.  (1.2)

∗supported by DFG Schwerpunktprogramm “Ergodentheorie, Analysis und effiziente Simulation dynamischer Systeme”.

Here σ_{l} contains a few eigenvalues (usually close to or on the imaginary axis) and the corresponding invariant subspace E_{l} ⊂ ℝ^{m} is known, whereas σ_{h} contains the remaining (usually large) part of the spectrum with an invariant subspace E_{h} ⊂ ℝ^{m} which is generally unknown. The indices ‘l’ and ‘h’ will be used to indicate low and high dimensional objects.

Under some well–known gap conditions on the spectrum, such as the following (cf. [38])

Re σ_{h} < α < β < Re σ_{l} and α < r · β  (GC)_{r}

for some r ∈ ℕ, there exists a locally invariant C^{r}–manifold which may be written as a graph over E_{l}

W = {(x, g(x)) : x ∈ V_{l} ⊂ E_{l}}, V_{l} a suitable neighborhood in E_{l},  (1.3)

for some C^{r}–function g : V_{l} → E_{h}. Depending on the position of the spectral set σ_{l}, this manifold is called strongly unstable (0 < Re σ_{l}), center–unstable (0 ≤ Re σ_{l}), center (0 = Re σ_{l}) or pseudo unstable (β < Re σ_{l}, β < 0), and further types occur under time reversal of system (1.1) (cf. [3], [38]).

Assuming that the remaining spectrum is stable, the dynamics of (1.1) close to u_{0} can be reduced to the low dimensional system

ẋ = P_{l}F (x, g(x)),  (1.4)

where P_{l} is the projection onto E_{l} along E_{h}. Reduction processes of this type are very well studied analytically and the applications to ordinary and partial differential equations are abundant (for example [7], [16], [21], [26], [29], [41]).

In this paper we develop a numerical method for computing the higher derivatives of the function g under the sole assumption that σ_{l}, E_{l} are known, that a linear solver for the translates DF (u_{0})−λI is provided and that the same information on the transpose DF (u_{0})^{T} is available. Such a situation typically occurs for spatial discretizations of time–dependent partial differential equations. In some of these applications even the requirement on the transpose might be too restrictive.

In section 2 we summarize some well–known facts from invariant manifold theory and we formally set up the equations for the higher derivatives D^{n}g, n = 1, …, r at the equilibrium. These are linear equations for the multilinear forms and, since for the first derivative one obtains the well–known matrix equation of Sylvester (see [40]), we call them the multilinear Sylvester equations. We show that a suitable reformulation of the multilinear Sylvester equation does not only determine D^{n}g but at the same time the n–th derivative of the reduced vector field in (1.4).

In section 3 we show that the multilinear Sylvester equations can be solved via a generalization of Rosenblum’s formula (cf. [36]) under the nonresonance conditions

{ ∑_{ν=1}^{n} λ_{ν} : λ_{ν} ∈ σ_{l}} ∩ σ_{h} = ∅  (NRC)_{n}

for n = 1, …, r. These conditions are weaker than the gap condition (GC)_{r}. In the linearization theory the nonresonance conditions are imposed on any spectral decomposition σ_{l} ∪ σ_{h} and this leads to the well–known Sternberg conditions (see [32], [37], [17], [23]).

While the complex integrals in the generalized Rosenblum formula have the advantage of being generalizable to an infinite dimensional setting, the formula itself seems to be of little practical use.

As an alternative we set up a solution method which generalizes the well–known Bartels– Stewart algorithm for the linear Sylvester equation, (cf. [11]). However, this algorithm still does not meet the above requirements because it works with a representation of the Jacobian DF (u_{0}) on the large space E_{h}.

Therefore, in section 4 we develop a new version of this algorithm in which all computations are based on solving linear systems with bordered Jacobians of the form

Γ = ( \begin{matrix}DF (u_{0}) − λI & φ \\ ψ^{T} & 0 \end{matrix} ).  (1.5)

Here the number of columns in φ and ψ equals the dimension of the low dimensional space E_{l} and λ is a shift parameter of the form λ = ∑_{ν=1}^{n} λ_{ν}, λ_{ν} ∈ σ_{l}.

For these linear systems we use the so–called mixed block elimination method developed by Govaerts and Pryce ([12], [14]). The remarkable feature of their algorithm is that it works with a black box solver for the large matrix DF (u_{0}) − λI. The only requirements are that the complete matrix Γ in (1.5) is well–conditioned and that the solver for the large matrix DF (u_{0}) − λI is normwise backward stable in the sense of [18]; the matrix DF (u_{0}) − λI itself however may be almost singular.

Some generalizations of our algorithm to discrete dynamical systems are indicated in sections 3 and 4, and an application to the Brusselator system is provided in section 5.

Let us finally remark that the overall algorithm is not a completely new approach but is inherent in many concrete analytical computations (see for example the projector method in [26]). But we think it is useful to show that there is a general theory for the multilinear Sylvester equation and that this leads to a general solution algorithm which works under the linear algebra restrictions described above.

2 Prerequisites from invariant manifold theory

2.1 Locally invariant manifolds in the case of general spectral decomposition

Consider the autonomous first order differential equation

\frac{du}{dt} =^{def} u̇ = F (u),  (2.1)

where F is in C^{r}(Ω, ℝ^{m}), Ω is an open subset of ℝ^{m} and r is a natural number. Throughout, we assume that u_{0} ∈ Ω is a steady state of (2.1) and we use A to denote the Jacobian DF (u_{0}) at u_{0}. Locally invariant C^{r}–manifolds near u_{0} are known to exist under certain gap conditions on the spectrum σ(A). For the following theorem, adopted from Shub (1987) [38] (see also [20]), we assume that real numbers α < β < 0 < γ < δ are given, such that the sets

\begin{aligned}σ_{s} & := & \{λ ∈ σ(A) : Re λ < α\}, \\ σ_{c} & := & \{λ ∈ σ(A) : β < Re λ < γ\}, \\ σ_{u} & := & \{λ ∈ σ(A) : δ < Re λ\} \end{aligned}

Fig. 1 form a disjoint decomposition of the spectrum σ(A) (Figure 1).

Let E_{s}, E_{c} and E_{u} denote the corresponding invariant subspaces of ℝ^{m} with respect to this decomposition, i.e. for ∗ ∈ {s, c, u} the space E_{∗} is the maximal subspace of ℝ^{m} which satisfies

A(E_{∗}) ⊂ E_{∗} and σ(A|_{E_{∗}}) = σ_{∗}.

Theorem 2.1: Under the above assumptions there exist zero neighborhoods V_{s} ⊂ E_{s}, V_{c} ⊂ E_{c}, V_{u} ⊂ E_{u} and Lipschitz functions g_{ss} : V_{s} → V_{c} ⊕ V_{u}, g_{sc} : V_{s} ⊕ V_{c} → V_{u}, g_{c} : V_{c} → V_{s} ⊕ V_{u}, g_{cu} : V_{c} ⊕ V_{u} → V_{s} as well as g_{uu} : V_{u} → V_{s} ⊕ V_{c}, such that each of these five functions maps 0 onto 0 and has a vanishing derivative at the origin. The shifted graphs of these functions

\begin{aligned}W_{ss}^{V}(u_{0}) & := & \{u_{0} + v_{ss} + g_{ss}(v_{ss}) : v_{ss} ∈ V_{s}\} , \\ W_{sc}^{V}(u_{0}) & := & \{u_{0} + v_{sc} + g_{sc}(v_{sc}) : v_{sc} ∈ V_{s} ⊕ V_{c}\} , \\ W_{c}^{V} (u_{0}) & := & \{u_{0} + v_{c} + g_{c}(v_{c}) : v_{c} ∈ V_{c}\} , \\ W_{cu}^{V}(u_{0}) & := & \{u_{0} + v_{cu} + g_{cu}(v_{cu}) : v_{cu} ∈ V_{c} ⊕ V_{u}\} and \\ W_{uu}^{V}(u_{0}) & := & \{u_{0} + v_{uu} + g_{uu}(v_{uu}) : v_{uu} ∈ V_{u}\} \end{aligned}

are locally invariant manifolds for the flow of (2.1).

Moreover, on suitable neighborhoods one can choose g_{ss} and g_{uu} as C^{r}–functions. For the C^{r} smoothness of the three remaining functions the following spectral gap conditions are sufficient: g_{cu} is in C^{r} if

α < r · β  (2.2)

and g_{sc} is in C^{r} under the condition

r · γ < δ.  (2.3)

If both gap conditions (2.2) and (2.3) are satisfied, then also g_{c} may be chosen as a C^{r}– function.

We recall that the local invariance of any of the manifolds W in Theorem 2.1 holds with respect to the neighborhood V := u_{0} + V_{s} ⊕ V_{c} ⊕ V_{u}, which means that the flow ϕ of the system (2.1) satisfies

v ∈ M, ϕ(t, v) ∈ V for t_{0} < t < t_{1}, t_{0} < 0 < t_{1} ⟹ ϕ(t, v) ∈ M for t_{0} < t < t_{1}.

The existence and smoothness of the strong stable and strong unstable manifold have been known for a long time (for example Irwin (1972) [22] and 1980 [23] in the discrete case; see also [6]). For the smoothness of the generalized center stable, center and center unstable manifold, see [20] and [38]. A short overview of various methods of proof can be found in [3], and a reference to Hölder continuity is [19].

2.2 Taylor expansion

In the next sections we assume that the spectrum of the Jacobian A splits into two parts, σ_{l} and σ_{h} = σ(A) \ σ_{l}, with

σ_{l} = \overline{σ_{l}} and σ_{h} = \overline{σ_{h}},  (2.4)

where the bar denotes the conjugate complex. Let ℝ^{m} = E_{l} ⊕ E_{h} be the associated decomposition into invariant subspaces of A and let P_{l}, P_{h} be the corresponding projectors. Moreover, we set

F_{l}(u) := P_{l}F (u) and F_{h}(u) := P_{h}F (u).  (2.5)

We assume the existence of a locally invariant C^{r} manifold

W = {(x, g(x)) : x ∈ V_{l}},  (2.6)

where g ∈ C^{r}(V_{l}, E_{h}) is defined on an open set V_{l} ⊂ E_{l} containing P_{l}u_{0} =: u_{0,l} and satisfies g(u_{0,l}) = P_{h}u_{0} =: u_{0,h}.

It is well–known that W is locally invariant if and only if g(x) satisfies the differential equation

g′(x)F_{l}(x, g(x)) = F_{h}(x, g(x)) = P_{h}F (x, g(x)), x ∈ V_{l}.  (2.7)

This means that the vector field restricted to W is tangent to the manifold. On this manifold the dynamics is then determined by the reduced system

ẋ = F_{l}(x, g(x)) = P_{l}F (x, g(x)).  (2.8)

We are interested in the local approximation of the manifold W and the reduced vector field, hence we expand F and g in their Taylor-series at the steady state. In the following it is convenient to introduce the function

g̃(x) := (x, g(x)), x ∈ V_{l},

and to abbreviate evaluations at u_{0}, u_{0,l} by an upper index ‘^{0}’, for example

D^{ν}F^{0} = D^{ν}F (u_{0}), D^{ν}g^{0} = D^{ν}g(u_{0,l}), D^{ν}g^{0} = D^{ν}g̃(u_{0,l}).

In order to express higher derivatives of equation (2.7) we employ the higher order Leibniz and chain rule (cf. [2], [32], [1], notice the inconsistencies in all of them). To avoid factorials and to reduce the number of indices we write down simplified versions of the formulae. If B is some bilinear form and f_{1}, f_{2} are smooth maps, then the Leibniz rule for the n–th derivative of B(f, g)(x) reads as follows

D^{n}B(f_{1}, f_{2})(x)(h_{1}, …, h_{n}) = ∑_{π∈P_{n}} B(D^{|π_{1}|}f_{1}(x)h_{π_{1}}, D^{|π_{2}|}f_{2}(x)h_{π_{2}}),  (2.9)

where the sum is taken over all partitions π = (π_{1}, π_{2}) ∈ P_{n} of {1, …, n} into two disjoint subsets π_{1} and π_{2} (with π_{1} = ∅ or π_{2} = ∅ included). For any finite subset ω = {ω_{1}, …, ω_{n}} of ℕ, let n = |ω| ∈ ℕ_{0} denote its cardinality, and ordering the elements according to ω_{1} < ⋯ < ω_{n} we use the notation

h_{ω} = (h_{ω_{1}}, …, h_{ω_{n}}), h_{∅} = 1.

The precise definition of P_{n} is

P_{n} = {π = (π_{1}, π_{2}) : π_{1}, π_{2} ⊂ {1, …, n}, π_{1} ∪ π_{2} = {1, …, n}, π_{1} ∩ π_{2} = ∅}.

For the chain rule of order n ≥ 1 it is convenient to consider nontrivial, ordered partitions π = (π_{1}, …, π_{k}) of a given subset ω ⊂ ℕ, n = |ω| as follows

∅ ≠ π_{i} ⊂ ω, ⋃_{i=1}^{k} π_{i} = ω, π_{i} ∩ π_{j} = ∅ and min π_{i} < min π_{j} for all i < j.

Let Q_{k}(ω) denote the set of all such partitions π = (π_{1}, …, π_{k}) and let |π| = k denote their length. Moreover, we set Q(ω) = ⋃_{k=1}^{n} Q_{k}(ω), Q(∅) = {∅} for formal reasons, and we abbreviate Q(n) = Q({1, …, n}). With these settings we may write the chain rule of order n for given maps f_{1}, f_{2} and a finite subset ω ⊂ ℕ, |ω| = n in the following form

D^{n}(f_{1} ∘ f_{2})(x)h_{ω} = ∑_{π∈Q(ω)_{(k:=|π|)}} D^{k}f_{1}(f_{2}(x)) (D^{|π_{1}|}f_{2}(x)h_{π_{1}}, …, D^{|π_{k}|}f_{2}(x)h_{π_{k}}) . (2.10)

The proof is easily accomplished by an induction argument on n = |ω|.

Applying the Leibniz rule (2.9) and the chain rule (2.10) the n–th derivative of equation (2.7) reads

∑_{π∈P_{n}_{(π=(π_{1},π_{2}))}} ∑_{τ∈Q(π_{2})_{(k:=|τ|)}} D^{|π_{1}|+1}g^{0}[h_{π_{1}}, D^{k}F_{l}^{0}(D^{|τ_{1}|}g^{0}h_{τ_{1}}, …, D^{|τ_{k}|}g^{0}h_{τ_{k}})]

= ∑_{π∈P_{n}} D^{|π_{1}|+1}g^{0}[h_{π_{1}}, D^{|π_{2}|}(F_{l} ∘ g̃)^{0}h_{π_{2}}] = D^{n}((g′) · (F_{l} ∘ g̃))^{0}(h_{1}, …, h_{n}) (2.11)

= D^{n}(F_{h} ∘ g̃)^{0}(h_{1}, …, h_{n}) = ∑_{π∈Q(n)_{(k:=|π|)}} D^{k}F_{h}^{0}(D^{|π_{1}|}g^{0}h_{π_{1}}, …, D^{|π_{k}|}g^{0}h_{π_{k}})

for (h_{1}, …, h_{n}) ∈ E_{l}^{n}.

Let L^{n}(E_{1}, …, E_{n}; E_{0}) denote the Banach space of n–multilinear continuous maps from the Banach spaces E_{1}, …, E_{n} into the Banach space E_{0}. If E_{ν} = E_{1}, 2 ≤ ν ≤ n, this space is denoted by L^{n}(E_{1}, E_{0}) and we set L(E_{1}, E_{0}) := L^{1}(E_{1}; E_{0}).

Because of the invariance of the subspaces E_{l} and E_{h} we get

D_{h}F_{l}(u_{0}) = 0 and D_{l}F_{h}(u_{0}) = 0,

where D_{l}, D_{h} denotes differentiation with respect to E_{l} and E_{h}, respectively. By setting

A_{l} := D_{l}F_{l}(u_{0,l}, u_{0,h}) ∈ L(E_{l}, E_{l}),
A_{h} := D_{h}F_{h}(u_{0,l}, u_{0,h}) ∈ L(E_{h}, E_{h}),

we may write

A = A_{l} ⊕ A_{h},

where

σ(A_{l}) = σ_{l} and σ(A_{h}) = σ_{h}.

Equation (2.11) looks almost intractable for general n. But the first cases are easily written down. For n = 1 we have the linear Sylvester equation

Dg^{0} · A_{l} − A_{h} · Dg^{0} = 0.  (2.12)

Since σ(A_{l}) ∩ σ(A_{h}) = σ_{l} ∩ σ_{h} = ∅, this has the unique solution

Dg^{0} = 0  (2.13)

(cf. [40]). Using this the next equations are

D^{2}g^{0}(A_{l}h_{1}, h_{2}) + D^{2}g^{0}(h_{1}, A_{l}h_{2}) − A_{h} · D^{2}g^{0}(h_{1}, h_{2})
= D_{l}^{2}F_{h}^{0}(h_{1}, h_{2})  for all (h_{1}, h_{2}) ∈ E_{l} × E_{l},

D^{3}g^{0}(A_{l}h_{1}, h_{2}, h_{3}) + D^{3}g^{0}(h_{1}, A_{l}h_{2}, h_{3}) + D^{3}g^{0}(h_{1}, h_{2}, A_{l}h_{3}) − A_{h}D^{3}g^{0}(h_{1}, h_{2}, h_{3})

= D_{l}^{3}F_{h}^{0}(h_{1}, h_{2}, h_{3}) + D_{h}D_{l}F_{h}^{0} [h_{1}, D^{2}g^{0}(h_{2}, h_{3})] + D_{h}D_{l}F_{h}^{0} [h_{2}, D^{2}g^{0}(h_{1}, h_{3})]

+ D_{h}D_{l}F_{h}^{0} [h_{3}, D^{2}g^{0}(h_{1}, h_{2})] − D^{2}g^{0} [h_{1}, D_{l}^{2}F_{l}^{0}(h_{2}, h_{3})] − D^{2}g^{0} [h_{2}, D_{l}^{2}F_{l}^{0}(h_{1}, h_{3})]

− D^{2}g^{0} [h_{3}, D_{l}^{2}F_{l}^{0}(h_{1}, h_{2})] for all (h_{1}, h_{2}, h_{3}) ∈ E_{l} × E_{l} × E_{l}.

In the general case by application of (2.13) and F (u_{0}) = 0 we delete in (2.11) the summands corresponding to π_{1} = ∅ and π_{2} = ∅ on the left–hand side. Moreover, we extract the n summands corresponding to |π_{1}| = n − 1, |π_{2}| = 1 from the left–hand side and the summand with |π| = 1, i. e. π = (π_{1}) = ({1, …, n}), from the right–hand side. This yields the following result.

Theorem 2.2: For n = 1, …, r − 1 the n–th derivative X := D^{n}g^{0} ∈ L^{n}(E_{l}, E_{h}) at u_{0,l} satisfies the multilinear Sylvester equation

\begin{aligned}(L(X))(h_{1}, …, h_{n}) & := & ∑_{ν=1}^{n} X(h_{1}, . . , h_{ν−1}, A_{l}h_{ν}, h_{ν+1}, . . , h_{n}) − A_{h}X(h_{1}, …, h_{n}) \\ & = & C_{h,n}(h_{1}, …, h_{n})  for all (h_{1}, …, h_{n}) ∈ E_{l}^{n}, \end{aligned} (2.14)

where C_{h,n} ∈ L^{n}(E_{l}, E_{h}) is determined by the derivatives D^{ν}F_{l}^{0}, ν = 1, …, n−1, and D^{ν}F_{h}^{0}, ν = 1, …, n, as well as by the previous derivatives D^{2}g^{0}, …, D^{n−1}g^{0}.

In the following corollary we reformulate the multilinear Sylvester equation (2.14) in a form useful for later purposes (see section 4.2) and write down an explicit formula for the right– hand side.

Corollary 2.3: (i) Setting Y := Y_{n} := \begin{cases}0 & for n = 1 \\ D^{n}(F_{l} ∘ g̃)^{0} & for n ≥ 2 \end{cases}\} ∈ L^{n}(E_{l}, E_{l}) the multilin-

ear forms X and Y satisfy

L(X) + Y = C_{n},  (2.15)

where C_{n} ∈ L^{n}(E_{l}, ℝ^{m}) is determined by the derivatives D^{ν}F_{l}^{0}, ν = 1, …, n − 1, and D^{ν}F^{0}, ν = 1, …, n, as well as by the previous derivatives D^{2}g^{0}, …, D^{n−1}g^{0}. The explicit form is

\begin{aligned}C_{n}(h_{1}, …, h_{n}) & = & ∑_{π∈Q^{∗}(n)_{(k:=|π|)}} D^{k}F^{0}(D^{|π_{1}|}g^{0}h_{π_{1}}, …, D^{|π_{k}|}g^{0}h_{π_{k}}) \\ & − & ∑_{π∈P_{n}^{∗}} ∑_{τ∈Q(π_{2})_{(k:=|τ|)}}D^{|π_{1}|+1}g^{0}[h_{π_{1}}, D^{k}F_{l}^{0}(D^{|τ_{1}|}g^{0}h_{τ_{1}}, …, D^{|τ_{k}|}g^{0}h_{τ_{k}})], \end{aligned} (2.16)

where we use the abbreviations

P_{n}^{∗} := {π = (π_{1}, π_{2}) ∈ P_{n} : π_{1}, π_{2} ≠ ∅, |π_{2}| ≠ 1} and Q^{∗}(n) := {π ∈ Q(n) : |π| > 1}

The derivatives of g̃ are given by

Dg^{0} = ( \begin{matrix}I \\ Dg^{0} \end{matrix} ) = ( \begin{matrix}I \\ 0 \end{matrix} ) and D^{ν}g^{0} = ( \begin{matrix}0 \\ D^{ν}g^{0} \end{matrix} ) for ν ≥ 2,

and none of the summands in (2.16) contains D^{n}g^{0} or D^{n}g^{0}.

(ii) The right–hand side C_{h,n} in equation (2.14) fulfils

C_{h,n} = P_{h}C_{n}.

C_{h,n} has the same form as C_{n} in (2.16), but with D^{k}F^{0} replaced by D^{k}F_{h}^{0}.

Proof of the corollary: In case n = 1 we have L(X) = 0, see (2.12). For n ≥ 2 using the n-th derivative of equation (2.7) we find

D^{n}((g′) · (F_{l} ∘ g̃))^{0} = D^{n}(F ∘ g̃)^{0} − Y,

and hence the multilinear form C_{n} equals

D^{n}(F ∘ g̃)^{0} − D^{n}((g′) · (F_{l} ∘ g̃))^{0} + L(X).

Now we use F (u_{0}) = 0, (2.13) and the chain rule (2.10) on the first term D^{n}(F ∘ g̃)^{0} and the Leibniz rule (2.9) on the second term (see (2.11)). Extracting the summands as above our conclusion follows. □

If equation (2.14) has a unique solution X ∈ L^{n}(E_{l}, E_{h}) for each C_{h,n} ∈ L^{n}(E_{l}, E_{h}), then obviously equation (2.15) has a unique solution (X, Y) ∈ L^{n}(E_{l}, E_{h}) × L^{n}(E_{l}, E_{l}) for every right–hand side C_{n} ∈ L^{n}(E_{l}, ℝ^{m}) and vice versa. In case of unique solvability, both Theorem 2.2 and Corollary 2.3 (i) define an iterative process for computing the derivatives

D^{1}g^{0}, …, D^{n−1}g^{0}, D^{n}g^{0},

from which the formal Taylor expansion of the function g is obtained. Therefore, in the following sections we analyze the solvability of (2.14) as well as methods for its numerical solution.

There are two ways of computing the derivatives of the reduced vector field F_{l} ∘ g̃ at the steady state u_{0,l}.

First, if instead of (2.14) we solve the enlarged system (2.15) (this will be done in section 4) the derivative Y = D^{n}(F_{l} ∘ g̃)^{0} is obtained as a by–product of the computation of the n–th derivative D^{n}g^{0}.

Second, if the derivative D^{n}g^{0} has been computed we can use it to evaluate the next derivative D^{n+1}(F_{l} ∘ g̃)^{0} as follows. From equation (2.10) we have

D^{n+1}(F_{l} ∘ g̃)(x) = G(x) + D_{h}F_{l}(x, g(x)) · D^{n+1}g(x), x ∈ V_{l}, n ∈ {1, …, r − 1},

with a function G(x) ∈ L^{n+1}(E_{l}, E_{h}) that can be computed from D^{1}g(x), …, D^{n}g(x) and D^{1}F_{l}(x, g(x)), …, D^{n+1}F_{l}(x, g(x)). Since D_{h}F_{l}^{0} = 0, we can evaluate the (n + 1)–th derivative at the steady state through

D^{n+1}(F_{l} ∘ g̃)^{0} = G^{0}.  (2.17)

The first terms in a Taylor expansion of the reduced system (2.8) are given by

\begin{matrix}ξ̇ = A_{l}ξ + \frac{1}{2}D_{l}^{2}F_{l}^{0}ξ^{2} + \frac{1}{6}(D_{l}^{3}F_{l}^{0}ξ^{3} + 3D_{h}D_{l}F_{l}^{0}[ξ, D^{2}g^{0}ξ^{2}]) + \frac{1}{24}(D_{l}^{4}F_{l}^{0}ξ^{4} + \\ 6D_{h}D_{l}^{2}F_{l}^{0}[ξ^{2}, D^{2}g^{0}ξ^{2}] + 3D_{h}^{2}F_{l}^{0}[D^{2}g^{0}ξ^{2}]^{2} + 4D_{h}D_{l}F_{l}^{0}[ξ, D^{3}g^{0}ξ^{3}]) + … \end{matrix} (2.18)

with ξ = x − u_{0,l}.

For a fair comparison of the second method with the first we would have to go up to order n + 1 in the first method. This seems to be more costly than the second approach and, moreover, it requires the conditions of solvability to hold up to order n + 1 rather than n.

In summary, we propose to evaluate the Taylor expansion of order r for the reduced vector field by solving (2.15) successively for n = 2, …, r−1. This automatically yields D^{n}(F_{l} ∘g̃)^{0}, n = 2, …, r − 1 and D^{r}(F_{l} ∘ g̃)^{0} is then computed from (2.17) with n = r − 1.

2.3 The mapping case

We briefly report the analogs of the previous results for a discrete dynamical system

u^{n+1} = F (u^{n}), n ∈ ℤ,  (2.19)

with a regular fixed point u_{0}. Now the spectrum of A := DF (u_{0}) is decomposed into the annuli 0 < |λ| < α and β < |λ| < γ and δ < |λ|, where β < 1 < γ. The analog of Theorem 2.1 then asserts the existence of locally invariant manifolds tangent to the corresponding invariant subspaces, and the gap conditions for C^{r}–smoothness read:

α < β^{r} or γ^{r} < δ,  (2.20)

respectively (see Shub (1987) [38]).

The local invariance of the manifold W = {(x, g(x)) : x ∈ V_{l}} can now be characterized by the functional equation

g(F_{l}(x, g(x))) = F_{h}(x, g(x))  (2.21)

for all x ∈ V_{l}, such that F_{l}(x, g(x)) ∈ V_{l}, and the reduced system is

x^{n+1} = F_{l}(x^{n}, g(x^{n})).

As above, for the expansion of g and the reduced map the higher order derivatives of equation (2.21) are computed and evaluated at the fixed point. One obtains the following analog of Theorem 2.2.

Theorem 2.4: For n ∈ {1, …, r}, the multilinear form X := D^{n}g^{0} ∈ L^{n}(E_{l}, E_{h}) satisfies the equation

\begin{matrix}X(A_{l}h_{1}, …, A_{l}h_{n}) − A_{h} · X(h_{1}, …, h_{n}) \\ = C_{h,n}(h_{1}, …, h_{n})  for all (h_{1}, …, h_{n}) ∈ E_{l}^{n}, \end{matrix}  (2.22)

where C_{h,n} ∈ L^{n}(E_{l}, E_{h}) is determined by derivatives of F_{l} and F_{h} of order less than or equal to n evaluated at u_{0} and by D^{1}g^{0}, …, D^{n−1}g^{0}.

In the case n = 1, one finds C_{h,1} = 0, therefore, Dg^{0} = 0 is the unique solution of the linear Sylvester equation

Dg^{0} · A_{l} − A_{h} · Dg^{0} = 0.

Using this result the next equations for D^{2}g^{0} and D^{3}g^{0} are given by

D^{2}g^{0}(A_{l}h_{1}, A_{l}h_{2}) − A_{h} · D^{2}g^{0}(h_{1}, h_{2}) = D_{l}^{2}F_{h}^{0}(h_{1}, h_{2})

and

D^{3}g^{0}(A_{l}h_{1}, A_{l}h_{2}, A_{l}h_{3}) − A_{h} · D^{3}g^{0}(h_{1}, h_{2}, h_{3})
= D_{l}^{3}F_{h}^{0}(h_{1}, h_{2}, h_{3}) + D_{h}D_{l}F_{h}^{0} [h_{1}, D^{2}g^{0}(h_{2}, h_{3})] + D_{h}D_{l}F_{h}^{0} [h_{2}, D^{2}g^{0}(h_{1}, h_{3})]
+ D_{h}D_{l}F_{h}^{0} [h_{3}, D^{2}g^{0}(h_{1}, h_{2})] − D^{2}g^{0} [A_{l}h_{1}, D_{l}^{2}F_{l}^{0}(h_{2}, h_{3})]
− D^{2}g^{0} [A_{l}h_{2}, D_{l}^{2}F_{l}^{0}(h_{1}, h_{3})] − D^{2}g^{0} [A_{l}h_{3}, D_{l}^{2}F_{l}^{0}(h_{1}, h_{2})].

3 The multilinear Sylvester equation

In this section E_{0} and E_{ν} (ν = 1, …, n) denote subspaces of ℂ^{m} or ℝ^{m}, respectively. For given maps A_{ν} ∈ L(E_{ν}, E_{ν}), B ∈ L(E_{0}, E_{0}) and C ∈ L^{n}(E_{1}, …, E_{n}; E_{0}) we seek a solution X ∈ L^{n}(E_{1}, …, E_{n}; E_{0}) of the multilinear Sylvester equation

\begin{matrix}∑_{ν=1}^{n} X(h_{1}, …, h_{ν−1}, A_{ν}h_{ν}, h_{ν+1}, …, h_{n}) − BX(h_{1}, …, h_{n}) \\ = C(h_{1}, …, h_{n})  for all (h_{1}, …, h_{n}) ∈ E_{1} × ⋯ × E_{n}. \end{matrix}  (3.1)

Notice that this is a linear system for X ∈ L^{n}(E_{1}, …, E_{n}; E_{0}). We start with a theoretical result.

3.1 A generalization of Rosenblum’s formula

In case n = 1 equation (3.1) turns into the linear Sylvester equation

X · A − B · X = C.  (3.2)

Under the spectral condition

σ(A) ∩ σ(B) = ∅

this equation has a unique solution X ∈ L(E_{1}, E_{0}) given by Rosenblum’s formula

X = \frac{1}{2πi} ∫_{G}(ζI − B)^{−1}C(ζI − A)^{−1} dζ,  (3.3)

where G is a simple closed curve enclosing σ(A) and excluding σ(B) (see Rosenblum [36] and [40]). This result can be generalized as follows.

Theorem 3.1: The following statements (i) and (ii) are equivalent:

(i) For every right–hand side C ∈ L^{n}(E_{1}, …, E_{n}; E_{0}), the equation (3.1) has a unique solution X ∈ L^{n}(E_{1}, …, E_{n}; E_{0}).

(ii) The linear maps A_{ν} (ν = 1, …, n) and B fulfil the nonresonance condition of order n:

{ ∑_{ν=1}^{n} λ_{ν} : λ_{ν} ∈ σ(A_{ν})} ∩ σ(B) = ∅.  (3.4)

Under the nonresonance condition (3.4) the unique solution X ∈ L^{n}(E_{1}, …, E_{n}; E_{0}) has the representation

\begin{matrix}X(h_{1}, …, h_{n}) = \frac{1}{(2πi)^{n}} ∫_{G_{n}}⋯ ∫_{G_{1}} ( ∑_{ν=1}^{n} ζ_{ν}I − B)^{−1} \\ · C((ζ_{1}I − A_{1})^{−1}h_{1}, …, (ζ_{n}I − A_{n})^{−1}h_{n}) dζ_{1} ⋯ dζ_{n}, \end{matrix}  (3.5)

where for ν = 1, …, n, G_{ν} is a simple closed curve enclosing σ(A_{ν}) and excluding

σ(B − ∑_{{μ=1 \atop μ≠ν}}^{n} ζ_{μ}I)

for all ζ_{μ} ∈ G_{μ} (μ = 1, …, n, μ ≠ ν).

Remark: The statements of the theorem are true in the real and in the complex case. Moreover, in the special case E_{1} = E_{ν} and A_{1} = A_{ν} for ν = 2, …, n it is easy to see that the symmetry of C implies the symmetry of X under the assumption (3.4). In this case it is also obvious that the gap condition (GC)_{r} implies the nonresonance (3.4) for any n = 1, …, r. Finally, for a different integral representation of the solution X under the assumption of such a spectral gap we refer to Iooss and Adelmeyer [21]; the precise relationship to our approach is worked out in [25].

Proof of the theorem: Let us first assume that (i) holds and that (ii) is false. Then, there exist λ_{ν} ∈ σ(A_{ν}), ν = 1, …, n, such that

β := ∑_{ν=1}^{n} λ_{ν} ∈ σ(B).  (3.6)

Consider first the complex case E_{ν} ⊂ ℂ^{m}. Then we have eigenvectors x_{ν} ∈ E_{ν} and y ∈ E_{0} satisfying x_{ν}^{H}A_{ν} = λ_{ν}x_{ν}^{H} and By = βy. Define X ∈ L^{n}(E_{1}, …, E_{n}; E_{0}) \ {0} by

X(h_{1}, …, h_{n}) := y · ∏_{ν=1}^{n} (x_{ν}^{H} · h_{ν}) for h_{ν} ∈ E_{ν} (ν = 1, …, n)

and find that X is a nontrivial solution of the homogeneous Sylvester equation:

∑_{ν=1}^{n} X(h_{1}, …, h_{ν−1}, A_{ν}h_{ν}, h_{ν+1}, …, h_{n}) = y · ∑_{ν=1}^{n} (λ_{ν}x_{ν}^{H}h_{ν} · ∏_{{μ=1 \atop μ≠ν}}^{n} (x_{μ}^{H} · h_{μ}) )

= βy ∏_{μ=1}^{n} (x_{μ}^{H} · h_{μ}) = BX(h_{1}, …, h_{n}).

In the real case one proceeds in the standard way. First, consider the complexifications E_{ν,ℂ}, A_{ν,ℂ}, B_{ℂ}, C_{ℂ} of E_{ν}, A_{ν}, B, C and construct as above the nontrivial solution X_{ℂ} ∈ L^{n}(E_{1,ℂ}, …, E_{n,ℂ}; E_{0,ℂ}) of the complexified equation. Then, both Re(X_{ℂ}) and Im(X_{ℂ}) restricted to the real spaces E_{ν} are in L^{n}(E_{1}, …, E_{n}; E_{0}), they solve the real equation and they cannot vanish simultaneously.

It remains to show that under condition (ii) formula (3.5) yields a solution for any C so that uniqueness follows from linear algebra. First, notice that by the nonresonance condition we may choose suitable curves G_{ν} ⊂ ℂ which satisfy the requirements of the theorem. For example take small circles enclosing the eigenvalues of A_{ν}. Then, inserting (3.5) into the multilinear Sylvester equation and splitting the terms A_{ν}h_{ν} into (A_{ν} −ζ_{ν}I)h_{ν} +ζ_{ν}h_{ν} implies

∑_{ν=1}^{n} X(h_{1}, …, h_{ν−1}, A_{ν}h_{ν}, h_{ν+1}, …, h_{n}) − BX(h_{1}, …, h_{n}) =

\frac{1}{(2πi)^{n}} ∫_{G_{n}} ⋯∫_{G_{1}} ∑_{ν=1}^{n} (∑_{μ=1}^{n} ζ_{μ}I − B)^{−1}C((ζ_{1}I − A_{1})^{−1}h_{1}, . . , −h_{ν}, . . , (ζ_{n}I − A_{n})^{−1}h_{n})

+ C((ζ_{1}I − A_{1})^{−1}h_{1}, …, (ζ_{n}I − A_{n})^{−1}h_{n}) dζ_{1} ⋯ dζ_{n}.

The first integral terms on the right–hand side vanish as can be seen by evaluating the integral with respect to G_{ν}. The second term yields C(h_{1}, …, h_{n}) due to Cauchy’s integral formula.

In the real case, the spectra are symmetric with respect to the real axis, and, by taking symmetric curves G_{ν}, our formula defines a real multilinear form X. □

We briefly mention the analog of Theorem 3.1 for the mapping case. The multilinear Sylvester equation is now of the form (cf. (2.22))

X(A_{1}h_{1}, …, A_{n}h_{n}) − B · X(h_{1}, …, h_{n}) = C(h_{1}, …, h_{n})  for h_{ν} ∈ E_{ν}  (3.7)

and the nonresonance condition reads

{ ∏_{ν=1}^{n} λ_{ν} : λ_{ν} ∈ σ(A_{ν})} ∩ σ(B) = ∅.  (3.8)

Theorem 3.2: Equation (3.7) has a unique solution X for any right–hand side C if and only if the nonresonance condition (3.8) holds. In that case X has the integral representation

\begin{matrix}X(h_{1}, …, h_{n}) = \frac{1}{(2πi)^{n}} ∫_{G_{n}}⋯ ∫_{G_{1}} ( ∏_{ν=1}^{n} ζ_{ν}I − B)^{−1} \\ · C((ζ_{1}I − A_{1})^{−1}h_{1}, …, (ζ_{n}I − A_{n})^{−1}h_{n}) dζ_{1} ⋯ dζ_{n}. \end{matrix}  (3.9)

Here each curve G_{ν} does not contain 0, has σ(A_{ν}) in its interior and

σ(( ∏_{{μ=1 \atop μ≠ν}}^{n} ζ_{μ})^{−1} · B)

in its exterior for all ζ_{μ} ∈ G_{μ} (μ ≠ ν).

3.2 The tensor approach

In this subsection we reconsider the multilinear Sylvester equation from the point of view of tensor analysis, and we indicate some generalization of our previous result. We use the fact that the pair (L^{n}(E_{1}, …, E_{n}; E_{0}), ⊗) is a tensor product of the vector spaces E_{n}^{∗}, …, E_{1}^{∗}, E_{0} (see e.g. [15], [28]), if the multilinear function ⊗ : E_{n}^{∗} × ⋯ × E_{1}^{∗} × E_{0} → L^{n}(E_{1}, …, E_{n}; E_{0}), (u_{n}, …, u_{0}) ↦ u_{n} ⊗ ⋯ ⊗ u_{0}, is defined by

(u_{n} ⊗ ⋯ ⊗ u_{0})(h_{1}, …, h_{n}) := u_{0} · ∏_{ν=1}^{n} u_{ν}h_{ν} for h_{ν} ∈ E_{ν} (ν = 1, …, n).

Here E_{ν}^{∗} denotes the dual space of the finite dimensional 𝕂–vector space E_{ν} and L_{ν}^{∗} ∈ L(E_{ν}^{∗}, E_{ν}^{∗}) denotes the dual of a linear mapping L_{ν} ∈ L(E_{ν}, E_{ν}) defined by (L_{ν}^{∗}u)(h) = u(L_{ν}h), h ∈ E_{ν}, u ∈ E_{ν}^{∗}. The tensor product L_{n}^{∗} ⊗ ⋯ L_{1}^{∗} ⊗ L_{0} of the linear transformations L_{ν}^{∗} ∈ L(E_{ν}^{∗}, E_{ν}^{∗}), ν = 1, …, n, and L_{0} ∈ (E_{0}, E_{0}) is given by

((L_{n}^{∗} ⊗ ⋯ L_{1}^{∗} ⊗ L_{0})(u_{n} ⊗ ⋯ ⊗ u_{0}))(h_{1}, …, h_{n})
= (L_{n}^{∗}u_{n} ⊗ ⋯ ⊗ L_{1}^{∗}u_{1} ⊗ L_{0}u_{0})(h_{1}, …, h_{n}) = L_{0}u_{0} · ∏_{ν=1}^{n}(L_{ν}^{∗}u_{ν})h_{ν} (3.10)

= L_{0}u_{0} · ∏_{ν=1}^{n} u_{ν}(L_{ν}h_{ν}) = L_{0}(u_{n} ⊗ ⋯ ⊗ u_{0})(L_{1}h_{1}, …, L_{n}h_{n})

for any (u_{n} ⊗ ⋯ ⊗ u_{0}) ∈ L^{n}(E_{1}, …, E_{n}; E_{0}) and (h_{1}, …, h_{n}) ∈ E_{1} × ⋯ × E_{n}.

Let L_{0}, …, L_{n} be maps as above and let p be a polynomial in the variables x_{0}, …, x_{n} with coefficients c_{k_{0},...,k_{n}} ∈ 𝕂

p(x_{0}, …, x_{n}) = ∑_{k_{0},...,k_{n}=0}^{l} c_{k_{0},...,k_{n}} · ∏_{ν=0}^{n} x_{ν}^{k_{ν}}.  (3.11)

With these we associate linear operators in L(L^{n}(E_{1}, . . , E_{n}; E_{0}), L^{n}(E_{1}, . . , E_{n}; E_{0})) through

p(L_{0}, L_{1}^{∗} …, L_{n}^{∗}) = ∑_{k_{0},...,k_{n}=0}^{l} c_{k_{0},...,k_{n}} · L_{n}^{∗ k_{n}} ⊗ ⋯ ⊗ L_{1}^{∗ k_{1}} ⊗ L_{0}^{k_{0}}

and

(T_{L_{0},...,L_{n}}(X))(h_{1}, …, h_{n}) := ∑_{k_{0},...,k_{n}=0}^{l} c_{k_{0},...,k_{n}}L_{0}^{k_{0}}X(L_{1}^{k_{1}}h_{1}, …, L_{n}^{k_{n}}h_{n}), h_{ν} ∈ E_{ν} .

Using equation (3.10) and (L_{ν}^{∗})^{k_{ν}} = (L_{ν}^{k_{ν}})^{∗} implies

p(L_{0}, L_{1}^{∗}, …, L_{n}^{∗}) = T_{L_{0},...,L_{n}} .

With this approach we can express the equations (3.1), (3.7) with the aid of tensor products as

p(B, A_{1}^{∗}, …, A_{n}^{∗}) · X = C,  (3.12)

where in case of (3.1) we have

p(x_{0}, …, x_{n}) = −x_{0} + ∑_{ν=1}^{n} x_{ν} ,  (3.13)

and in case (3.7)

p(x_{0}, …, x_{n}) = −x_{0} + ∏_{ν=1}^{n} x_{ν} .  (3.14)

This reformulation of the problem is helpful in two different respects.

First, we may generalize the Theorem of Stéphanos ([27], [39]) as follows (the proof is by extending the proof of Theorem 2.1 in [28, chapter 2], cf. [25] for details).

Theorem 3.3: Let U_{0}, …, U_{n} be finite dimensional 𝕂–vector spaces where 𝕂 ∈ {ℝ, ℂ}, m_{ν} := dim U_{ν}. Let p be the polynomial from equation (3.11). If the maps L_{ν} ∈ L(U_{ν}, U_{ν}), ν = 0, …, n, have eigenvalues μ_{ν,1}, …, μ_{ν,m_{ν}} (repeated according to their algebraic multiplicity) then the m_{0} ⋯ m_{n} numbers

p(μ_{0,r_{0}}, …, μ_{n,r_{n}}), r_{ν} = 1, …, m_{ν} (ν = 0, …, n),

are the eigenvalues of the linear map p(L_{0}, …, L_{n}) repeated according to their algebraic multiplicity.

Theorems 3.1 and 3.2 are then simple consequences of Theorem 3.3 when applied to (3.12), (3.13), (3.14).

The second aspect is the relationship between tensor products of linear maps and Kronecker products of matrices (see [28]):

For ν = 1, …, n let {e_{ν,1}, …, e_{ν,m_{ν}}} be an ordered basis of E_{ν}, and let {e_{ν,1}^{∗}, …, e_{ν,m_{ν}}^{∗}} be the corresponding dual basis. If M_{ν} ∈ 𝕂^{m_{ν},m_{ν}} denotes the matrix representation of L_{ν} ∈ L(E_{ν}, E_{ν}) then M_{ν}^{T} is the matrix representation of L_{ν}^{∗}. Further, if one chooses the following basis of the space L^{n}(E_{1}, …, E_{n}; E_{0}) = E_{n}^{∗} ⊗ ⋯ ⊗ E_{1}^{∗} ⊗ E_{0}

{e_{n,r_{n}}^{∗} ⊗ ⋯ ⊗ e_{1,r_{1}}^{∗} ⊗ e_{0,r_{0}} : r_{ν} ∈ {1, …, m_{ν}}, ν = 0, …, n}  (3.15)

ordered lexicographically with priority to the first components of (r_{n}, …, r_{0}), then the matrix representation of the tensor product L_{n}^{∗} ⊗ ⋯ L_{1}^{∗} ⊗ L_{0} is given by the Kronecker product

M_{n}^{T} ⊗ ⋯ M_{1}^{T} ⊗ M_{0} ∈ 𝕂^{m̃,m̃},

where m̃ := ∏_{ν=0}^{n} m_{ν}.

Similarly, for the polynomial p defined in (3.11) we use the Kronecker product and define

p(M_{0}, …, M_{n}) = ∑_{k_{0},...,k_{n}=0}^{l} c_{k_{0},...,k_{n}} · M_{n}^{k_{n}} ⊗ ⋯ ⊗ M_{0}^{k_{0}}.

Let vec : L^{n}(E_{1}, …, E_{n}; E_{0}) → 𝕂m̃ be the isomorphism which assigns to each multilinear form X its coordinate vector with respect to the basis (3.15). A recursive definition of this map is given in [25]. It generalizes the well–known vec–function of a matrix as used for example in [33], [27]. With these assumptions and settings, and by identifying the linear maps A_{ν} and B with their matrix representations, equation (3.12) may be written as

p(B, A_{1}^{T} , …, A_{n}^{T}) · vec(X) = vec(C).  (3.16)

For equation (3.13) we use

p(B, A_{1}^{T} , …, A_{n}^{T}) = ∑_{ν=1}^{n} (I_{m_{n}···m_{ν+1}} ⊗ A_{ν}^{T} ⊗ I_{m_{ν−1}···m_{0}}) − I_{m_{n}···m_{1}} ⊗ B,

and for (3.14)

p(B, A_{1}^{T} , …, A_{n}^{T}) = A_{n}^{T} ⊗ ⋯ ⊗ A_{1}^{T} ⊗ I_{m_{0}} − I_{m_{n}···m_{1}} ⊗ B.

Since the matrix p(B, A_{1}^{T} , …, A_{n}^{T}) is of dimension m_{0} ⋯ m_{n}, this approach seems feasible for numerical computations only in case the number n and the dimensions m_{ν} are very small.

3.3 A multilinear Bartels–Stewart algorithm

The integral representation of the solution X in Theorem 3.1 is certainly of limited practical value. Its evaluation by a quadrature rule would require solving a multitude of linear systems. Furthermore, the tensor approach seems to be practicable for very small dimensions only. We therefore give an alternative proof of existence and uniqueness under the nonresonance condition (3.4). Our method of proof can be formulated as an algorithm. It is in fact a generalization of the Bartels–Stewart algorithm for solving the linear Sylvester equation (see [4], [11]).

From now on we work with bases in finite dimensional spaces, hence, we identify E_{ν} with 𝕂^{m_{ν}}, 𝕂 ∈ {ℝ, ℂ}, for some m_{ν} and regard A_{ν}, B as matrices in 𝕂^{m_{ν},m_{ν}} and 𝕂^{m_{0},m_{0}}, respectively.

In the first step we construct unitary matrices U_{ν} to form the complex Schur decompositions

U_{ν}^{H}A_{ν}U_{ν} = A′_{ν} = ( \begin{matrix}A′_{ν,11} & ⋯ & A′_{ν,1m_{ν}} \\ & ... & ⋮ \\ ◯ & & A′_{ν,m_{ν}m_{ν}} \end{matrix} ) .  (3.17)

Notice that in theory this form always exists and in practice is numerically stable ([11]). Furthermore, note that in the later application it requires the solution of a low dimensional eigenvalue problem (m_{ν} ≪ m_{0}).

Then, we perform a linear change of bases in the domain of X and C. Thus equation (3.1) is equivalent to

∑_{ν=1}^{n} X(U_{1}h_{1}, …, U_{ν−1}h_{ν−1}, U_{ν}A′_{ν}h_{ν}, U_{ν+1}h_{ν+1}, …, U_{n}h_{n}) − BX(U_{1}h_{1}, …, U_{n}h_{n})

= C(U_{1}h_{1}, …, U_{n}h_{n}) for all h_{ν} ∈ ℂ^{m_{ν}} (ν = 1, …, n),

or to

∑_{ν=1}^{n} X′(h_{1}, …, h_{ν−1}, A′_{ν}h_{ν}, h_{ν+1}, …, h_{n}) − BX′(h_{1}, …, h_{n}) = C′(h_{1}, …, h_{n}), (3.18)

where we define

X′(h_{1}, …, h_{n}) := X(U_{1}h_{1}, …, U_{n}h_{n}) and C′(h_{1}, …, h_{n}) := C(U_{1}h_{1}, …, U_{n}h_{n}). (3.19)

We consider multi–indices (α_{1}, …, α_{n}) = α ∈ ℕ^{n}, α_{ν} ≤ m_{ν}, in a given strict order α^{1} ≺ … ≺ α^{M}, M := m_{1} ⋯ m_{n}, such that the following holds for every index α and every ν ∈ {1, …, n}:

(α_{1}, …, α_{ν−1}, k, α_{ν+1}, …, α_{n}) ≺ α  for all k ∈ {1, …, α_{ν} − 1}.  (3.20)

For example, take a lexicographic order with arbitrarily chosen priority of the components. Let e_{i} denote the i–th coordinate vector. From equation (3.18) we successively compute

X′_{α^{j}} := X′(e_{α_{1}^{j}}, …, e_{α_{n}^{j}}), j = 1, …, M,

as follows:

\begin{aligned}C′_{α^{j}} & = & ∑_{ν=1}^{n} X′(e_{α_{1}^{j}}, …, A′_{ν}e_{α_{ν}^{j}}, …, e_{α_{n}^{j}}) − BX′_{α^{j}} \\ & = & [( ∑_{ν=1}^{n} A′_{ν,α_{ν}^{j}α_{ν}^{j}})I − B]X′_{α^{j}} + ∑_{ν=1}^{n} ∑_{k=1}^{α_{ν}^{j}−1} A′_{ν,kα_{ν}^{j}}X′_{β,β=(α_{1}^{j},...,α_{ν−1}^{j},k,α_{ν+1}^{j},...,α_{n}^{j}).} \end{aligned} (3.21)

The matrix [(∑_{ν=1}^{n} A′_{ν,α_{ν}^{j}α_{ν}^{j}}) I − B] is regular due to the nonresonance condition (3.4) and because of assumption (3.20) the terms X′_{β} in the second sum on the right–hand side satisfy β ≺ α^{j} and hence are known at stage j. Therefore, X′_{α^{j}} can be computed and is uniquely determined by the above equation. In the complex case we set

X(h_{1}, …, h_{n}) := X′(U_{1}^{H}h_{1}, …, U_{n}^{H}h_{n}) for all h_{ν} ∈ ℂ^{m_{ν}} (ν = 1, …, n)

and find the solution of the original multilinear Sylvester equation (3.1). In the real case we work with the complexifications of A_{ν}. Since (3.1) has a unique solution, the above X provided by the algorithm is then a real multilinear form. As an alternative we can set up an algorithm based on the real Schur decomposition, see section 4.4. Equation (3.21) then becomes a low dimensional Sylvester equation for which the tensor approach is appropriate.

Further remarks in [4], [10] and [11] about the overall algorithm for the linear Sylvester equation apply to the multilinear case in an analogous fashion. For example, it may be appropriate to pretransform the matrix B before the linear systems (3.21) are solved.

As in section 3.1 we sketch here the analog of the above algorithm for the mapping case. Starting with (3.7) we form the Schur decompositions (3.17), define X′ and C′ as in (3.19) and determine the coefficients X′_{α^{j}} from the equations

\begin{aligned}C′_{α^{j}} & = & X′(A′_{1}e_{α_{1}^{j}}, …, A′_{n}e_{α_{n}^{j}}) − BX′_{α^{j}} \\ & = & [( ∏_{ν=1}^{n} A′_{ν,α_{ν}^{j}α_{ν}^{j}})I − B]X′_{α^{j}} + ∑_{β≺α^{j} ; β_{ν} ≤α_{ν}^{j}_{(ν=1,...,n)}} ( ∏_{ν=1}^{n} A′_{ν,β_{ν}α_{ν}^{j}})X′_{β.} \end{aligned}  (3.22)

Each of these equations has a unique solution and the back transformation works as above.

4 Numerical Approximation

4.1 General approach

According to Theorem 3.1 we can determine the derivatives D^{ν}g^{0} = D^{ν}g(u_{0,l}) ∈ L^{ν}(E_{l}, E_{h}), ν = 1, …, n, recursively from (2.14) if the nonresonance conditions (NRC)_{ν}, ν = 1, …, n, are satisfied.

For the practical implementation of the multilinear Bartels–Stewart algorithm from 3.3 we must choose bases φ ∈ ℝ^{m,m_{l}} of E_{l} and Φ ∈ ℝ^{m,m_{h}} of E_{h}. We define

S := (φ, Φ), S^{−1} =: ( \begin{matrix}ψ^{T} \\ Ψ^{T} \end{matrix} ) , ψ ∈ ℝ^{m,m_{l}}, Ψ ∈ ℝ^{m,m_{h}}

and obtain the block diagonal form

S^{−1}AS = ( \begin{matrix}ψ^{T} Aφ & ψ^{T} AΦ \\ Ψ^{T} Aφ & Ψ^{T} AΦ \end{matrix} ) =: ( \begin{matrix}Â_{l} & 0 \\ 0 & Â_{h} \end{matrix} ) ,  (4.1)

where Â_{l}, Â_{h} are the basis representations of A_{l}|_{E_{l}}, A_{h}|_{E_{h}}. The projectors P_{l}, P_{h} are given by

P_{l} = φψ^{T} and P_{h} = ΦΨ^{T} = I − P_{l},

and the function ĝ representing the invariant manifold in the new coordinate system assumes the form

ĝ(x) = Ψ^{T} g(φx), x ∈ U_{0} ⊂ ℝ^{m_{l}},

where U_{0} is a suitable neighborhood of x_{0} := ψ^{T} u_{0}.

With these transformations and the setting D^{n}g^{0} = D^{n}ĝ(x_{0}) the multilinear Sylvester equation (2.14) may now be written as

\begin{matrix}∑_{ν=1}^{n} D^{n}g^{0}(h_{1}, …, h_{ν−1}, Â_{l}h_{ν}, h_{ν+1}, …, h_{n}) − Â_{h}D^{n}g^{0}(h_{1}, …, h_{n}) \\ = Ĉ_{h,n}(h_{1}, …, h_{n})  for all h_{1}, …, h_{n} ∈ ℝ^{m_{l}}, \end{matrix}  (4.2)

with Ĉ_{h,n} given by

Ĉ_{h,n}(h_{1}, …, h_{n}) = ∑_{π∈Q^{∗}(n)_{(k:=|π|)}} Ψ^{T} D^{k}F^{0}(G_{|π_{1}|}h_{π_{1}}, …, G_{|π_{k}|}h_{π_{k}})

− ∑_{π∈P_{n}^{∗}} ∑_{τ∈Q(π_{2})_{(k:=|τ|)}}D^{|π_{1}|+1}g^{0}[h_{π_{1}}, ψ^{T} D^{k}F^{0}(G_{|τ_{1}|}h_{τ_{1}}, …, G_{|τ_{k}|}h_{τ_{k}})].

Here we have used the functions G_{ν} ∈ L^{ν}(ℝ^{m_{l}}, ℝ^{m}) defined by

G_{1} = φ, G_{ν} = ΦD^{ν}g^{0} for ν ≥ 2.  (4.3)

We can now apply the multilinear Bartels–Stewart algorithm to (4.2), but notice that this still involves the matrices Φ, Ψ, Â_{h} and hence is only suitable for small systems.

4.2 The modified algorithm for large systems

In large systems it is reasonable to assume that the matrices Â_{l}, φ, ψ satisfying

Aφ = φÂ_{l}, ψ^{T} A = Â_{l}ψ^{T} , ψ^{T} φ = I

are available, but that any explicit use of Â_{h}, Φ, Ψ should be avoided. It is then useful to introduce the partially transformed mapping

X := D^{n}(g ∘ φ)(x_{0}) = ΦD^{n}ĝ(x_{0}) ∈ L^{n}(ℝ^{m_{l}}, ℝ^{m})  (4.4)

which equals G_{n} from (4.3) in case n ≥ 2 and is the unique solution of the blown up system

\begin{aligned}∑_{ν=1}^{n} X(h_{1}, …, h_{ν−1}, Â_{l}h_{ν}, h_{ν+1}, …, h_{n}) − AX(h_{1}, …, h_{n}) & = & ΦĈ_{h,n}(h_{1}, …, h_{n}) \\ ψ^{T} X(h_{1}, …, h_{n}) & = & 0. \end{aligned} (4.5)

Notice that the second equation enforces the range of X to be in E_{h}. Moreover, in the expression ΦĈ_{h,n} the terms involving ΦΨ^{T} can be replaced by I − φψ^{T} and the quantities ΦD^{ν}g^{0} = G_{ν}, 2 ≤ ν ≤ n − 1, are known at this stage.

Equation (4.5) is an overdetermined system which by construction has a unique solution. In practice, however, rounding errors or inaccuracies in the matrices A′_{l}, φ, ψ may destroy this property. Therefore, we propose to increase not only the number of equations but also the number unknowns and to solve an equivalent square system

\begin{aligned}∑_{ν=1}^{n} X(h_{1}, . . , h_{ν−1}, Â_{l}h_{ν}, h_{ν+1}, . . , h_{n}) + (φY −AX)(h_{1}, …, h_{n}) & = & ΦĈ_{h,n}(h_{1}, …, h_{n}) \\ ψ^{T} X(h_{1}, …, h_{n}) & = & 0. \end{aligned} (4.6)

It is easy to see that this system is regular and that the unique solution is given by the above X and Y = 0. Therefore it is possible to use ∥Y ∥ for an a–posteriori test of the accuracy.

As an alternative we can simply replace the leading terms (I − φψ^{T})D^{k}F^{0} in ΦĈ_{h,n} on the right–hand side of (4.6) by D^{k}F^{0}, which corresponds to the change from C_{h,n} to C_{n} in Corollary 2.3. This does not affect X but changes the multilinear form Y to

Y (h_{1}, …, h_{n}) := ∑_{π∈Q^{∗}(n)_{(k:=|π|)}} ψ^{T} D^{k}F^{0}(G_{|π_{1}|}h_{π_{1}}, …, G_{|π_{k}|}h_{π_{k}}).  (4.7)

In case n = 1 we have Y = 0. For n ≥ 2 using the chain rule (2.10) we find

Y = D^{n}(ψ^{T} F ∘ (φ + g ∘ φ))(x_{0}).  (4.8)

Hence the multilinear form Y equals the term Y_{n} of Corollary 2.3 in our current coordinate system and the pair (X, Y) ∈ L^{n}(ℝ^{m_{l}}, ℝ^{m}) × L^{n}(ℝ^{m_{l}}, ℝ^{m_{l}}) as in (4.4) and (4.7) is the unique solution of the system

\begin{aligned}∑_{ν=1}^{n} X(h_{1}, . . , h_{ν−1}, Â_{l}h_{ν}, h_{ν+1}, . . , h_{n}) + (φY −AX)(h_{1}, …, h_{n}) & = & Ĉ_{n}(h_{1}, …, h_{n}) \\ ψ^{T} X(h_{1}, …, h_{n}) & = & 0, \end{aligned} (4.9)

with Ĉ_{n} given by

Ĉ_{n}(h_{1}, …, h_{n}) = ∑_{π∈Q^{∗}(n)_{(k:=|π|)}} D^{k}F^{0}(G_{|π_{1}|}h_{π_{1}}, …, G_{|π_{k}|}h_{π_{k}})

− ∑_{π∈P_{n}^{∗}} ∑_{τ∈Q(π_{2})_{(k:=|τ|)}}G_{|π_{1}|+1}[h_{π_{1}}, ψ^{T} D^{k}F^{0}(G_{|τ_{1}|}h_{τ_{1}}, …, G_{|τ_{k}|}h_{τ_{k}})].

As in section 3.3 we proceed with the Schur decomposition of the small matrix Â_{l}:

U^{H} Â_{l}U = A′_{l}, A′_{l} upper triangular.  (4.10)

Then, introduce the transformed mappings

X′(h_{1}, …, h_{n}) := X(Uh_{1}, …, Uh_{n}),
Y′(h_{1}, …, h_{n}) := Y (Uh_{1}, …, Uh_{n}),  (4.11)
C′_{n}(h_{1}, …, h_{n}) := Ĉ_{n}(Uh_{1}, …, Uh_{n}).

With these settings (4.9) is equivalent to

\begin{aligned}∑_{ν=1}^{n} X′(h_{1}, . . , h_{ν−1}, A′_{l}h_{ν}, h_{ν+1}, . . , h_{n}) + (φY′−AX′)(h_{1}, …, h_{n}) & = & C′_{n}(h_{1}, …, h_{n}) \\ ψ^{T} X′(h_{1}, …, h_{n}) & = & 0. \end{aligned} (4.12)

Ordering the multi–indices (α_{1}^{j}, …, α_{n}^{j}) = α^{j} ∈ {1, …, m_{l}}^{n}, j = 1, …, M := (m_{l})^{n}, according to (3.20) we can successively compute

X′_{α^{j}} := X′(e_{α_{1}^{j}}, …, e_{α_{n}^{j}}), j = 1, …, M, e_{i} = i–th unit vector in ℂ^{m_{l}},

from the linear systems

( \begin{matrix}[( ∑_{ν=1}^{n}(A′_{l})_{α_{ν}^{j}}_{α_{ν}^{j}})I − A] & φ \\ ψ^{T} & 0 \end{matrix} ) · ( \begin{matrix}X′_{α^{j}} \\ Y′_{α^{j}} \end{matrix} ) = ( \begin{matrix}R_{α^{j}} \\ 0 \end{matrix} ) .  (4.13)

with

R_{α^{j}} = C′_{n,α^{j}} − ∑_{ν=1}^{n} ∑_{k=1}^{α_{ν}^{j}−1}(A′_{l})_{kα_{ν}^{j}}X′_{β,β=(α_{1}^{j},...,α_{ν−1}^{j},k,α_{ν+1}^{j},...,α_{n}^{j}).}  (4.14)

The C′_{n}–terms in (4.14) can be computed with the help of G_{1} = φ and the previous derivatives G_{ν} = ΦD^{ν}g^{0}, ν = 2, …, n − 1, via

\begin{aligned}C′_{n}(h_{1}, . . , h_{n}) & = & ∑_{π∈Q^{∗}(n)_{(k:=|π|)}} D^{k}F^{0}(G_{|π_{1}|}Uh_{π_{1}}, …, G_{|π_{k}|}Uh_{π_{k}}) \\ & − & ∑_{π∈P_{n}^{∗}} ∑_{τ∈Q(π_{2})_{(k:=|τ|)}}G_{|π_{1}|+1}[Uh_{π_{1}}, ψ^{T} D^{k}F^{0}(G_{|τ_{1}|}Uh_{τ_{1}}, . . , G_{|τ_{k}|}Uh_{τ_{k}})], \end{aligned} (4.15)

where we use the notation Uh_{ω} = (Uh_{ω_{1}}, …, Uh_{ω_{l}}), ω = (ω_{1}, …, ω_{l}).

Finally, we transform back to obtain the n–th derivative in the original coordinates

X(h_{1}, …, h_{n}) := X′(U^{H}h_{1}, …, U^{H}h_{n})  for h_{1}, …, h_{n} ∈ ℝ^{m_{l}}.  (4.16)

From the unique solvability of equation (4.5) it follows that this X is indeed a real multilinear form. Our complete algorithm is now defined by the formulae (4.10) and (4.13) together with (4.14), (4.15) and (4.16).

Several remarks concerning the practical implementation follow:

First, as in (4.16) we obtain by

Y (h_{1}, …, h_{n}) := Y′(U^{H}h_{1}, …, U^{H}h_{n}) for h_{1}, …, h_{n} ∈ ℝ^{m_{l}}

a real multilinear form. According to equation (4.8) for n ≥ 2 (the only case of interest) this multilinear form is the n-th derivative of the right–hand side of the reduced system in our current coordinate system. Thus our algorithm computes the unique solution of equation (2.15) and, as remarked at the end of section 2.2, we obtain as a by–product of the computation of D^{n}(g ∘ φ) the n-th derivative of the reduced system.

Second, due to the symmetry of the multilinear forms C′_{h,n} and X′ the number of coefficients X′_{α^{j}} to be computed reduces considerably: If the algorithm has reached some stage j ∈ {1, …, M} then the computation of X′_{α^{j}} may be skipped if α^{j} is a permutation of some β ≺ α^{j}. The number of systems (4.13) to be solved then reduces from

(m_{l})^{n} to ( \begin{matrix}n + m_{l} − 1 \\ m_{l} − 1 \end{matrix} ) .

We do not organize the index arrangement here but we show the savings for the bilinear case in section 4.4.

Finally, through the triangular Schur decomposition we have switched to a complex calculation even in case of a real multilinear Sylvester equation. One possibility which avoids complex numbers will be indicated for the bilinear case in section 4.4.

4.3 Solution of bordered systems

Linear systems with bordered matrices and large but sparse upper left block are quite typical in numerical continuation and bifurcation calculations (see [5], [8], [9], [13], [24], [30], [35]).

The upper left block is of the form λI − A, λ ∈ ℝ or λ ∈ ℂ, and it may be singular or almost singular for example if A has a zero eigenvalue in σ_{l}.

Various approaches for solving such systems have been developed and we propose to use the mixed block elimination method of Govaerts and Pryce (see [12], [14]).

In [12] linear systems of the type

m m_{l}  1  1
\begin{matrix}m \\ m & _{l} \end{matrix} ( \begin{matrix}A \\ C \end{matrix} \begin{matrix}B \\ D \end{matrix} ) · ( \begin{matrix}x \\ y \end{matrix} ) = ( \begin{matrix}f \\ g \end{matrix} ) \begin{matrix}m \\ m & _{l} \end{matrix} or Γz = h for short

(4.17) are studied for the case m_{l} = 1. The assumption is that Γ is well–conditioned and that black box solvers S for A and S^{∗} for A^{T} are given which are normwise backward stable (see [18]). This means there is a modest constant C_{S}, such that, for any RHS r ∈ ℝ^{m}, S(r) is an approximation of the solution of As = r satisfying

(A + E)S(r) = r + e

for some perturbation matrix E and vector e such that

∥E∥ ≤ eps C_{S}∥A∥ and ∥e∥ ≤ eps C_{S}∥r∥,

eps being the floating point roundoff unit (∼ 2.2 ∗ 10^{−16} for our example in section 5).

The method mixes the Block Crout and Doolittle factorization with one iterative refinement and it needs one call of S^{∗} and two calls of S. It is shown in [12] that this yields a stable procedure BEM for solving (4.17) provided the condition number of A is at most of the order O(\frac{1}{eps}). This algorithm BEM is generalized in a recursive way to systems with wider borders and then is called BEMW in [14]. Typically, 1 ≤ m_{l} ≪ m holds and the method is stable (with a constant growing only linearly in m_{l}) if the leading principal sub–matrices of Γ of orders m + 1, …, m + m_{l} are not too ill–conditioned. Altogether, BEMW requires m_{l} solves with A^{T} and m_{l} + 1 solves with A and an additional number of flops of the order

4mm_{l}(2m_{l} + 1) + \frac{8}{3}m_{l}^{3}

if real matrices are considered.

4.4 The bilinear case

We now exploit the details of the case n = 2 and show how the real Schur decomposition can be used to compute the bilinear form X = D^{2}(g ∘ φ)^{0} = ΦD^{2}g^{0} ∈ L^{2}(ℝ^{m_{l}}, ℝ^{m}). Of course, this can be generalized to arbitrary n.

For n = 2 equation (4.9) reads

\begin{aligned}X(Â_{l}h_{1}, h_{2}) + X(h_{1}, Â_{l}h_{2}) − AX(h_{1}, h_{2}) + φY (h_{1}, h_{2}) & = & Ĉ_{2}(h_{1}, h_{2}) \\ & = & D^{2}F^{0}(φh_{1}, φh_{2}) \\ ψ^{T} X(h_{1}, h_{2}) & = & 0. \end{aligned}  (4.18)

According to (4.7), (4.8) the solution fulfils Y (h_{1}, h_{2}) = ψ^{T} D^{2}F^{0}(φh_{1}, φh_{2}).

For the small matrix Â_{l} we compute the real Schur decomposition

U^{T} Â_{l}U = A′_{l} = ( \begin{matrix}a_{11} & ⋯ & a_{1p} \\ & ... & ⋮ \\ ◯ & & a_{pp} \end{matrix} ) .  (4.19)

Here U ∈ ℝ^{m_{l},m_{l}} is orthogonal and A′_{l} is upper quasi–triangular, where the blocks a_{ii} either denote a 1–by–1 matrix or a 2–by–2 matrix corresponding to a pair of complex conjugate eigenvalues. Defining X′, Y′ and C′_{2} as in (4.11) and (4.15) we may rewrite (4.18) as

\begin{aligned}X′(A′_{l}h_{1}, h_{2}) + X′(h_{1}, A′_{l}h_{2}) − AX′(h_{1}, h_{2}) + φY′(h_{1}, h_{2}) & = & C′_{2}(h_{1}, h_{2}) \\ & = & D^{2}F^{0}(φ Uh_{1}, φ Uh_{2}) \\ ψ^{T} X′(h_{1}, h_{2}) & = & 0. \end{aligned} (4.20)

Let the mapping τ : {1, …, p} → {1, 2} define the sizes of the blocks a_{ik} via a_{ik} ∈ ℝ^{τ(i),τ(k)}. Then we partition the bilinear forms X′ and C′_{2} into small bilinear forms by using the same pattern as in (4.19).

For example, in the case τ(1) = 2, τ(2) = 1 and τ(3) = 2 the bilinear form X′ may be pictured as follows

X′ = (X′_{ik}) =

[ \begin{matrix}X′_{11} & X′_{12} & X′_{13} & X′_{14} & X′_{15} & ⋯ & X′_{1m_{l}} \\ X′_{21} & X′_{22} & X′_{23} & X′_{24} & X′_{25} & ⋯ & X′_{2m_{l}} \\ ⋮ & ⋮ & X′_{33} & X′_{34} & X′_{35} & ⋯ & X′_{3m_{l}} \\ ⋮ & ⋮ & ⋮ & X′_{44} & X′_{45} & ⋯ & X′_{4m_{l}} \\ ⋮ & ⋮ & ⋮ & X′_{54} & X′_{55} & ⋯ & X′_{5m_{l}} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ · & · & · & · & · & ⋯ & X′_{m_{l}m_{l}} \end{matrix}] = [ \begin{matrix}x_{11} & x_{12} & x_{13} & ⋯ & x_{1p} \\ ⋮ & x_{22} & x_{23} & ⋯ & x_{2p} \\ ⋮ & ⋮ & x_{33} & ⋯ & x_{3p} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ · & · & · & ⋯ & x_{pp} \end{matrix} ] ,

where X′_{ik} := X′(e_{i}, e_{k}) and e_{i} is the i–th unit vector in ℝ^{m_{l}}. Elements in the left lower part which can be obtained by symmetry are replaced by dots.

The blocks x_{α^{j}} ∈ L^{2}(ℝ^{τ(α_{1}^{j})}, ℝ^{τ(α_{2}^{j})}; ℝ^{m}), α^{j} = (α_{1}^{j}, α_{2}^{j}) ∈ {1, …, p}^{2}, are defined through

x_{α^{j}}(η_{1}, η_{2}) := X′(E_{α_{1}^{j}} · η_{1}, E_{α_{2}^{j}} · η_{2}), j = 1, . . , M := p^{2}, η_{1} ∈ ℝ^{τ(α_{1}^{j})}, η_{2} ∈ ℝ^{τ(α_{2}^{j})},

where the i–th column of E_{k} ∈ ℝ^{m_{l},τ(k)} is the [ ∑_{μ=1}^{k−1} τ(μ) + i]–th unit vector in ℝ^{m_{l}} (i ∈ {1, τ(k)}). Similarly, we construct y_{α^{j}} by partitioning Y′ and c_{α^{j}} by partitioning C′_{2}:

y_{α^{j}}(η_{1}, η_{2}) = Y′(E_{α_{1}^{j}} · η_{1}, E_{α_{2}^{j}} · η_{2}) = ψ^{T} D^{2}F^{0}(φ UE_{α_{1}^{j}} · η_{1}, φ UE_{α_{2}^{j}} · η_{2}),
c_{α^{j}}(η_{1}, η_{2}) := C′_{2}(E_{α_{1}^{j}} · η_{1}, E_{α_{2}^{j}} · η_{2}) = D^{2}F^{0}(φ UE_{α_{1}^{j}} · η_{1}, φ UE_{α_{2}^{j}} · η_{2}).

For k_{1}, k_{2} ∈ {1, …, p} we have A′_{l}E_{k_{1}} = ∑_{i=1}^{k_{1}} E_{i}a_{ik_{1}} and hence

X′(A′_{l}E_{k_{1}}η_{1}, E_{k_{2}}η_{2}) = ∑_{i=1}^{k_{1}} X′(E_{i}a_{ik_{1}}η_{1}, E_{k_{2}}η_{2}) = ∑_{i=1}^{k_{1}} x_{(i,k_{2})}(a_{ik_{1}}η_{1}, η_{2})

as well as

X′(E_{k_{1}}η_{1}, A′_{l}E_{k_{2}}η_{2}) = ∑_{i=1}^{k_{2}} X′ (E_{k_{1}}η_{1}, E_{i}a_{ik_{2}}η_{2}) = ∑_{i=1}^{k_{2}} x_{(k_{1},i)}(η_{1}, a_{ik_{2}}η_{2}).

Summarizing, we have transformed (4.20) into a system of equations, j = 1, …, M, given by

x_{α^{j}}(a_{α_{1}^{j}}_{α_{1}^{j}}η_{1}, η_{2}) + x_{α^{j}}(η_{1}, a_{α_{2}^{j}}_{α_{2}^{j}}η_{2}) − Ax_{α^{j}}(η_{1}, η_{2}) + φy_{α^{j}}(η_{1}, η_{2})

= c_{α^{j}}(η_{1}, η_{2}) − ∑_{i=1}^{α_{1}^{j}−1} x_{(i,α_{2}^{j})}(a_{iα_{1}^{j}}η_{1}, η_{2}) − ∑_{i=1}^{α_{2}^{j}−1} x_{(α_{1}^{j},i)}(η_{1}, a_{iα_{2}^{j}}η_{2}) (4.21)

ψ^{T} x_{α^{j}}(η_{1}, η_{2}) = 0.

If we now order the multi–indices α^{j}, j = 1, …, M, according to (3.20), then as in the complex case in section 4.2 the terms x_{β} on the right–hand side of (4.21) satisfy β ≺ α^{j} and the c_{α^{j}} can be computed from known quantities.

By the symmetry of X′ we have x_{(i,k)} (η_{1}, η_{2}) = x_{(k,i)}(η_{2}, η_{1}), and the symmetry of Y′, C′_{2} yields the same property for the y_{α^{j}}, c_{α^{j}}. Therefore, we can omit the calculation of those pairs (x_{α^{j}}, y_{α^{j}}) for which α^{j} is a permutation of some β ≺ α^{j}. Moreover, the symmetry of the forms x_{(i,i)}, y_{(i,i)} and c_{(i,i)} on the block diagonal will be used in case 4 below.

We notice that the equation (4.21) has the same structure as the original bilinear Sylvester equation (4.18), but now η_{i} is in ℝ or ℝ^{2} and hence the total size of (4.21) can be estimated by m · τ(α_{1}^{j}) · τ(α_{2}^{j}) ≤ 4m.

For the calculation of x_{α^{j}} we transform (4.21) into an equivalent square system with a bordered matrix which will then be solved by BEMW.

In the following we consider the cases (τ(α_{1}^{j}), τ(α_{2}^{j})) = (1, 1), (1, 2), (2, 1), (2, 2), where, in addition, we have to distinguish between α_{1}^{j} = α_{2}^{j} and α_{1}^{j} ≠ α_{2}^{j} in the case τ(α_{1}^{j}) = τ(α_{2}^{j}) = 2. The cases are easily recognized in our example above. We abbreviate the right–hand side of (4.21) as r_{α^{j}} and put e_{1} := \binom{1}{0}, e_{2} := \binom{0}{1}.

Case 1: τ(α_{1}^{j}) = τ(α_{2}^{j}) = 1. Equation (4.21) is equivalent to the bordered system

( \begin{matrix}[(a_{α_{1}^{j}}_{α_{1}^{j}} + a_{α_{2}^{j}}_{α_{2}^{j}})I − A] & φ \\ ψ^{T} & 0 \end{matrix} ) · ( \begin{matrix}x_{α^{j}} \\ y_{α^{j}} \end{matrix} ) = ( \begin{matrix}r_{α^{j}} \\ 0 \end{matrix} ) .  (4.22)

The solution fulfils y_{α^{j}} = ψ^{T} D^{2}F^{0}(φ UE_{α_{1}^{j}}, φ UE_{α_{2}^{j}}).

Case 2: τ(α_{1}^{j}) = 2, τ(α_{2}^{j}) = 1. Now we abbreviate x_{i} := x_{α^{j}}(e_{i}, 1), y_{i} := y_{α^{j}}(e_{i}, 1), i = 1, 2, and extract u = ({u_{11} \atop u_{21}} {u_{12} \atop u_{22}} ) := a_{α_{1}^{j}}_{α_{1}^{j}}, v := a_{α_{2}^{j}}_{α_{2}^{j}} from (4.19). In (4.21) we successively insert (e_{1}, 1) and (e_{2}, 1) and find the following systems

( \begin{matrix}[(u_{11} + v)I − A] & u_{21}I & φ \\ ψ^{T} & 0 & 0 \end{matrix} ) · ( \begin{matrix}x_{1} \\ x_{2} \\ y_{1} \end{matrix} ) = ( \begin{matrix}r_{α^{j}}(e_{1}, 1) \\ 0 \end{matrix} )  (4.23)

and

( \begin{matrix}u_{12}I & [(u_{22} + v)I − A] & φ \\ 0 & ψ^{T} & 0 \end{matrix} ) · ( \begin{matrix}x_{1} \\ x_{2} \\ y_{2} \end{matrix} ) = ( \begin{matrix}r_{α^{j}}(e_{2}, 1) \\ 0 \end{matrix} ) .  (4.24)

This may be written as a regular square system of order 2m + 2m_{l}:

( \begin{matrix}[(u_{11} + v)I − A] & u_{21}I & φ & 0 \\ u_{12}I & [(u_{22} + v)I − A] & 0 & φ \\ ψ^{T} & 0 & 0 & 0 \\ 0 & ψ^{T} & 0 & 0 \end{matrix} ) · ( \begin{matrix}x_{1} \\ x_{2} \\ y_{1} \\ y_{2} \end{matrix} ) = ( \begin{matrix}r_{α^{j}}(e_{1}, 1) \\ r_{α^{j}}(e_{2}, 1) \\ 0 \\ 0 \end{matrix} )  (4.25)

Case 3: τ(α_{1}^{j}) = 1, τ(α_{2}^{j}) = 2. Analogous to case 2.

Case 4: τ(α_{1}^{j}) = τ(α_{2}^{j}) = 2 and α_{1}^{j} = α_{2}^{j}. Define x_{ik} := x_{α^{j}}(e_{i}, e_{k}), y_{ik} := y_{α^{j}}(e_{i}, e_{k}), i, k = 1, 2, and u = ({u_{11} \atop u_{21}} {u_{12} \atop u_{22}} ) := a_{α_{1}^{j}}_{α_{1}^{j}} = a_{α_{2}^{j}}_{α_{2}^{j}}. From the symmetry of X′, Y′ and α_{1}^{j} = α_{2}^{j} we find x_{12} = x_{21}, y_{12} = y_{21} and r_{α^{j}}(e_{1}, e_{2}) = r_{α^{j}}(e_{2}, e_{1}). In (4.21) we successively insert (η_{1}, η_{2}) = (e_{i}, e_{k}) with (i, k) = (1, 1), (1, 2), (2, 2) and end up with a bordered system of dimension 3(m + m_{l}) given by

[\begin{matrix}[2u_{11}I − A] & 2u_{21}I & 0 & φ & 0 & 0 \\ u_{12}I & [(u_{11} + u_{22})I − A] & u_{21}I & 0 & φ & 0 \\ 0 & 2u_{12}I & [2u_{22}I − A] & 0 & 0 & φ \\ ψ^{T} & 0 & 0 & 0 & 0 & 0 \\ 0 & ψ^{T} & 0 & 0 & 0 & 0 \\ 0 & 0 & ψ^{T} & 0 & 0 & 0 \end{matrix} ] · [ \begin{matrix}x_{11} \\ x_{12} \\ x_{22} \\ y_{1} \\ y_{2} \\ y_{3} \end{matrix} ] = [ \begin{matrix}r_{α^{j}}(e_{1}, e_{1}) \\ r_{α^{j}}(e_{1}, e_{2}) \\ r_{α^{j}}(e_{2}, e_{2}) \\ 0 \\ 0 \\ 0 \end{matrix} ]

(4.26)

Case 5: τ(α_{1}^{j}) = τ(α_{2}^{j}) = 2 and α_{1}^{j} ≠ α_{2}^{j}. We proceed as in case 4 but now the matrices u := a_{α_{1}^{j}}_{α_{1}^{j}} and v := a_{α_{2}^{j}}_{α_{2}^{j}} are different. In contrast to case 4 we no longer have x_{12} = x_{21} and hence we find a regular bordered system of dimension 4(m + m_{l}) with the left upper block of the form

[\begin{matrix}[(u_{11} + v_{11})I − A] & u_{21}I & v_{21}I & 0 \\ u_{12}I & [(u_{22} + v_{11})I − A] & 0 & v_{21}I \\ v_{12} & 0 & [(u_{11} + v_{22})I − A] & u_{21}I \\ 0 & v_{12}I & u_{12}I & [(u_{22} + v_{22})I − A] \end{matrix}] (4.27)

In all five cases the upper left block of the square bordered matrices can be singular or almost singular for example if A has a zero eigenvalue in σ_{l}.

Another observation is that, if complex eigenvalues occur, but a real computation is required then it is of course not sufficient to have a linear solver for real translates A − λI, λ ∈ ℝ. Rather, certain bordered block matrices as in (4.25)–(4.27) appear, where A is on the diagonal and multiples of the identity form the off-diagonal elements. For a direct sparse linear solver (as we will use it in the example below) this structure is still acceptable. In this case one can also use the solver on (4.25)–(4.27) directly without recourse to BEMW.

After the computation of the coefficients x_{α^{j}}, j = 1, …, M, the bilinear form X′ can be constructed from these via

X′(e_{i_{1}}, e_{i_{2}}) = x_{(k_{1},k_{2})}(ϵ_{1,j_{1}}, ϵ_{2,j_{2}}) for i_{1}, i_{2} = 1, . . , m_{l}, e_{i} the i–th unit vector in ℝ^{m_{l}},

where, for ν = 1, 2, the integers k_{ν} ∈ {1, …, p} as well as j_{ν} ∈ {1, τ(k_{ν})} are uniquely determined by

i_{ν} = ∑_{μ=1}^{k_{ν}−1} τ(μ) + j_{ν}

and ϵ_{ν,j_{ν}} is the j_{ν}–th unit vector in ℝ^{τ(k_{ν})}. Now, analogous to (4.16), we transform back to obtain the n–th derivative in the original coordinates

X(h_{1}, …, h_{n}) := X′(U^{T} h_{1}, …, U^{T} h_{n})  for h_{1}, …, h_{n} ∈ ℝ^{m_{l}}.  (4.28)

Clearly, we can avoid this final transformation (4.28), if φ and ψ are chosen in such a way that Â_{l} in (4.1) already has the real Schur form (4.19).

5 Numerical Results

As a typical example we consider a spatial discretization of the well–known Brusselator ([34], [7], [26], [31])

\begin{aligned}u_{t} & = & D_{1}Δu + λ (A − (B + 1)u + u^{2}v) \\ v_{t} & = & D_{2}Δv + λ (Bu − u^{2}v) \end{aligned}  \begin{aligned}t & ≥ & 0 \\ x & ∈ & X \end{aligned}  (5.1)

first in one then in two space dimensions (X = [0, 1] and X = [0, 1]^{2}), in both cases subject to Dirichlet boundary conditions

u = A, v = \frac{B}{A} at x ∈ ∂X.  (5.2)

A, B, D_{1}, D_{2}, λ are positive parameters and \sqrt{λ} is the original length scale in the underlying trimolecular reactions. With the boundary conditions (5.2), equation (5.1) admits for all λ > 0 a single uniform steady–state solution:

u_{0} = A, v_{0} = \frac{B}{A} .  (5.3)

In the one dimensional example at parameter values A = 2, B = 6, D_{1} = 2D_{2}, D_{2} = 10^{−3} we discretize the second derivatives u_{xx}, v_{xx} by the standard three point difference method with Δx = \frac{1}{200}. After eliminating the boundary conditions, the resulting line system is of dimension 398. The eigenvalues of the Jacobian can still be computed explicitly and at

λ_{0} = 240 ∗ (1 − cos (π/200)) ≈ 0.0296082

a Hopf bifurcation occurs.

For the values λ = 0.01, λ = λ_{0} and λ = 0.08 we compute the quadratic approximation of the invariant manifold which belongs to the leading complex conjugate pair of eigenvalues. It is a weakly stable manifold for λ = 0.01, a center manifold at λ = λ_{0} and the unstable manifold for λ_{0} < λ < \overline{λ} = 0.11842. Notice that this is covered by case 4 in section 4.4. The computations were done by combining Govaerts’ solver [14] with a sparse linear solver from the package MATLAB, which applies a nonsymmetric minimum degree ordering algorithm.

Fig. 2 Fig. 3 Fig. 4 Fig. 5

In Figures 2–4 we show in the same scale the linear and quadratic approximation for these three parameter values, while Figure 5 contains all quadratic approximations in one picture. In each figure we plot the first component of both

φ s and φ s + X(s, s), s ∈ [−1, 1]^{2},

where φ = (φ_{1}, φ_{2}) spans the linear invariant subspace (with φ_{1}, φ_{2} normalized with respect to the L_{2}–norm) and X is the computed bilinear form.

Table 1 shows the growth of the number of flops as given by MATLAB with increasing dimension for the case λ = λ_{0}. In each column we show the factor by which the number of flops increases from the previous column. For these computations we use the “tight” settings of the minimum degree ordering parameters as defined by the MATLAB function spparms, which may lead to orderings with less fill-in.

Table 1: One dimensional case: Growth of the number of flops

dim(line sys.) = m 48  98  198  398  798 1598 3198

dim(system (4.26)) 150 300 600 1200 2400 4800 9600

factor  1 2.055 2.008 2.013 2.006 1.999 1.999

The left upper 3m×3m block of the bordered matrix in (4.26) has the bandwidth 2m + 1 which grows linearly in m. From Table 1 it can be seen that in the current example the quadratic approximations can still be computed with an amount of flops which grows only linearly in m.

In the two–dimensional example we compute in an analogous way the quadratic approximations of the invariant manifolds belonging to the leading complex conjugate pair of eigenvalues using the same parameter values as above and discretizing the Laplace operator by the standard five point difference method. It turned out here, however, that the MATLAB solver for sparse matrices, which is used as the black box solver in Govaerts’ algorithm BEMW, is efficient when combined with a reverse Cuthill-McKee ordering to reduce bandwidth. Table 2 shows the growth of the amount of flops and is arranged like Table 1.

Table 2: Two dimensional case: Growth of the number of flops

dim(line sys.) = m 72 162  338  648 1250 2450

dim(system (4.26)) 222 492 1020 1950 3756 7356

factor  1 3.729 4.012 2.938 3.423 3.603

Since the Jacobian has an upper and lower bandwidth of \sqrt{2m} in this example, the factors 3 up to 4 lie within the expected range. These factors certainly reduce further if some general iterative solver or a multigrid method is combined with BEMW.

References

[1] Abraham, R., Marsden, J. E., Ratiu, T. (1988): Manifolds, Tensor Analysis, and Applications. 2. ed., Appl. Math. Sci. 75, Springer, Berlin, New York.

[2] Abraham, R., Robbin, J. (1967): Transversal Mappings and Flows. W. A. Benjamin, New York, Amsterdam.

[3] Aulbach, B., Garay, B. M. (1993): Linearizing the expanding part of noninvertible mappings. Z. Angew. Math. Phys. 44, 469–494.

[4] Bartels, R. H., Stewart, G. W. (1972): Algorithm 432: Solution of the matrix equation AX + XB = C. Commun. ACM 15, 820–826.

[5] Beyn, W.–J. (1991): Numerical methods for dynamical systems. In: W. Light, ed., Nonlinear Partial Differential Equations and Dynamical Systems. Advances in Numerical Analysis 1, pp. 175–236. Clarendon, Oxford.

[6] Brunovský, P., Fiedler, B. (1986): Numbers of zeros on invariant manifolds in reaction– diffusion equations. Nonlinear Anal., Theory Methods Appl. 10, 179–193.

[7] Carr, J. (1981): Applications of Centre Manifold Theory. Appl. Math. Sci. 35, Springer, Berlin, New York.

[8] Chan, T. F. (1984a): Deflation techniques and block–elimination algorithms for solving bordered singular systems. SIAM J. Sci. Stat. Comput. 5, 121–134.

[9] Chan, T. F. (1984b): Techniques for large sparse systems arising from continuation methods. In: T. Küpper, H. Mittelmann, H. Weber, eds., Numerical Methods for Bifurcation Problems. ISNM 70, pp. 116-128. Birkhäuser, Basel, Boston.

[10] Golub, G. H., Nash, S., Van Loan, C. (1979): A Hessenberg–Schur method for the problem AX + XB = C. IEEE Trans. Autom. Control 24, 909–913.

[11] Golub, G. H., Van Loan, C. F. (1989): Matrix Computations. 2. ed., The Johns Hopkins University Press, Baltimore, Maryland.

[12] Govaerts, W. (1991): Stable solvers and block elimination for bordered systems. SIAM J. Matrix Anal. Appl. 12, 469–483.

[13] Govaerts, W. (1995): Bordered matrices and singularities of large nonlinear systems. Int. J. Bifurcation Chaos Appl. Sci. Eng. 5, 243–250.

[14] Govaerts, W., Pryce, J. D. (1993): Mixed block elimination for linear systems with wider borders. IMA J. Numer. Anal. 13, 161–180.

[15] Greub, W. (1978): Multilinear Algebra. 2. ed., Springer, Berlin, New York.

[16] Guckenheimer, J., Holmes, P. (1990): Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. 3. print., rev. and corr., Appl. Math. Sci. 42, Springer, Berlin, New York.

[17] Hartman, P. (1982): Ordinary Differential Equations. Repr. of the 2. ed., Birkhäuser, Boston, Basel.

[18] Higham, N. J. (1996): Accuracy and Stability of Numerical Algorithms. SIAM, Philadelphia.

[19] Hilger, S. (1992): Smoothness of invariant manifolds. J. Funct. Anal. 106, 95–129.

[20] Hirsch, M. W., Pugh, C. C., Shub, M. (1977): Invariant Manifolds. Lect. Notes Math. 583, Springer, Berlin, New York.

[21] Iooss, G., Adelmeyer, M. (1992): Topics in Bifurcation Theory and Applications. Advanced Series in Nonlinear Dynamics 3, World Scientific, Singapore, London.

[22] Irwin, M. C. (1972): On the smoothness of the composition map. Quart. J. Math., Oxford II. Ser. 23, 113–133.

[23] Irwin, M. C. (1980): Smooth Dynamical Systems. Academic Press, London, New York.

[24] Keller, H. B. (1983): The bordering algorithm and path following near singular points of higher nullity. SIAM J. Sci. Stat. Comput. 4, 573–582.

[25] Kleß, W. (1997): Numerische Approximation lokal invarianter Mannigfaltigkeiten in großen dynamischen Systemen. PhD thesis, University of Bielefeld.

[26] Kuznetsov, Y. A. (1995): Elements of Applied Bifurcation Theory. Appl. Math. Sci. 112, Springer, Berlin, New York.

[27] Lancaster, P., Tismenetsky, M. (1985): The Theory of Matrices. 2. ed., Academic Press, Orlando, London, New York.

[28] Marcus, M. (1973): Finite Dimensional Multilinear Algebra, Part I. Pure and Appl. Math. 23, Marcel Dekker, New York.

[29] Mielke, A. (1991): Hamiltonian and Lagrangian Flows on Center Manifolds with Applications to Elliptic Variational Problems. Lect. Notes Math. 1489, Springer, Berlin, New York.

[30] Moore, G. (1987): Some remarks on the deflated block elimination method. In: T. Küpper, R. Seydel, H. Troger, eds., Bifurcation: Analysis, Algorithms, Applications. ISNM 79, pp. 222–234. Birkhäuser, Basel, Boston.

[31] Murray, J. D. (1993): Mathematical Biology. 2., corr. ed., Biomathematics 19, Springer, Berlin, New York.

[32] Nelson, E. (1969), Topics in Dynamics I: Flows. Mathematical Notes, Princeton University Press, Princeton, New Jersey.

[33] Neudecker, H. (1969): Some theorems on matrix differentiation with special reference to Kronecker matrix products. J. Am. Stat. Assoc. 64, 953–963.

[34] Nicolis, G., Prigogine, I. (1977): Self–Organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations. John Wiley & Sons, New York.

[35] Rheinboldt, W. C. (1986): Numerical Analysis of Parametrized Nonlinear Equations. Univ. Arkansas Lect. Notes Math. Sci. 7, John Wiley & Sons, New York.

[36] Rosenblum, M. (1956): On the operator equation BX − XA = Q. Duke Math. J. 23, 263–269.

[37] Sell, G. R. (1985): Smooth linearization near a fixed point. Am. J. Math. 107, 1035– 1091.

[38] Shub, M. (1987): Global Stability of Dynamical Systems. Translated by J. Christy, Springer, Berlin, New York. Original edition: Stabilité globale des systèmes dynamiques. Astérisque 56, Société Mathématique de France (1978).

[39] Stéphanos, C. (1900): Sur une extension du calcul des substitutions linéaires. J. Math. Pures Appl., V. Sér. 6, 73–128.

[40] Stewart, G. W., Sun, Ji–guang (1990): Matrix Perturbation Theory. Academic Press, San Diego, Boston, New York.

[41] Wiggins, S. (1990): Introduction to Applied Nonlinear Dynamical Systems and Chaos. Texts Appl. Math. 2, Springer, Berlin, New York.

ℂ  σ_{s}  ^{α}  ^{β} σ_{c} ^{γ}  ^{δ} σ_{u}_{Figure 1: σ(A) = σ_{s} ∪ σ_{c} ∪ σ_{u}}

31

Figure 2: Linear and quadratic approximation at λ = 0.01 Figure 3: Linear and quadratic approximation at λ = λ_{0} (Hopf–point)

Figure 4: Linear and quadratic approximation at λ = 0.08 Figure 5: Quadratic approximations at λ = 0.01, λ_{0}, 0.08