We investigate the use of Krylov subspace methods to solve linear, oscillatory ODEs. When we apply a Krylov subspace method to a properly formulated equation, we retain the asymptotic accuracy of the asymptotic expansion whilst converging to the exact solution. We will demonstrate the effectiveness of this method by computing Error and Mathieu functions

Olver, Sheehan

English

Oxford University Research Archive

GMRES for oscillatory matrix−valued differential equations

GMRES for oscillatory matrix-valued differential equations

Mathematical Institute Eprints Archive

Report no. 09/3
GMRES for oscillatory matrix-valued differential equations
Sheehan Olver
Oxford University Computing Laboratory, Wolfson Building, Parks Road, Oxford, UK
Sheehan.Olver@sjc.ox.ac.uk
Abstract We investigate the use of Krylov subspace methods to solve linear,
oscillatory ODEs. When we apply a Krylov subspace method to a properly
formulated equation, we retain the asymptotic accuracy of the asymptotic
expansion whilst converging to the exact solution. We will demonstrate the
effectiveness of this method by computing Error and Mathieu functions.
Oxford University Computing Laboratory
Numerical Analysis Group
Wolfson Building
Parks Road
Oxford, England OX1 3QD
September, 2009
1. Introduction
Our aim is to compute the fundamental solution to the differential equation
Y ′(t) = (B(t) + ωA(t))Y (t), t ∈ (a, b) (1.1)
where A and B are d×d matrix-valued functions and ω is large. Applications for such equa-
tions include the computation of special functions (such as Airy, Bessel, Hypergeometric and
Matthieu functions [1]) the time-independent Schro¨dinger equation and semi-discretizations
of the linear time-dependent Schro¨dinger equation.
When the eigenvalues of A are imaginary, the solutions to (1.1) become more oscillatory
as ω →∞. Thus traditional time-stepping methods are inefficient for large ω. The accuracy
of modified Magnus expansions [5] does not degenerate as ω increases when used to compute
(1.1). On the other hand, the approach we construct actually improves with accuracy as
ω increases, and at an arbitrarily high asymptotic order. Moreover, we solve the equation
globally, allowing us to compute over unbounded domains and to higher accuracy than a
time stepping approach can achieve.
The simplest form of (1.1) is when d = 1, in which case we obtain the solution exactly:
Y = exp
∫
(B + ωA) dt.
The nonhomogenous form of (1.1) is the Levin differential equation [6] (changing notation
to emphasize that these are scalar functions)
Ly = y′ + iωg′y = f. (1.2)
A particular solution is
y = e−iωg
∫
feiωg dt.
In other words, solving (1.2) allows us to compute oscillatory integrals:
∫ b
a
feiωg dt = y(b)eiωg(b) − y(a)eiωg(a).
Many methods have been developed in recent years for computing oscillatory integrals,
with a recent review in [4]. One particular approach is to apply the GMRES method [10]
directly to the differential equation [7, 8, 9], which we refer to as differential GMRES. By
reformulating (1.2) as a shifted linear system
Mu = (M0 + iω)u = f, (1.3)
where M0 represents a linear operator, we achieve an asymptotic error in residual of
O
(
ω−n−1
)
,
1
where n is the number of GMRES iterations. This is the same asymptotic order as an
asymptotic expansion, however, differential GMRES actually converges for fixed ω, subject
to a condition on the growth of f in the complex plane. In Section 2, we go into the details
of this approach.
The goal, then, is to generalize this approach to the higher dimensional case of (1.1).
This is accomplished by reforming the equation as a (matrix-valued) shifted linear system
(1.3). By doing so, we obtain a method which also simultaneously achieves high asymptotic
order and (based on numerical results) convergence.
Differential GMRES in its pure form requires taking derivatives (and integrals in the
higher dimensional case) of the functions involved. In the general case, this is impractical.
To avoid this, while we develop the framework in the infinite-dimensional setting, in practice
we represent non-oscillatory functions as Chebyshev polynomials. This could be handled
automatically and adaptive by the chebfun system [2], however, for concreteness and speed
we use fixed order Chebyshev polynomials in our examples. More precisely, we represent
functions by there values at Chebyshev–Lobatto points, and the FFT can be used to compute
derivatives and anti-derivatives.
2. Oscillatory integrals
GMRES [10] is an iterative algorithm originally developed for solving finite-dimensional
linear systems
Av = b for A ∈ Cd×d and b ∈ Cd.
The Krylov subspace is defined as
Kn[A, b] = span
{
b, . . . , An−1b
}
.
GMRES finds an element v ∈ Kn[A, b] such that the norm
‖Av − b‖2
is minimized. This is accomplished through Arnoldi iteration [11], which finds an orthonor-
mal matrix Q = (q1, . . . , qn) whose columns span Kn[A, b] and an upper Hessenberg (only
zeros below the first subdiagonal) matrix H ∈ C(n+1)×n such that
AQ = QH and q1 =
b
‖b‖2
.
Then c ∈ Cn is chosen to minimize the norm
‖Hc− ‖b‖ e1‖2 .
The GMRES approximation is now v = Qc.
2
In [7, 8], this was generalized for unbounded, infinite-dimensional operators such as the
differentiation operator D. In exactly the same manner as the finite-dimensional case, given
the linear operator L, function f and a semi-inner product 〈·, ·〉, Arnoldi iteration determines
a row-vector q = (q1, . . . , qn) whose columns are functions which span the Krylov subspace
Kn[L, f ] such that (where we use the convention Lq = (Lq1, . . . ,Lqn))
Lq = qH and q1 = f‖f‖ .
Similarly, differential GMRES finds a function v ∈ Kn[L, f ] that minimizes the seminorm
‖Lv − f‖ .
This is accomplished by finding c ∈ Cn that minimizes the finite-dimensional norm
‖Hc− ‖b‖ e1‖2 ,
thence v = qc.
Now consider the case of Arnoldi iteration applied to a shifted linear system of the form
(1.3). We will denote the Hessenberg matrix produced by Arnoldi iteration for a particular
value of ω as Hω. A fact known from the finite-dimensional case which is also true in the
infinite-dimensional case is that the orthonormal basis Q is independent of ω and
Hω = H0 + iωIn,n+1. [3]
In other words, we only need to compute the Arnoldi iteration for one choice of ω to determine
the GMRES approximation for all choices of ω. Furthermore, GMRES satisfies the property
that the error in residual is
‖Mv − f‖ = O
(
ω−n
)
. [8]
In our particular case we wish to solve the Levin differential equation (1.2). But L is
not in the form of a shifted linear operator. If we assume that g′ does not vanish, we can
trivially put it into the required form:
M = L 1
g′
= D 1
g′
+ iω.
We thus apply differential GMRES toM, f and a suitable inner product to obtain v. Then
y = vg′ and hence we approximate the oscillatory integral by
v(b)
g′(b)
eiωg(a) − v(a)
g′(a)
eiωg(a).
since we represent functions by their values, we use the standard dot product on the sample
vector at Chebyshev–Lobatto points.
3
2 4 6 8 10 12 14 n
10-13
10-10
10-7
10-4
0.1
10 20 30 40 50 n
10-13
10-10
10-7
10-4
0.1
Figure 1: The error in computing
√
pierfc
√−iω
2
√−iω using differential GMRES with 15 (left figure) and
50 (right figure) mapped Chebyshev–Lobatto points, for ω = 1 (plain), 10 (dotted), 100 (dashed)
and 1000 (thick).
If the integral does contain a stationary point, i.e., g′ vanishes in [a, b], modifying the
operator to take the form of a shifted linear operator is more complicated, and detailed in
[9].
Consider the integral
∫ ∞
1
eiωt
2
dt =
√
pierfc
√−iω
2
√−iω for ω > 0.
Since the interval is unbounded, we represent functions by rational Chebyshev series — i.e.,
in terms of the basis Tk(
t−2
t ) — or, more precisely, by the values they take at the mapped
Chebyshev–Lobatto points. In Figure 1, we compute the absolute error of our approximation
for several choices of ω. As can be seen, the rate of convergence as n → ∞ increases with
the frequency, and the number of mapped Chebyshev–Lobatto points required to achieve
machine precision accuracy decreases with the frequency.
3. Oscillatory differential equations
We will use the notation exp to denote the matrix exponential, though we only apply it
to diagonal matrices, where
exp diag (d1, . . . , dn) = diag (exp d1, . . . , exp dn).
We will use the indefinite integral notation to denote
∫
A(t) dt =
∫ t
a
A(t) dt.
We now consider the solution of (1.1). Motivated by the preceding section, our goal is to
transform this equation into a nonhomogenuous shifted linear system. We will only consider
4
the case where A is diagonalizable with distinct eigenvalues. Thus assume that there exists
a matrix-valued function V that is nonsingular for all t ∈ [a, b] and diagonal matrix-valued
function Λ = diag (λ1, . . . , λd) where λi 6= λj for any i and j, so that
AV = V Λ. (3.1)
The requirement that V is smooth and nonsingular prevents application of this expansion
to coalescing eigenvalues. This is similar to the case of stationary points for oscillatory
integrals, and likewise outside the scope of this discussion. Note also that V is not determined
uniquely, however, our approach works for any choice of V .
We apply the transformation Y = VW to obtain
(VW )′ = (B + ωA)W ⇔
V ′W + VW ′ = (BV + ωV Λ)W ⇔
W ′ = (H + ωΛ)Y for H = V −1BV − V −1V ′.
We now apply the transformation
W = (I + U)e
R
(diagH+ωΛ) dt.
Therefore
U ′ + (I + U)(diagH + ωΛ) = (H + ωΛ)(I + U).
We can rephrase this as
LU = F for LU = U ′ + UdiagH −HU + ω[U,Λ],
where [U,Λ] is the standard matrix commutator [U,Λ] = UΛ− ΛU and F = H − diagH.
Our goal now is to premultiply this operator by an inverse to the commutator opera-
tor. Because of our choice of transformations, we have ensured that F has zeros along the
diagonal. Thus we can utilize the following inverse to the commutator:
Definition 3.1 For Λ = diag (λ1, . . . , λn) with distinct entries and M with zeros along the
diagonal,
QM =

