Riccati-transformations in connection with boundary value problems

Riccati-transformations in connection with boundary value

problems

Citation for published version (APA):

Loon, van, P. M. (1982). Riccati-transformations in connection with boundary value problems. (Eindhoven University of Technology : Dept of Mathematics : memorandum; Vol. 8213). Technische Hogeschool Eindhoven.

Document status and date: Published: 01/01/1982

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

•

EINDHOVEN UNIVERSITY OF TECHNOLOGY

Department of Mathematics and Computing Science

Memorandum 1982-13 August 1982

Riccati-transformations in connection with Boundary Value Problems

by P. van Loon

University of Technology

Department of Mathematics and Computing Science P.O. Box 513, Eindhoven

Internal note

Title: Riccati-transformations in connection with Boundary Value Problems

In deze interne notitie doe ik verslag van datgene waarmee ik mij het eerste half jaar van mijn promotie-onderzoek o.l.v. Prof. dr. G.W. Veltkamp heb beziggehouden.

Ret ligt in de bedoeling dit rapport te zijner tijd uit te breiden tot een volwaardig THE-rapport. Met name zullen de uitbreidingen te vinden zijn op het terre in van speciale toepassingen, zoals randwaardenproblemen startend in een singulier punt. Verder zal zeker ook de vraag van effi-ciency aan de orde komen.

Eindhoven, 14 juni 1982

I Riccati-transformations

1.1. For t E [0,1] let x(t) ,f(t), b E 1Rn and A(t) ,BO,B

1 E 1R nxn

and

f(t),A(t) E C[O,I].

The problem of interest is to find a numerical solution of the linear differential equation

(I . I a) dt x(t) d = A(t)x(t) + f(t) ,

subject to the boundary conditions

(l.lb)

If we express x(t) in terms of its initial values, viz.,

(1.2) x(t)

=

M(t)x(O) + get) ,

nXn n

(M(t) E 1R ,g(t) E 1R ) , then we obtain

(I. 3a) dt M(t) d = A(t)M(t)

(1 • 3b) M(O)

=

I n and d dt get)

=

A(t)g(t) + f(t) 0.4) g(O) = 0 .

