Next: Stepsize control Up: Correction Previous: Pseudo-arclength continuation   Contents

Moore-Penrose continuation

CL_MatCont implements a continuation method that is slightly different from the pseudo-arclength continuation.

Definition 1 Let A be an N×(N+1) matrix with maximal rank. Then the Moore-Penrose inverse of A is defined by A+ = AT(AAT)-1.

Let A be an N×(N+1) matrix with maximal rank. Consider the following linear system with x,v Î RN+1,b Î RN:

Ax
=
b
(8)
vTx
=
0
(9)

where x is a point on the curve and v its tangent vector with respect to A, i.e. Av=0. Since AA+b = b and vTA+b = áAv,(AAT)-1bñ = 0, a solution of this system is

x = A+b.
(10)

Suppose we have a predicted point X0 using (1). We want to find the point x on the curve which is nearest to X0, i.e. we are trying to solve the optimization problem:


min
x 
{ ||x-X0|| | F(x)=0 }
(11)

So, the system we need to solve is:

F(x)
=
0
(12)
wT(x-X0)
=
0
(13)

where w is the tangent vector at point x. In Newton's method this system is solved using a linearization about X0. Taylor expansion about X0 gives:

F(x)
=
F(X0) + Fx(X0)(x-X0) + O(||x-X0||2)
(14)
wT(x-X0)
=
vT(x-X0) + O(||x-X0||2) .
(15)

So when we discard the higher order terms we can see using (8) and (10) that the solution of this system is:

x = X0 - Fx+(X0)F(X0) .
(16)

However, the null vector of Fx(X0) is not known, therefore we approximate it by V0 = vi, the tangent vector at xi. Geometrically this means we are solving F(x)=0 in a hyperplane perpendicular to the previous tangent vector. In other words, the extra function g(x) in (2) becomes:

gk(x) = áx-Xk,Vk ñ,
(17)

where Fx(Xk-1)Vk=0 for k=1,2,¼. Thus, the Newton iteration we are doing is:

Xk+1
=
Xk - H-1x(Xk,Vk)H(Xk,Vk)
(18)
Vk+1
=
Vk - H-1x(Xk,Vk)R(Xk,Vk)
(19)


H(X,V) = æ
ç
ç
è
F(X)
0
ö
÷
÷
ø
,  Hx(X,V) = æ
ç
ç
è
Fx(X)
VT
ö
÷
÷
ø
(20)


R(X,V) = æ
ç
ç
è
Fx(X)V
0
ö
÷
÷
ø
 .
(21)

One can prove that under the same conditions as for the pseudo-arclength continuation, the Newton iterations (18) and (19) converge to a point on the curve xi+1 and the corresponding tangent vector vi+1, respectively. In the pseudo-arclength continuation, we had to compute a tangent vector when a new point was found. In this case however, we already compute the tangent vectors Vk at each iterate (19), so we only need to normalize the computed tangent vectors.