0 m12λ2−λ1 . . .
m1n
λn−λ1
m21
λ1−λ2 0
. . .
...
...
. . . . . .
m(n−1)n
λn−λn−1
mn1
λ1−λn . . .
mn(n−1)
λn−1−λn 0
 .
From inspection, it is clear then that [QF,Λ] = Q[F ]Λ− ΛQ[F ] = F. However, the term
LQF = (QF )′ + (QF )diagH −HQF + ωF
5
will not necessarily have zeros along the diagonal, Because of the term HQF . In other words,
we cannot generate the Krylov subspace for LQ and F . Fortunately, diagonal matrices lie in
the kernel of the commutator. Hence we use the following, alternative commutator inverse:
[·,Λ]−1U = QU +
∫
diag (HQU) dt.
Then (using D =
∫
diag (HQF ) dt)
MF = L[·,Λ]−1F = (QF )′ + diag (HQF )−HQF
+ (QF )diagH +DdiagH −HD + ωF.
Since QF has zeros along the diagonal, the first and fourth terms also have zeros along
the diagonal. The second term cancels the diagonal of the third term. Finally, since D is
diagonal
DdiagH −HD = (diagH −H)D
also has zeros on the diagonal. ThusM successfully maps the set of infinitely differentiable
matrix-valued functions with zeros along the diagonal to itself.
Without modification, we can now construct a differential GMRES method for M and
F , provided an appropriate semi-inner product is used. We will use the Frobenius inner
product, where the dot product of two functions remains the dot product of their values at
Chebyshev–Lobatto points. This returns a function G which satisfies
F ≈MG = L[·,Λ]−1G.
Therefore, Y ′ ≈ (B + ωA)Y for
Y = V (I + [·,Λ]−1G)e
R
(diagH+ωΛ) dt.
The fundamental solution is then
YF(t) = Y (t)Y (0)
−1.
4. Example: Mathieu functions
Consider the system
Y ′ = (Bα + ωA)Y,
Bα =
(
0 1
α 0
)
and A =
(
0 1
− cosh 2t 0
)
.
The eigenvalues of A are ±√cosh 2t, hence the solution is highly oscillatory. Indeed, as
t→∞, it becomes exponentially more oscillatory. The fundamental solution to this equation
6
1 2 3 4 5
-1.0
-0.5
0.5
1.0
0.5 1.0 1.5 2.0 2.5
-3
-2
-1
1
2
3
Figure 2: In the left graph, the first row of the differential GMRES approximation to the funda-
mental solution of the Mathieu equation with α = 0 and ω = 10. In the right graph, a comparison
of the (1, 1)-entry of the GMRES approximation (solid) to the real part of Mathematica’s
built-in routine (dotted).
can be written in terms of Mathieu functions [1] asC
(
α(1 + ω), 12ω(1 + ω), it
)
−S
(
α(1 + ω), 12ω(1 + ω), it
)
i
C′(α(1+ω), 12ω(1+ω),it)
1+ω −i
S′(α(1+ω), 12ω(1+ω),it)
1+ω
 .