(Any nonsingular matrix X satisfying (1.3a) is called a fundamental matrix

The solution of the boundary value problem (BVP) (1.1) can in this way be found from the solutions of the two linear initial value problems (IVP)

(1.3) and (1.4). Namely, from (l.lb) and (1.2) we get

(1. 5)

As direct consequence we have

Theorem 1.1 (Keller [2])

Let A(t),f(t) € CP[O,l]. Then (1.1) has a unique solution x(t) € cp+1[0,J]

if and only if B_O + BIM(I) is nonsingular.

o

Finding the solution of (1.1) by solving (1.2) - (1.5) is called the shooting

method. This method is analytically nice and simple. However, forward inte-gration is in many cases numerically unstable, due to the presence of both rapidly increasing (dominant) and non-increasing (dominated) solutions. Therefore a mUltiple shooting technique has been developed (see f.i. Stoer,

Bulirsch [5J). There the interval [0,1

J

is partioned in, say m, subintervals

(1. 6)

On each subinterval [t.,t.

IJ

(i

=

O,-,m-I) one solves the IVP's:

1. 1.+ (I. 7) and (1.8) M (t . , t .)

=

In 1. 1. A(t)M(t, t.) 1. A(t)g(t,t.) + f(t) 1. g(t.,t.) == 0 . 1. 1.

Now

(1. 9) x(t)

=

M(t,t#)x(t~) + g(t,t.)

1 1 1 (t E [t.,t. IJ) 1 1+

and therefore we have to solve the system

(1. 10) Hx == d , where -T T T _{xTCt »T} x == (x (0) ,x (t l)

,

...

,

m dT == (b ,g etl,tO) T T

,

...

,

g (t ,t 1» T _{m m-} T and BO B_J M(tl,t_O) -I _n 0 H

₌

M(t 2,t1) -I n 0 Met ,t 1) -I m m- n

Now the internal points t. I are chosen such that the non-increasing 1+

solutions of (1.7) are not yet too strongly dominated by the rapidly

increasing solutions at t == t. l'

1+

1.2. Anoth~r possibility is that of decoupling the dominant and dominated

solutions, a technique which is also used in combination with multiple shooting (cf. Matthey [3J).

Assumption 1.2

The matrix A(t) has on the whole interval [O,IJ k-separated eigenvalues.

i. e. k eigenvalues A. with Re A. > 0 and n - k eigenvalues L with

L L J

Re A. < O.

J

o

Remark 1.3

Note that this separation is only for ease of speaking centered around O.

If it is anywhere else we can find a simple scaling which transforms it to

O.

1.2.1. Suppose A(t) is block-uppertriangular and x(O) is known. Then

(1.11)

can be solved by first solving

(1. 12a)

and hereafter

(I. 12b) Xl (0) given.

Such a decoupling is only of numerical interest if the partitioning of A(t)

is such that cr(A_2Z(t» is in the lefthalfplane (~-) of the complex plane

and cr(A11(t» in the righthalfplane (~+). With assumption 1.2 we therefore

k n-k kxk

have to demand that Xl (t) E'. 1R , x

2 (t) € 1R • All (t) € 1R and

A

22(t)

E: 1R(n-k)x(n-k) •

Moreover. numerical methods for IVP's are most stable if the eigenvalues of the Jacobian have negative real parts. Because of this itts more natural to integrate (1.12b) backward, i.e. to express xt(t) in terms of x}(l) rather than x1(O).

(Notice that the solution of (1.11) can always uniquely be expressed in

terms of x1(I) and K

2(O), since the homogeneous system corresponding to

(1.11). subject to x

2(O)

=

0 and x,(O)

=

0 has only the trivial solution.)

Thus: If for the differential equations (1.11) x

2(O) and xI (1) are given

and cr(A₂₂(t» c ~- and cr(A

11(t» c ~+ for all t E [O.IJ, then they can be

solved ~n a numerical~stable way.

1.2.2. How can we transform the system (1.1), for which assumption 1.2 is sat

is-fied, such that for the resulting system

(1.13) _{1t yet)}

₌

_{A(t)y(t) + f(t)} the matrix ~

[:II(t)

~12(t)

r

k

A(t)

₌

is A21 (t) A_Z2(t) } n-k k n-k

(I. 14a) i) block-uppertriangular, i.e. A

21(t) - 0

and

An undoubtly useful tool hereby is the real form of the Schur-decomposition of A(t).

nXn

Let U

_o

E ~ be orthogonal and such that the matrix

(1.15) A (t)

o

=

U-I

_o

A(t)U

_O

is, at t = 0, both (quasi)-uppertriangular and its eigenvalues are in

cor-rect order, i.e.

cr(A~2(O»

c

~-

and

cr(A~l(O»

c

~+.

Such a U

_o

exists by

assumption 1.2 (see f.i. van Dooren [IJ).

Then with

o

-1

x (t) = U

_o

x(t) ~

we obtain for

xO

a system which at t

o

has the properties of (1.14~.

nXn I

Now let Vet) E ~ ,V(t) E C [0,1

J

and regular.

With the transformation

(1.16) yet) = V -) (t)x 0 (t)

(1.la) changes in

(1.17) dt yet) d

=

[-1 V (t)A (t)V(t) 0 + V (t)· dt V(t)Jy(t) -) d

1

+ V (t)U- 1 - )

_o

f(t) So

(I.18a)

and

Now we want to choose Vet) such that A(t) satisfies (1.14) and (for

nume-rical purposes) H(V(t»

D

I!V(t) 11 II

v-I

(t) II not too large.

However, the expression for

A

₂₁ (t) is rather complicated" Therefore we

restrict ourselves to Riccati-transformations. i.e. V(t)'s of the form

(1.19) Vet)

(Note that (by continuity), at least on a subinterval [O,t_l), this suffices

since, for instance with V(O)

=

I , (1.14) is satisfied at t

=

0).

n Since we now have (1. 20) - :tR21(t) +

A~I(t)

+

o

0 + A 22(t)R21 (t) - R21 (t)Al1(t)

o

- R2l(t)AI2(t)R21(t)

So (1.14a) is satisfied in some interval [O,t

l) if R21 (t) is a solution of the

•

(I.2l)

As initial value we may take RZI (0) = 0 (in that case (1.14b) is satisfied,

at least for t

=

0). although other choices are possible.

Remark 1.4

One can easily check that for any fundamental matrix X, with X11(t)

inver-tible, the matrix X_2I

(t)X~!(t)

satisfies (1.21). Hence, the solution of

(I.ZI) subject to RZ1(0)

=

0 ~s . equal to M₂₁(t)M-I _I1(t).

Remark 1.5

As we have seen in 1.2.1, the solution of

(1.22) ddt yet)

=

A(t)y(t) ,

for which Y2(0) and YI (I) are pven, is uniquely defined. Therefore, trans-formation of

(l .23)

to block-uppertriangular form is possible if and only if (I.Z3) has a unique solution. If (1.23) has a non-trivial solution, then we say it has a

critical length at t

=

I, (Notice that a critical length depends on the

choice of R

21 (0).

As soon as, for some given RZI (0), (1.23) has a critical length at

t

=

tl E (0,1), then R

21 (t) tends to infinity if t goes to tl and therefore

R21 (t) does not exist for t ~ tl'

This is the main argument used by criticasters of Riccati-transformations to reject them. In 1.3 we will discuss how this problem can be overcome.

Till then it is assumed that our choice of R₂₁ (0) is such that (1.23) has

no critical length in [O,IJ.

Using (1.21), (1.19) and (1.17) we obtain

~

[y

I

(t)]

= dt Y2(t) Remark 1. 6 For p

=

I or p 00 we have II Vet) II II V-I (t) II p p 2 ( 1 + II R2 I (t) II p )

o

Hence, if IIR₂₁ (t)lip is not too large then Vet) of (1.19) can be used.

0

1.2.3. Similarly to (1.3) we may write the solution of (1.24) as

(I. 25) yet)

=

M(t)y(O) + get)

with

(1.26) dt H(t) d -

₌

A(t)M(t)

-

MO

=

I _n

and

We readily see that M21 (t) :: 0 d ~ = _{A22 (t)M22 (t)} M₂₂(0) dt M22 (t)

,

= I _n-k (1. 28) d ~ = All (t)M 1I (t) MI t (0) dt MIl (t)

,

=

I _k d ...

=

A 1I (t)MI2(t) + AI2(t)M22(t) MI2 (0)

=

0 . dt M12 (t)

•

Thus (1. 29)

However, we want yI(t) be expressed in terms of yI(I) and Y2(0) instead of

yeO), so let

(I. 30)

Then Nil' NI2 and

hI

satisfy

d ... = All (t)N] I (t) _{NII ()}

=

Ik dt NIl (t)

_•

(1.31) d ... _{'" All (t)N 12 (t)} + A_I2(t)N₂₂(t) N 12(I)

=

0 dt _{NI2 (t)}

_•

d ~ dt hi (t)

=

All (t)h_l (t) + A I2(t)g2(t) +

f

1(t) , hi (I) = 0 .

•

Remark 1.7

From (1.29) and (1.30) we obtain

One observes that [N11(t)

I

N

12(t)] is nothing but the first k rows of the

....

fundamental matrix N(t) defined by

N(l)

i.e.

=

I n

Hence, if Y2(O) and Yl(l) are known, then yet) can (numerically stable) be

computed from

(1.32)

The vector

(Y~(l). y~(O»T

can be determined from the boundary conditions,

viz. ,

I~J

y(l)

=

b

•

(1.33)

T T T

So, to determine (y_I(I)'Y2(O» we need:

RZ1(O) (see 1.2.4)

R21 (I),

M

₂₂(1) and g2(1). which are computed by a forward sweep, i.e.

from the left to the right.

N

₁₁ (O).

N

12(O) and hI(O), which are computed by a backward sweep, i.e.

from the right to the left.

1.2.4. What is a good choice for R2I (D)?

1.2.4.1. That the choice of the initial value R

21(O) is of importance can be seen

from the simple 2-dimensional case, where

•

Then

For different starting values we get the following trajectories (b < 0):

c-a

b

o

figure 1.1.

Clearly, if RZJ(O) ~ c~a then RZ1(t) exists for all t ~

O.

On the other hand,

c-a

z

(c-a)R_Zl (t) - b RZI (t) -I In c-a _-

_~

_{R '(0)} c-a ZI

Note that the conditions c < 0 and a > 0 can be replaced by a > c, which

corresponds to remark 1.3.

1.2.4.Z. As already indicated we want to have

and

o

0 +

a

Since A (0) is (quasi)-uppertriangular and a(A]](O» c ~ and a(A_2Z(O» c ~

a possible choice is RZl(O)

=

0 .

If with this choice there is no critical length in [O,IJ, then (1.33)

reduces to the more simple form:

( [N11(O) N:

2 (O)

1

[ Ik

M

₂: (I)

l)

[y

1

(I) : (1. 34) _BJIa + _BIUa = 0 _n-k RZI (I ) Y2(O) =b - B U _{o 0}

ll(O)]-BU [_

0 ] o I a qz(l)

o

T

I.Z.4.3. Another possibility is for instance R₂₁ (0) = AIZ(O) • Then the separation

between aeAl! (0» and a(AZZ(O» is at least as good as the separation between

o

0

creAl! (0» and a(A₂₂(O». lVhy we have chosen for R

Z1 (O) = 0 will be explained

e

1.2.4.4. One often chooses an initial value for R21 that depends on the boundary

condition at t

=

O. Suppose we have separated boundary conditions, i.e.

(1. 35)

2

and let BO be invertible.

If we choose

then the matrix in the system (1.33) is block-uppertriangular, namely

[ Ik 0

] [Nil

(0)

N:

2_{(0)] [ Ik}

M

_{2: (I)}

J

(1. 36) _BOUO + _BIUO R21 (0) _{I n- k} 0 _n-k _ R21 (I)

[I 2

_2-

1

_ BI + _{B:R21 (I)} B} M22 (}) • B2 0 So 'j,,-(o·) c.~ b~ Ce>1.V\.p" ..

".J ...

d ':J~p J V\o."'V\.~~ _{. - - I ·} -(,61

(I\);:;LB;

-+

i~'R21(1')l

(b,,-

R::d~) ~:h(1))

.

Y2(0)

=

(B

_o

2)-1 b .t.

However, we are not sure that (1.14b) is satisfied, not even at t

=

O.

Moreover, the solution of (1.21) would now depend on the boundary condition (1.lb), which implies that if we have to solve (1.la) for different boundary conditions, then each time (1.21) has to be solved.

e

1.3. Till now on we have assumed that R2l (t) exists for all tE [O,IJ and

more-- +

and cr(AI1(t» c ~ for all t E [O,IJ. However,

this is not quite likely.

Example Let 1.8 t

[ cos

ACt)

=

.

-s~n t =

r

cos 2t L-sin 2t sin

:J [:

-I

OJ

[cos

sin t cos t -sin

2tJ .

2t -cos

Then the Riccati-differential equation (1.21) becomes

-sin

:J

cos 2 sin 2t(R₂₁ (t) - 1) - 2 cos 2t R Z1 (t) •

In figure 1.2 trajectories for some different initial values are sho and

in figure 1.3 the corresponding eigenvalues of

All

(t) (Note that

A

₂₂(t) =

-All

(t».

As one sees the eigenva1uestructure of the solutions corresponding to

R21 (0)=10 t R21 (0)=3 10 8 6

e

4

II

~

it

1'\ _II,

i,

2

1\

jl

1\

\ i '\, ! '\, "-0 ' : I'

I

"-,

_i

"

\ \ -2 \ \

e

-4 -6 R21 (0)=0 -8 R21 (0)=-10 -10 figure 1.2. \,

"

_\, "-'\ \,

On the dotted lines (-tan t and - cotan t)

d

we have that dtR21 (t)=O.

\, ~'\, -;. t

i

'\, _" I

I

\

\i

~I

I: R 21(0)=2.5

2 o~--~~----~---~~~~~--~---0.4 0.6 0.8 1.0 1.4 1.6 -)0 t R21 (0)=0 R (0)=-10 21 R₂₁ (0)=2.5 figure 1.3.

1.3.1. What can we do if such things happen?

The only possibility is to restart the whole process, which means that we actually use a multiple shooting technique.

Suppose at t

=

t. a restart has to be done. Then the first thing we have to

1

o

do is to transform ACt.) in such a way that (1.14) is satisfied at t = t .• To

1 1

Assume that (m-I) restarts are necessary. Then

is a partitioning of [O,IJ •

m-I i m-I

We may define the sequences of matrix {Ui}i=O and {A (t)}i=O by

(1.37)

where U. is orthogonal and such that Ai(t.) is (quasi)-uppertriangular

1 L

and

cr(A~2(ti»

c (- and

cr(A~I(ti»

c t+.

Define for t E [ti,ti+1J . (i = O,-,m-l)

(1. 38) i x (t) U:1f(t) 1 _ U:1x(t) 1 (x(t) is the solution of (1.1». Then at [t_i,t_i+

1

J

we obtain the i-th Riccati-differential equation

d i dt R21 (t) (1. 39)

o .

With ( 1 .40)

e

"

-(l.la) is transformedon[ti,t i+l] into (1.41) d

[yt(t)]

[A~:

(t)

A~

2

(t) ] [yt(t)]

+

[ft(t)]

dt := Y2(t) "'1 Y 2 (t) f~ (t) A 2Z(t) where - i All (t) := All(t) i i i + A_I2(t)R ZI (t) , "'i A22 (t)

=

Ai Z2(t) i i - RZI (t)A_I2(t) and "'1 AI2 (t) := Al2 (t) i

Now we get the following equalities. From (I.tb): (I .42) = b • i-I i From U. IX (t.) = x(t.)

=

U.x (t.) 1- 1 1 1 1 (i = I,-,m-I): (1.43)

=

U.y i (t.) 1 1 (i :: I. -. m-l) •

Since the systems (1.40) are decoup1ed we may say that for i

=

O,-,m-l

(1 .44)

... i --i .-i -i -i

where M

22(t), gZ(t), NI1(t), N1Z(t) and h1(t) are defined analoguous1y to

(1.28) and (1.31). From (1.44) we get (1.45)

-[~

i-I Y j (t₁) i-I Y 2 (ti )

a

and (1.46)

N,:(t

~~

1

i ) n-k

Combination of (1.42), (1.45) and (1.46) yields

~

=: b

and from (1.43), (1.45) and (1.46)

[ i-I

1

y (t.) l. i-I Y2 (t_{i _l )} (1.48) • N1

~l.

:<t

i )

1

n-k

In matrix-vector notation this is

~ (1.49) Hy

=

d . h ~ nmXnm d d ~ JR nm Wl. t H E JR an y , Co

,

...

,

,

...

,

~ and H =

[~I"O)

_{~:"O) ~}

ROUO 0 _n-k

[

~

_~,:J

[-I

N -I -l1(tl)

'1:"1)

J

U_o -U 1 R~I(tl) 0 _n-k

o

o

One observes at once this system is similar to the multiple shooting

system (1. 10).

Note: the dimension of the system (1.49) is n less than (1.10) since

o

m-l

YI(t_O) and Y₂ (t_m) are not determined.

We have the following

Theorem 1.9

The matrix

is a fundamental matrix corresponding to (1.51) Proof: d i i i -- x (t) = A (t)x (t) dt ( i 0, - ,m-l) . Let

J

(t) E lRnxn satisfy ( i .. 0, .. , m-l ) and [ I . 0 /(t)

₌

I

}iCt)

-R~

I

(~)

_n-k Then (cf. (I. 44»

[ytCt) ]

=

[N~

I :t)

N~2Ct)]

rytCti+I)]

Y2(t) ~~ y 2(t.) M 22(t) l ~ So i ~ xl (t) Qi(t) yl(t i+l ) i

=

~ x₂(t) y 2(ti) i i We know that xl (t) Y

1 (t), 'it E [ti,ti+l], and

i i i i i

,...,

Since NIl (t) is invertible we obtain

NJ: (til

""i

]-1

n-k

Hence,

From the fact that this relation holds for all solutions xi(t) of (1.51), it follows that

Since Qi(t.) is regular the theorem is a direct consequence of this.

D

~

With this the.orem we see even better the connection between mUltiple shooting and repeated Riccati-transformati9ns, since

NJ: (til ""i

]-1

n-k and U.Mi(t)

=

M(t,t.)U. , 1 ~ 1

•

•

where M(t,t.) is defined in (1.7). 1

-

~

Our hope, however, is that, since the matrices R_{21, M22 ,} NIl and NI2

are obtained by integration in stable directions, with Riccati-transformations a smaller number of internal points are necessary, such that the dimension of the system (1.49) is more then n less than (1.10). In the next chapter we will discuss under what conditions this will be so •

1.4. So far we have concealed one important disadvantage of (repeated)

Riccati-1.4.1.

transformations. For the calculations of the backward-sweep. we have to

know the values of R₂₁(t). Hence, we have to store these values, which is, in

general, not feasible.

1.4.1.1. One possibility to overcome this problem is by making the backward-sweep independant of the forward-sweep. This can be done by another Riccati-transformation. In that case we want to find a transformation such that for the resulting system the first k coordinates can be computed directly.

Let us concentrate on the interval [t.,t. IJ (i=O,-,m-I). Then

we

first bring

1 1+

Define

(1.52) and

on (quasi)-undertriangular form by a transformation

D.,

like in (1.37).

1

Ai(t) _- - -1 U. A(t)U.

-1 1

ii(t) _- U. - -1 x(t) 1

•

Then we look for a Riccati-transformation of the form

(1.53) l. where R I2(t) satisfies (1. 54) (1.55) where and

Supposing

R~2(t)

(and

R~l(t»

exists for all t E [ti,ti+1J

z~(t)

can be

determined as soon as

z~(ti+l)

is known.

Now we have the following relations:

(1. 56)

(1. 57) (1.58)

-i

Combination of (1.56) - (1.58) yields

[z1(t)]_

[ [ Ik -Ri 1Z(t)]U-i

-I]

[XI

(t)]

_i

(1. 59) =: W (t)x(t) • Y2 (t)

_[-R~I

_(t) -1 x 2(t) I _n-k ]U. _~ Let, analogously to (1.30), ,... "" (1.60) z;(t)

=

N;I(t)z~(ti+l)

+

h;(t)

(i=O,-,m-l) •

~

Using (1.59) the boundary condition (l.lb) becomes:

(1.61)

=

b •

By (1.60) and (1.44) this is equal to

(1. 62)

•

•

•

Theorem 1. 10

Wi(t) is invertible for all t E [t.,t.

IJ

1. 1.+ ( i =0. - , m-l) •

Proof:

We will show that Zi(t) = Wi(t)U. is invertible. Therefore it suffices that

1.

Let

for some x E

Then, with y get

Since

are the invariant subspaces corresponding to the eigenvalues of A(t) with

In the same way as is done in (J.43) some continuity relations at t = t. ~

(i = 1,-,m-l) can be derived.

The resulting system is similar to (1.48).

[

z

~

(t . + 1)

1

Thus: x(t) can be found by first determining i ~

Y2{t.)

the continuity relations at t. and then solving ~

~

from (1.62) and

(I. 63)

Notice that (1.63) has a unique solution by Theorem 1.10.

i

1.4.1.2. Notice that the assumption that R_I2(t) exists for all t E [ti,ti+IJ is not

in all cases satisfied.

~

-If we find a critical length for R_I2(t) at t E (ti,t

i+1), or

cr(A1I(t»

¢

~+,

then a restart has to be done.

i

-Moreover, in that case, we also have to know the value for R₂₁(t), which,

i

-in practice, means that we have to calculate R₂₁ (t) a second time.

A second disadvantage of this method, called (repeated) modified Riccati-transformation, is that two (quite independent) Riccati-differential equations have to be solved.

Together with the necessity of inversion of the matrix Wi(t) in (1.62) and (1.63) it seems to me that the repeated modified Riccati-transformation is too complicated if some restarts are necessary.

1.4.2. Probably the best alternative is a variation on the solutionmethod proposed by Scott [4J • We have, as already seen, the relation

(!. 64)

Moreover it's simple to evaluate that

(1. 65)

i i i

where 5₁₁, 5

12 and ql are the solutions of, resp.,

d i _{i ( i i i )}

with i

dt Sll(t) = -Sll(t) AIICt) - AI2(t)R₂₁ (t) 5

11 (ti) Ik '

(1. 66) _{dt SI2(t)} d i -Sll(t)Ai i _I2(t)M... i ₂₂(t) with _{SI2(t i )} i

=

0

and

d i _{i ( i ... i}

f!(t») with i

o •

dt ql (t) = _{-Sll(t) A12 (t)g2(t)} + _{ql (t i )} =

,

Then, as soon as

x~(t.)

is known, xi(t) (t

~

[t.,t. IJ) can be found by

~ ~ 1+

solving

(1.67)

[

~!(

t)

1

g2(t)

i i

However, the derivation of x (t) from x (t.) is, in general, not well-~

conditioned, unless the intervals [t.,t. I J are very small. Hence, we _{~ ~+}

rewrite (1.65) in such a way that xi(t) can be expressed in terms of

i i

•

This is quite easy, since from (1.65) we have

(1.68)

So, if

x~(ti)

and

x~(ti+l)

are known, then

x~(t)

can be found with (1.68)

and

x~(t)

with (1.64). The values for

x~(ti)

and

x~(ti+l)

are obtained from the mUltiple shooting system which we get from the boundary conditions

[

0 ] - B U

1 m-1 gr;-l (I )

and the continuity relations

[ Ik

~;~

(t

_i) ] [ i-l

1

[ gi -

~

(

t

i)

]1

x (t.) U. 1 [ i-I ~ + ~- i-I R21 (t i ) x2 (ti_l) U

i[[

~

S~:(ti+!)

_{] [Xi(ti+!)]}

_[q~(:i+!)ll

S11 (t_{i +1)}

+

Notice that with this method all the integrations are made 1n a forward-sweep, namely

[

A~~

(t)

-+ [ : 0

o

]

~ Ri( ) ₌ . Ri(t) dt t 0

A~Z(t)

(l. 71)

. [

A~

1

(t)

:]

_ Ri(t{:

A~2(t)]

.

- R1(t) R1(t) 0 0

[8:1

(t)

1 Ri(t) SIZ(t) where

₌

... i RZI (t) MZZ(t) and

~r(t)l

₌

_[[

:

0]

- Ri(t)

[0

A~:(t]

_[:1(t)

1

dt 1() _{gz t}

A~Z(t)

0 gZ(t) (l. 72)

[f~:t)l

.

[f~(t)l

+ R1(t) 0

This method is quite similar with the Riccati-transformation introduced in I. Z, since

and

Hence the associated multiple shooting systems «1.49) and (1.69) with (1.70»

1.5. Conclusions

With solving a linear boundary value problem by the Riccati-method we will mean the method obtained by (1.71). (1.72) and (1.49). Now we can say that the Riccati--method may be succesful if the following conditions are satisfied:

the matrix A(t) has k-separated eigenvalues only a small number of restarts are necessary.

About the separation we can say that the more the eigenvalues are separated, the better the Riccati--method will work in connection with, for instance, ordinary multiple shooting.

On which facts the number of restarts depends will be discussed in the next chapter.

References

[IJ Dooren, P. van, A generalized eigenvalue approach for solving

Riccati equations, SIAM J. Sci. Stat. Comp., 2, 1981.

[2J Keller, H.B., Numerical Solution of Two Point Boundary Value

Problems, CBMS Regional Conference Series in Applied

Mathe-matics, 24, SIAM, Philadelphia, Pa., 1976.

[3J Matthey, R.M.M. and Staarink. G.W.M., A method for solving

general linear boundary value problems, report 8029,

Mathematisch Instituut Katholieke Universiteit Nijmegen, 1980.

[4J Scott, M.R., Invariant imbedding and its applications ,to

ordi-nary differential equations: an introduction, Applied Mathematics and Computation, Addison-Wesley, 1973.

[5J Stoer, J. and Bulirsch R., Introduction to Numerical Analysis,

I I The algebraic Riccati-equation

As seen in chapter I we are interested in the behaviour of the solution of a Riccati-di££erentialequation

(2. I)

nxn.

where A(t) € JR ~s such that at t

=

0 the conditions (1.14) are satisfied.

Besides, it is still assumed that A(t) has k-separated eigenvalues on [O,IJ.

It will be shown that this solution depends highly on (real) solutions of the algebraic Riccati-equation

(2.2)

Theorem 2.1

Let pet) E JR(n-k)xk (t E [O,t

l), tl > 0). The subspace

LS in invariant subspace of A(t) if and only if pet) satisfies (2.2).

Proof:

Thus: only if: Suppose Then =

[Ik

]A11(t)'

p(t) and which ~s (2.2).

o

2.1. First we want to say something about existence of a solution of (2.2).

Therefore we use some theorems proven by Stewart [3J.

(n-k)xk (n-k)xk

Define T

t : :JR -+ :JR , t E: [0, t}) , by

(2.3)

(2.4)

Interlude

inf II TtCP) II if 0

l

_{cr(T t )}

II P 11=1

Let T(P)

=

PA - BP • An upperbound for sep(A,B) can be given by (Stewart

[3J) :

sep(A,B) ~ inf Icr(A) - cr(B)

I •

In some cases a lowerbound can also be found.

Define

A (A)

=

max{Re A I A E cr(A)}

max

Suppose A (A) < 0 and A (-B) < O.

max max

-\

Then for T we have

co (2.5) T -\ (P)

₌

-J

e -tBp e tAdt • 0 Proof:

-'" -\ T(T (P» "" -

_f

e -tap tA e A dt + 0 co 0

= -

[ e -tBpe tA ]

₌

P • 0

f

co B e-tBp _e tAd _t

o

Here we have used the relation tA

lim e

=

0 - A (A) < 0 ,

t~ max

which is for instance proven by van Loan [2J. Moreover, he shows that

(2.6)

where n is the order of A and Pn-l ~s a (n-I)-degree polynomial.

From (2.5) we have

00

(2.7) sep(A,B) ;:: liT-III-I;:: [

f

lIe-tBlilie tAli dt]-l •

o

In combination with (2.6) we see that the value sep(A,B) depends in the

first place on the separation of the eigenvalues of A and B (not

necessa-rily centered around 0) and in the second place (by p _n-I) on the

eigen-vector structure of each.

Theorem 2.2 (Stewart [3J) Let As long as K(t) I <

4' '

o

there exists a unique solution pet) of (2.2) satisfying

(2.8) pet) <

.riEl

2

- oCt)

J +

It -

4K(t)

Moreover. a(A(t» is the disjoint union

(2.9)

o

A direct consequence of this theorem is that pet) of (2.Z) is continuous

in t

=

0, since

yeO)

= 0, which implies P(O)

=

0, and the spectra of

A₁₁(0) and A

2Z(0) are well separated, so 8(0) > O.

Let Theorem 2.2 be valid at t

=

t. Then both

... A

II A (~)P(~)II s 2n(t~y(t)

lZ

oCt)

and

Since

we have that, at least for t sufficiently small,

and

and K(t)

Since 6(t) >

°

and yet) = 0 we can apply Theorem 2.2 on ACt). Hence,

there exists a pet) with pet) =

°

and

at least if t - t is sufficiently small. Because of the fact that

we have found

Theorem 2.3

Suppose Theorem 2.2 is valid for all t E [O,t

l). I f

[

Ik

P(t)+P(t)

for all t E [O,t}), then there exists a continuous solution of (2.2) on

[O,t) with PCO)::

o.

Definition 2.4

~(t) is the invariant subspace corresponding to the eigenvalues of ACt) with positive real parts at time t.

o

Corollary 2.5

I f A(t) E: lRnxn is in C[O, I ] and has k-separated eigenvalues in [0,1] then Mk(t) changes continuously in [0,1].

Proof:

As long as Theorem 2.3 is valid this is trivial.

Suppose for t

=

t

l , the conditions of Theorem 2.3 are not valid anymore.

Then bring A(t

l) by a real Schur-decomposition on (quasi)-uppertriangular

form and restart the whole process.

. nXn

Let C(t) E lR be a general matrix in C[O,I] and Nk(t) some k-dimensional

invariant subspace. Then, in general, we can't say that Nk(t) changes con-tinuously in [0,1] since there may be discontinuities at those points where one of the k eigenvalues corresponding to Nk(t) becomes equal to one of the other (n-k) eigenvalues.

o

Orthonormal basisvectors of Mk(t) are, for instance, given by the columns of

where pet) E lR(n-k)xk satisfies (2.2).

Corollary 2.6

For any U(t) =

r

U11 (t)] E lRnxk

-1

L

_U21 (t)

if U11(t) exists then pet)

such that R(U(t»

=

Mk(t) we have:

-I

=

U

Coro llary 2.7

Let U(t) be as in Corollary 2.6 with Ull(t) invertible. Then

and

Proof:

Since R(U(t» = Mk(t) we have

- kxk - +

where All (t) E lR and cr(A

ll (t» ca;. Now

Because A(t) hask-separated eigenvalues this corollary is proven by the

last part of Theorem 2.2.

o

From Theorem 2.2 and Corollaries 2.6 and 2.7 we see that for any U(t) as in Corollary 2.7 we get the same

(2.10) pet)

=

U -1

2l (t)Ul1 (t) ,

I f pet) exists on [O,t}) and II ddt PCt) II is not too large with respect to the separation of the eigenvalues of ACt), then we will see that this pet) is the unique (local) attractor belonging to the Riccati-differential equation (2.1).

At the end of this section we will show the connection between the solutions of the algebraic Riccati-equations corresponding to, resp., the original matrix A(t) of (1.1) and the transformed matrix AO(t) of (1.15).

Lemma 2.8

Let Al (t) E: lRnxn have k-separated eigenvalues for all t E: [0, I].

Define A 2(t)

-I

=

_{U A1(t)U, where U is such that (A2(O»21}

=

0 and

With P.(t) we mean the solution of the algebraic Riccati-equation corres-~

ponding to Ai(t), which generates the space (Mk)i(t) (i = 1,2).

Let T >

a

be such that P₂(t) exists for all t E: [O,T). Then for t < T we have:

as long as UtI + U

I2P2(t) is invertible we have the relation

Proof: Since

we have

So

is a k-dimensional invariant subspace of A1 (t). Moreover,

+

a«A₂(t))11 + (AZ(t))lZPZ(t) c ~ , which implies that

from which the lemma follows.

Z.Z. Notice that for U(t) E JRnxk with R(U(t))

=

Mk(t) we have U

Z1 (0)

=

0,

so U₁₁(0) is invertible and therefore U}I (t) is invertible for t E [0,t₁)

(t

l sufficiently small). What does it mean if U11(t) is not invertible

anymore, say at t

=

t}?

Obviously it means for instance that M

k(t1) can't be written as

Hence, there will be no continuous solution of (Z.Z) with P(O)

=

0 and

which does still exist at t = t

l. Of course, the other way around is also

true.

45

-From geometrically point of view we may even see better what happens.

Since U₂₁(0)

=

0 we have

Let t] be such that UI1(t) is invertible for all t E [O,t

l) and UII (tl)

is not invertible. Then there exist a

°

:f: a. E lRk such that

This means that Mk(t

l) and Mk(O)L have a nontrivial intersection. Now lR

n .I.

can't be written anymore as the union of Mk(t_l) and Mk(O) •

We know that a solution pet) of (2.2) with P(D)

=

D exists on [O,t_l).

Suppose lim II pet) II is finite. Then the matrix

t-+t 1

has full rank.

This is in contradiction with the nontrivial intersection of Mk(t

l) and

MkCO/. Hence, IIP(t)1I -+ <lO if t -+ t

46

-2.3. Reading 2.1 one may think that the transformation

yet)

= [

Ik -pet)

o

I n-k x(t)

is a good one for problem (1.1).

However, we then get

d ::: [AII(t)d+A12(t)P(t) dt yet) - dt pet) [ Ik -pet) I

a

]f(t) n-k

which is, in general, not a decoupled system.

That is the reason why we want the solution of (2.1) instead of (2.2).

Our hope is that R21 (t) of (2.1) has the same nice properties as pet),

which implies that R

21(t) must be close to pet). This is the main reason

why we have chosen R_2I(O)

=

O.

Define

(2. 11) E(t)

=

R₂₁(t) - pet) •

Obviously, we obtain that E(t) satisfies

d d

dt E(t)

= -

dt pet) + [A

22(t) - P(t)A12(t)JE(t)

Let us look at the behaviour of E(t).

First we concentrate on the more general problem d

dt E(t)

=

D₂(t)E(t) - E(t)D

1(t) + f(t)E(t)

(2.13)

E(O)

=

°

where E : [O,lJ + JR(n-k)xk , D2 : [O,IJ + JR(n-k)x(n-k),

D1 : [O,lJ + JRkxk and f[O,IJ x JR(n-k)xk + JR(n-k)xk, all continuous

on [O,l].

The solution of (2.13) satisfies

t

(2.14) E(t)

=

J

M

2(t,T)f(T,E(T»M1(t,T)dT,

°

where MI(t,T) and M₂(t,T) are solution of, resp.,

and

Suppose f has the following property:

(2. IS)

Let t be fixed.

Define the sequence of matrices

{Ei(t)~Jm

by

i=O

I k'

n-+

•

(2.16)

t

=

f

°

One observes that for all ~c [O,t] we have

~

II(Ei+l_Ei)(~)1I

:0;

f

IIM2(~'L)11

cCr) II (Ei_Ei-1) (L)II

IIM1(~,L)lldL.

With we obtain (2. 17)

°

~

g(~)

=

J

c(L)IIM2(~'L)1I IIM1(~'L)lIdL

°

II (Ei+I-Ei)(OIi :0; sup II (Ei_Ei-l)(iJ)11 •

g(~)

=>

O:O;iJ~~

i+l i

sup II(E -E) (~)II ~ ~c[O,t]

sup

~dO,t]

Let v. (t)

=

sup II (Ei+l_Ei)

(~)II

1. ~dO, t] (i

=

0,1,2, ••• ) and d(t)

=

sup g(~). ~dO, t] Lennna 2.9 sup C;dO,t]

If d(t) :0; 0 < I, then E(C;)

=

lim E1.(~) exists for all ~ c [O,t] and

Proof:

From (2.17) we observe that v.(t) $ v. I(t) • d(t). Hence, if d(t) $ 0 < 1,

1. 1.-then, 00

LV.

(t) i=O 1. and therefore i-I

lim

L

(Ek+I_Ek)(~) + EO(~)

i-+<x> k=O

exists and, since

EO(~)

=

0, its norm is bounded by vO(t)/l-o .

By continuity arguments we see that E(t) satisfies (2.14) and is therefore the solution of (2.13). Since ~ su p

J

II M2

(~

, • ) II IIf ( • , 0) II ~dO,

tJ

0

o

we have that the existence of E(t) and its norm depends on the first place on

the behaviour of the function IIM2(~,.)1I • IIMI (~,.)II and in the second place

on c(.) and II f(.,O)11 •

If we return to our specific problem (2.12) we get

(2.18) DI (t)

=

AII(t) + A_{I2(t)P 2I (t) ,}

d

The foregoing results can't be used directly since there is no uniform

c(t) like in (2.15). The best estimate we can give is

Before we can use Lemma 2.9 we have to find a uniform estimate of IIEi(t)1I t

defined by (2.16) and (2.18).

For E; E [0, tJ we have

~

II E i + 1 (E;)II

~

J

II M2 U;, T) II II MI (t;, T)II

o

i IIf(T,E (T»lIdT E;

~

sup II f

(~,Ei(~)

)11 ~€[O,O

f

IIM₂U;,T)11 IIMl (E;,T)lIdT •

With w. (t)

₌

l. and

d

(t) we obtain

o

sup II Ei(E;)11 E;E[O, tJ E;

sup

_f

IIM2(~,1")11 II HI (';,T)II dT

~E[O, t]

0

~

( sup II ddt; PU;)II +

•

Write

pet)

=

sup II d P(E;)II

E;€[O,t]

and

a

I2(t)

=

_1;€[O,t] sup II A12 (E;)II

,

then

(2.19) w. 1 (t) ~ d(t)(p(t) + a_I2(t)w2

i(t» 1+

Like in Theorem 2.2 we have that if d(t)2p (t)a₁₂(t) <

t.

then wi(t)

~

d(t)p(t)

for all i €

~,

which implies that under this condition sup II Ei(OIi is

uniformly bounded by d(t)p(t).

E;ECO, tJ

- 2 1

Hence, if d(t) p(t)a₁₂(t) <

4'

then we may take in (2.15) as

Lipschitz-function

(2.20)

where

Now we may use Lemma 2.9 to find that if

then

E;.

sup

J

IIM₂(E;.,T)11 11M} (E;.,'t)1I .c(T)dT :::; <5 < I ,

1;€[O,t]

a

E(I;)

=

lim Ei(O

i~

(2.21)

~

IIE(OIl ::;;

(SUp

J

IIM2(~.T)1I

11M)

(~'T}lIll'(-r)dTJ

/ )-6 !;E[O, tJ 0

where

d

1T(T) .. IIdT P(T)II "" II f (,r, 0)11

However, if we know that E(t) exists then sometimes a sharper estimate can be given.

Corresponding to II· II we may define the measure ~ of a matrix by

(2.22)

hA )..l (A) = lim II I + hA II - 1 .. lim In lie II

MO h MO h

(see Dahlquist [IJ)

Let (2.23) then (2.24) t a(t,.)"

J

()..l(D₂(S» + )..l(-D\(s»)ds, T

From (2.14) and (2.18) we get t

II E (t)1I

::;; J

-e

Look at the function e [O,T) ~ [0,00), some T > 0, defined by

t

(2.25) e(t)

₌

_f

o

Then e(t) is the solution of the Riccati-differential equation:

(2.26) where { ddt e(t)

=

e(O) ==

°

aCt (t) = at 2 aCt ~(t) + _{Ct 12 (t) • e (t)} + at(t)e(t)

If ~(D2(t» + ~(-Dl(t» > 0, then e(t) will grow fast. However, if

aCt 2

~(D2(t» + ~(-Dl(t» < 0 Vt ~ 0, then if 4Ct12(t)~(t)

s

(at

(t»~e(t)

will, in general, approach

pet)

a solution of the algebraic Riccati-equation, co('re~r--(nu::U\"\~

to

'-"2.'2-6) •

Hence, if this function P remains or becomes small then e(t) will be small,

which implies that II E(t)1I will be small too.

The best this is illustrated by

Example 2.10

Suppose that for some t we have A(t) == A(t) Vt ~ t. This means that

•

Hence, P(t) = 0 Vt ~ t.

From this it follows that i f E(t) is such that IIE(t)1I

area of pet) then we get lim E(t)

=

O.

t-+<>o

in the contraction

I f we had used (2.18) we at best had found that IIE(t)11 would stay bounded.

0

That the estimate (2.24) is not always realistic ~s shown by the following

example •

Example 2.11

Let

D2

=

-D

I

=

[-10

40] and II. II

=

II· II 2 •

o

-10 Then SO 20(,"-,) exp(a(~,T»'= e ~ • However.

J .

0

Similar to Theorem 2.2 we see that i f for all ~ E [O,t) a(D

2(0) and

a(DI(~» are separated such that ~(D2(~» + ~(DI(~» < 0 and

II

:~

P(Olt and: II Al2 (011 are small, relative to this separation, then

II E(t)1I will be small.

This is of importance since

and

Hence, as long as a(A11(t) + A

I2(t)P(t» and a(A22(t) - P(t)A I2 (t) are

well separated and IIE(t)11 and IIA

I2(t)11 are small, then the

Riccati-transformation decouples the system (1.1) in a stable way.

2.4. Example 2.9 Let [ COS t A(t) = -sin t sin t

o

-20 [ COS t sin t -sin t 20 0 cos t cos t

Then pet) of (2.2) with P(O) = 0 is equal to -tant.

•

•

t R 21(O)=10 TI/2 O~~---+--~---R (0)=-10 21 -+ t

Notice that this P(t), which is the same as in Example 1.7, is much more 'attractive' since the eigenvalues are more separated.

Although a critical length is already found close to

I'

we can come further

without a restart then all cases of Example 1.7. There we had to stop since

Since, ~n this example, cr(A11(t) + A_I2(t)P(t»

=

{20} for all t, and A

I2(0)

=

0 we have that the eigenvalues are in correct order for small

values of t. If t is not very near to 0 on

I

we have that E(t) is small

(in all cases with IR

21 (0)1 ~ 10 we had that ]R21 (t) - p(t)1 s 0.75

if t E [0.2,1.4]). So, the eigenvalues stay for a long time in correct

order.

~ 2.5. The conclusion of this chapter is the following: success of a

Riccati-transformation depends on four crucial factors:

)0 a sufficient separation of the eigenvalues with positive, resp.,

negative real parts,

and relative to this separation

20 pet) should not be too large, ~.e. the invariant subspace Mk(t)

should not contain vectors whose directions become (nearly) ortho-gonal to Mk(O) ,

30 ddt pet) should not be too large, which implies that the directions

of the basisvectors of ~(t) are only slowly varying,

40 A₁₂(t) should not be too large.

The main question we have not answered yet is: how strong are these ditions, i.e., are there problems, other than those where A(t) is a

con-stant matrix, for which these conditions are satisfi~d? We will return

•

•

References

[1] Dahlquist, G., Stability and error bounds in the numerical

integration of ODE's, Trans. Royal Inst. of Technology, No. 130, Stockholm, 1959.

[2J Loan, C. van, The sensitivity of the matrix exponential,

SIAM Journ. Numer. An., 14, 1977 •

[3] Stewart, G.W., Error bound for approximate invariant

sub-spaces of closed linear operators, SIAM Journ. Numer. An., 8,

III. The conditioning of a linear boundary value problem

(Most of the .results in this chapter can be found in Matthey [2]) •

. 3.1. A linear BVP is called well-conditioned if the solution of the problem is only (relatively) moderately perturbed by perturbations in the boundary condi tions.

As indicated earlier the BVP (1.1) is not always well-conditioned, even not

,.

if the problem has a unique solution, depending continuously on the data.

This is for instance the case if the boundary conditions are

n-k k

(b

2 € lR , b 1 € lR ) ,

and is very near to a critical length (see Remark 1.5).

Let

(3.1) dt x(t) d

=

A(t)x(t) + fCt)

subject to

•

BO x(O) + BI xCI)

=

b •

Suppose, instead of (3.1), we have the perturbed problem

(3.2) dt x(t) - A(t)x(t) d - - + fCt),

For the difference (3.3) ~x(t)

=

~(t) - x(t) one gets (3.4) dt d ~x(t) = A(t)~x(t) subject to

•

which has the solution

where M(.,.) is defined analogously to (1.3) by

(3.6) dt M(t,s) d

=

A(t)M(t,s) , M(s,s)

=

In Write (3.7a)

•

and (3.7b) Then (3.5) is identical to

~x(t)

=

[set)

+ oS(t) ]-1 (ob-oBox(O) - oB1x(1)}

3.2. Let II-II be any norm working on lRn. In connection with the solution x(t)

of the differential equation (3.1) we define in the space

A= lR n x lR nXn x lR nXn the norm II· IIA as follows:

if u

=

(a,A,B) € A, then lIul~ = II all + II All IIx(O) II + IIBII IIx(l)ll, where II A II and II B II are the induced matrix norms corresponding to II • 11 •

Let 8

=

(b,BO,B

I) and 08 = (ob,oBO,oB]).

Since ~x(t) is independent of the inhomogeneous term f(t) we concentrate

ourselves to homogeneous problems.

The relative conditionnumber eN of a linear BVP is now defined by

r

(3.9) := lim max

e:+0 lIo811("'e:

sup II ~x(t) II) II 08 IIA

(

s~p

IIx(t) II , /

lmiA .

t

In this definition we must have that sup II x(t) II 1= 0 and II aliA 1= o. This is

t

obviously true if we suppose that b 1=

o.

Theorem 3.1

The conditionnumber CN of a linear BVP is equal to

r

(3. ] 0)

II ell - sup II [BOM(O. t) +

Bl

M(1, t)

]-1

11

t

sup IIx(t) II t

Proof: Using (3.8) we have II S-I (t) II U 6b - OBOX(O) - oB 1x(1) II 1 - "OS(t)S-1 (t) II 115-1 (t)" II oS If A 1 - [I oS(t)S-l (t) "

• Notice that oS(t)S-1 (t) == [OB

O + OB1M(l ,O)][BO + BtM(1 ,0)]-1, which is independant of t. Hence, II oS(t)S-t (t) II == 0 (II oS II).

So,

-I

sup I[ S (t)" II 013 "A

sup lIc:'x(t) II ~ --'-t _ _ _ _ -:--_ _ t t - IIbS(t)S-l (t)11

Moreover. for' any ob with II ob II ~ E: there exist OBO and oB l ' such that

(3. II)

= lIobll + HBal! I1x(O) II + II oBI II IIx(l) II == E: •

Let t be such that sup II S-l (t) II = II S-1

(~)

II and u such that 1 A 1 A t

IIS- (t)ull == IIS- (t)1l lIuli •

Now choose oS such that (3.11) is satisfied and

-]

ob - oBox(O) - OB

lx(l)

=

[I + oS(t)S (t)]u •

•

IIS-I(~)11

lIeBII_A ~ 1 + II

eS(~)S-l (~)II

-1 sup II S (t) II II 1513

n

_A -I sup II S (t) II II oallA t 1 + II oS(t)S-l (t)11 --'-t _ _ _ _ -.,-_~ ::; sup IIllltl tlll . + lIeS(t)S-l(t) II t

Hence, for the absolute conditionnumber CN we have derived a (3. 12) CN a := lim max c:i-O II 01311 =8 sup 116x(t) II t = sup II S-I (t) II • t

For the relative conditionnumber we thus obtain

-1 sup II S (t) II • II 13 IIA t CN _r = -sup x(t) t

sup II [BOM(O, t) + B1M(I, t)

J-

111 II S IIA t

= -sup II x(t)1I

t

Corollary 3.2

'The initial value problem d

dt x(t)

=

A(t)x(t)

(3.13)

x(O) = b (b

1=

0)

is well-conditioned if and only

if

sup IIM(t,O)x(O)1I

t

II x(O)1I ~ sup IIM(t,O) II •

t

This means that the most dominating solutions must be represented in x(O).

o

'.

,

63

-Notice that in absolute sens an initial value problem is well-conditioned

if and only if cr(A(t» ~ ~-. (See Strom [3J)

3.3.

3.3.1. With (3.10) and (3.12) we have expressions for the conditionnumber of a linear BVP. However, these expressions are not very manageable. Some diffe-rent estimates are possible.

Using the fact that 11 Set) II f:: 0 and S(t)x(t) :: b for all t E [O,IJ we have

(3.14)

sup II S -1 (t) II t

CNr S; II aliA • UbI!

7

inf US(e}1I

t

II S IIA -1

=

iiblI .

sup II S (t)" inf II See) "

II ellA < -- II bll t

e

sup f{(S(t» • t

However, this upperbound is, in general, not sharp.

Another estimate is derived from the fact that for any fundamental matrix X(t) corresponding to (3.1) we have

•

This leads to the inequality

(3.16) CNa

~

II [BOX(O) + B 1 X(l )]-111 sup II X(t) II •

t

In Matthey [2J it is shown that the fundamental matrix M(t,O) of (3.6) will give, in general, not a sharp estimate.

We will investigate under what conditions a matrix X(t) will give a sharp estimate in (3.16).

Therefore we use that if A E JRnxn is invertible and B E JRnxm , then

liB II

=

II A-I AB II ::;; IIA-111I1ABII =I>

IIABII'~

IIBII / II A-I II •

Thus, for any fundamental matrix X(t) we have

II [BOX(O) + B] X(1 )]-1 II - - - - : - - -

~

IIX(t)[BOX(O) + B 1X(l)]-I 11 IIX-I(t)II II [BOX(O) + BI X(I )]-1 II inf II X-I (t) II t sup II X(t) II -1 ::;; II X ( t) II II [B OX ( 0 ) + B I (1 ) ] .' II ::;; sup 1\ X(t)[BOX(O) + B) X(1 )]-1 II t ::;; sup 11 X(t) II II [BOX(O) + Bl X(l )]-1 II • t

Since I > t we obtain, with C inf IIX-t (t) II - sup H(X(t»

t/sup H(X(t)} t t t C • sup t II X(t)1I II [BOX(O) + B_t X(1 )]-1 II (3.18) (I) (2) ::;; _CN ::;; _{sup II X(t) II II [BOX(O)} + B 1X(l)]-11i a t

•

However. since

sup H(M(t,O»

~

sup H(X(t» • H(x(O»

~

(sup H(X(t»)2

t t t

we have that C s (sup H(M(t.O»)-!, which is, in general, very small.

t

In the same way one derives the inequalities

(3.19)

sup IIX(t)IIIICBOX(O) + B1X(I)]-111

t

(1 ) (2)

sup IIX(t) III1CB

OX(O) + B1X(l)]-lU .

t

Hence, the estimate (3.16) is sharp if the conditionnumber of BOX(O) + BIX(I)

is small.

Because H(B_O + BIM(I,O» s H(BOX(O) + BIX(I»

H(X(D»

we see that if

M(t,O) gives a rough estimate and X(t) a realistic estimate, then H(X(O»

is large.

3.3.2. Suppose A(O) satisfies (1.14) and that one Riccati-transformation will

suffice through the whole interval [O,IJ. Now we will check if the matrix

(3.20)

Example 3.3 Let A(t)

= [:

J

O<aeJR. Then M ( t ,

0)

= [e ta

0]

and

°

e-ta __ [e(t-ol)a Q(t) So H(M(t,O» = e 2ta and

{

a-2ta _{, t} S; I e 2 H(Q(t)) = 2ta-a , t ~ 1 e 2 2 2a

Since (sup H(Q(t))

=

e = sup H(M(t,O» we see that this Q(t) gives in

t t

(3.16) the best lowerbound, but for large values of a the quantity C lW"l (~.\til) '15

-a

still very small (e

!).

Now suppose that BO = B

J = Then and '1-+

eQ.

I • n

o

-a

1

J+e

o

-a J+e

1 .

SO _H(BO + BIM(O)=~Te~ and H(BOQ(O) + B1Q(1») = 1. Hence, Q(t) is in

•

(Note that we may expect better results from (3.19) than from (3.18), since the lowerbound in (3.19) depends on the boundary conditions,

while in (3.19) it only depends on the fundamental matrix).

Suppose we have separated boundary conditions, i.e.,

Then

Assuming that R

21(1) and NI2(O) are of order I and the spectrum of A(t)

is well-separated then this matrix is of the form

l

-

order I small

small ] order 1

and therefore will have a small conditionnumber •

The best this is illustrated by

Example 3.4

If the boundary conditions are x

2(O)

=

b2 and xl(l)

=

bl, then I • Hence, eN

n a = sup IIQ(t) II •

t

o

If some restarts are necessary then we don't know the corresponding fundamental matrix anymore.

For general boundary conditions we can not prospect how large the

conditionnumber of BOQ(O) + B1Q(I) will be. However, if the estimate

is realistic, then Q(t) is a much better choice than M(t,O), since

H(Q(O»

,....,IIN~!(O)II,

which is, in general, large.

3.4.

3.4.1. The solution of (3.1) w\"th .{hao C.ClV\ bE'. -f0uV\C{ b~ <::.oiv'",~

(3.21 )

I f some fundamental matrix X(t) is known, then (3.21) is solved by

(3.22a)

(3.22b) x(t) == X(t)y •

Write B

=

BOX(O) + BIX(I). Suppose b, Band X(t) are known exactly.

Let y and x(t) be the computed solutions and n the machine precision.

0.23)

So

(3.24)

in By :: b : IIL\yll _ lIy - yll ~ nH(B)lIyll

in x(t)

=

X(t)y : (with ~(t)

=

X(t)y)

II L\x(t) II :: II x(t) - x(t) II

~ II x(t) - ~(t) II + II ~(t) - x(t) II

~ nil X(t) IIllyll + IIX(t) IIIIL\yll

~ nil X(t) II lIyll (l + H(B»

= nil X(t) IIIIB-1bll (1 + H(B» •

sup II6.x(t) II

~

n IIB-1bll (I + H(B» sup IIX(t) "

t t

S n lib I! (I + H(B)) sup IIX(t) II liB-III

t

~

n H(B)lIbll (l + H(B» sup IIS-I(t)II •

t

Hence, if X(t) is such that H(B) is not large, then the conditionnumber

of the solutionmethod (3.22) is equivalent to the conditionnumber of the

BVP, and is therefore called stable.

3.4.2. Comparing the conditionnumber of the original system (3.1) and the

trans-formed system (3.25)

1

d

°

dt x (t) BOUOXO(O)

° °

=

A (t)x (t)

o

-1 0 -1

where x (t)

=

U

_o

x(t) and A (t)

=

U

_o

A(t)U

O (UO orthogonal), we see that

in the II • 11₂-norm they are exactly equal. Hence, the conditionnumber of

a linear BVP does not change significantly by an orthogonal transformation.

In 3.3.2 and 3.4.1 we thus have seen that if (3.1) has separated boundary conditions then a single Riccati-transformation, if possible, gives a stable solution method.

Under the same conditions a modified Riccati-transformation will, in general,

also give a stable solution method. This is explained by the fact that X(O)

is like [small k

order lJ and X(I)_ like [order 1

n-k k

So BOX(O) + BIX(I) will be of the form

[ - small order 1 order 1] , small smallJ

-

n-k

which has a small conditionnumber. However, for solving (3.22bl we still have to invert an nXn-matrix.

3.4.3. For the repeated Riccati-transformations another technique is used. The equation (3.21) is not solved by (3.22), but more directly. Using (1.7)

and t:~. -Hj we have

.

h

tt,

t-;;) ':

Lt~

M

l Lt) Ll l-1 . - 1 _ " 1

Le t t € [t., t. 1 ] ~ ~+ Then M(t,O)-(3.26) (i € {O,-,m- I}) • M(t,t.)M(t.,t. I) ••• M(tl,t O) ~ ~~-i - i -) -1 i-I U.Q (t){Q) (t.)U.

u.

lQ (t.) ~ ~ ~ ~- ~

In the same way we can find an expression for M(l,t). If sufficient restarts are done then the matrices M(O,t) and M(l,t) are constructed in a stable way.

Now (3.21) can be solved. The conditionnumber of the solutionmethod is equal

to H(BOM(O,t) + BIM(l,t». BO and BI can be scaled such that this is

equi-valent to II [BOM(O, t) + Bl MO, t) III •

However, solving (3.21) by construction of the matrices M(O,t) and M(l,t) like in (3.26) is not recommendable.

Actually we then solve the system (1.49) by simple elimination, which is not the cutest way. We do see that the conditionnumber of the solution method (1.49) is equivalent to the conditionnumber of the problem, since this is even true for elimination. Hence, this method is stable.

The same can be said about the repeated modified Riccati-transformations,

although inver 50 ioV'l of the fundamental matrices

(3.27)