We apply our approach to approximate this function with α = 0 and ω = 10 over the in-
terval (0, 5). As can be seen in Figure 2, our solution is equal to the built-in Mathematica
routine (after scaling to obtain the fundamental solution) for computing Mathieu functions
when t is small. As t increases, the Mathematica routine quickly explodes, whereas our
method remains nicely behaved. Furthermore, the true solution must be real, as is our
approximate solution, whilst the Mathematica routine grows a nonzero imaginary com-
ponent. Comparison with a numerical ODE solver with a very small step size reveals that
our solution is indeed the correct one. We omit a graph of the convergence rate for different
values of ω, which is similar in behaviour to Figure 1: the larger ω is, the faster the rate of
convergence.
References
[1] Abramowitz, M. and Stegun, I., Handbook of Mathematical Functions, National Bureau
of Standards Appl. Math. Series, #55, U.S. Govt. Printing Office, Washington,
D.C., 1970.
[2] Battles, Z. and Trefethen, L.N., An extension of Matlab to continuous functions and
operators, SIAM J. Sci. Comput. 25 (2004), 1743–1770.
7
[3] Datta, B.N. and Saad, Y., Arnoldi methods for large Sylvester-like observer matrix
equations and an associated algorithm for partial pole assignment, Linear Algebra
Appl 154–156 (1991), 225–244.
[4] Huybrechs, D. and Olver, S. , Highly oscillatory quadrature , in: Highly Oscillatory
Problems: Computation, Theory and Applications, ed. B. Engquist et al,
Cambridge University Press, Cambridge, 2009.
[5] Iserles, A. , On the method of Neumann series for highly oscillatory equations , Bit
Numerical Mathematics 44 (2004), 473–488.
[6] Levin, D., Procedures for computing one and two-dimensional integrals of functions
with rapid irregular oscillations, Math. Comp. 38 (1982), 531–538.
[7] Olver, S., GMRES for the differentiation operator, SIAM J. Numer. Anal., to appear.
[8] Olver, S., Shifted GMRES for oscillatory integrals, Numer. Math., to appear.
[9] Olver, S. , Fast, numerically stable computation of oscillatory integrals with stationary
points , preprint, Technical Report NA-09/2, Computing Laboratory, Oxford
University, 2009.
[10] Saad, Y. and Schultz, M., GMRES: A generalized minimal residual algorithm for
solving nonsymmetric linear systems , SIAM J. Sci. Stat. Comput. 7 (1986),
856–869.
[11] Trefethen, L. and Bau III, D., Numerical Linear Algebra , SIAM,1997.
8


An extension of Matlab to continuous functions and operators,

Arnoldi methods for large Sylvester-like observer matrix equations and an associated algorithm for partial pole assignment, Linear Algebra Appl 154–156

GMRES for the differentiation operator,

GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems ,

Handbook of Mathematical Functions,

Highly oscillatory quadrature , in:

Numerical Linear Algebra ,

numerically stable computation of oscillatory integrals with stationary points ,

On the method of Neumann series for highly oscillatory equations ,

Procedures for computing one and two-dimensional integrals of functions with rapid irregular oscillations,

Shifted GMRES for oscillatory integrals,

Oxford University Research Archive (ORA)

Sheehan Olver

GMRES for oscillatory matrix-valued differential equations

Abstract

Similar works

Full text

Available Versions

Oxford University Research Archive

Oxford University Research Archive

Mathematical Institute Eprints Archive

Oxford University Research Archive (ORA)

Oxford University Research Archive (ORA)