Volume 2010, Article ID 539087,23pages doi:10.1155/2010/539087
Research Article
Stability of Nonlinear Autonomous Quadratic Discrete Systems in the Critical Case
Josef Dibl´ık,
1, 2Denys Ya. Khusainov,
3Irina V. Grytsay,
3and Zden ˘ek ˇSmarda
11Department of Mathematics, Faculty of Electrical Engineering and Communication, Brno University of Technology, Technick´a 8, 616 00 Brno, Czech Republic
2Department of Mathematics and Descriptive Geometry, Faculty of Civil Engineering, Brno University of Technology, Veveˇr´ı331/95, 60200 Brno, Czech Republic
3Department of Complex System Modeling, Faculty of Cybernetics, Taras,
Shevchenko National University of Kyiv, Vladimirskaya Str., 64, 01033 Kyiv, Ukraine
Correspondence should be addressed to Josef Dibl´ık,[email protected] Received 28 January 2010; Accepted 11 May 2010
Academic Editor: Elena Braverman
Copyrightq2010 Josef Dibl´ık et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Many processes are mathematically simulated by systems of discrete equations with quadratic right-hand sides. Their stability is thought of as a very important characterization of the process. In this paper, the method of Lyapunov functions is used to derive classes of stable quadratic discrete autonomous systems in a critical case in the presence of a simple eigenvalueλ1 of the matrix of linear terms. In addition to the stability investigation, we also estimate stability domains.
1. Introduction
The main results on the stability theory of difference equations are presented, for example, by Agarwal 1, Agarwal et al. 2, Chetaev 3, Elaydi 4, Halanay and R˘asvan 5, Lakshmikantham and Trigiante6, and Martynjuk7. Instability problems are considered, for example, in8–10by Slyusarchuk. Note that stability and instability results often have a local character and are usually obtained without any estimation of the stability domain, or without investigating the character of instability. Moreover, it should be emphasized that global instability questions have only been discussed for linear systems.
Many processes and phenomena are described by differential or difference systems with quadratic nonlinearities. Among others, let us mention epidemic and populations models, models of chemical reactions, and models for describing convection currents in the atmosphere.
The stability of a zero solution of difference systems
xk1 fxk, 1.1
wherek 0,1, . . ., and x x1, . . . , xnT with differentiablef f1, . . . , fnT : Rn → Rn, is very often investigated by linearly approximating system1.1in question by using the matrix of linear terms
xk1 Axk gxk, 1.2
whereAf0,0, . . . ,0is the Jacobian matrix offat0,0, . . . ,0, andgx fx−Ax. This approach becomes unsuitable in what is called a critical case, that is, when the spectral radius of the matrixρA 1 because, among all systems1.2, there are classes of stable systems as well as classes of unstable systems. Concerning this, we formulate the following known resultssee, e.g., Corollary 4.344, page 222and Theorem 4.384, page 226.
Theorem 1.1. 1IfρA<1, then the zero solution of 1.2is exponentially stable.
2IfρA 1, then the zero solution of 1.2may be stable or unstable.
3IfρA>1 andgxisoxasx → 0, then the zero solution of 1.2is unstable.
In this paper, we consider a particular critical case when there exists a simple eigenvalue λ 1 of the matrix of linear terms and the remaining eigenvalues lie inside a unit circle centered at origin. The purpose of this paper is to obtain using the method of Lyapunov functionsconditions for the stability of a zero solution of difference systems with quadratic nonlinearities in the above case and derive classes of stable systems. In addition to the stability investigation, we estimate the stability domains as well. The domains of stability obtained are also called guaranteed domains of stability. Preliminary results in this direction were published in11.
1.1. Quadratic System and Preliminary Consideration In the sequel, the norms used for vectors and matrices are defined as
x
n
i1
x2i 1/2
1.3
for a vectorx x1, . . . , xnTand
A
λmaxATA1/2
1.4
for anym×nmatrixA. Here and in the sequel,λmax· orλmin·is the maximalor minimal eigenvalue of the corresponding symmetric and positive-semi- definite matrixsee, e.g., 12.
Consider a nonlinear autonomous discrete system with a quadratic right-hand side
xik1 n
s1
aisxsk n
s,q1
bisqxskxqk, i1. . . , n, 1.5
wherek0,1, . . .and the coefficientsaisandbisqwe assume thatbisq biqsare constant. As emphasized, for example, in3,7,12, system1.5can be written in a general vector-matrix form
xk1 Axk XTkBxk, k0,1, . . . , 1.6 where
aA{ais}, i, s1,2. . . , n, is ann×nconstant square matrix,
bmatrixXT{X1T, XT2, . . . , XTn}isn×n2rectangular and all the elements of then×n matricesXiT, i 1, . . . , n, are equal to zero except the ith row with entriesxT x1, x2, . . . , xn, that is,
XTik
⎛
⎜⎜
⎜⎜
⎜⎝
0 0 · · · 0
· · · · x1 x2 · · · xn
· · · · 0 0 · · · 0
⎞
⎟⎟
⎟⎟
⎟⎠, 1.7
cmatrix BT {B1, B2, . . . , Bn}isn2×nrectangular and then×nconstant matrices Bi{bisq}, i, s, q1, . . . , n, are symmetric.
The stability of the zero solution of1.6depends on the stability of the matrixA. IfρA<1, then the zero solution of1.6is exponentially stable for an arbitrary matrixBbyTheorem 1.1.
In this case, matrixBonly impacts on the shape of the stability domain of the equilibrium state. If the zero solution of1.6is investigated on stability by the second Lyapunov method and an appropriate Lyapunov function is taken as the quadratic formVx xTHxwith a suitablen×nconstant real symmetric positive-definite matrixH, which is defined below, then the first differenceΔV along the trajectories of1.6equals
ΔVxk Vxk1−Vxk xTk1Hxk1−xTkHxk
Axk XTkBxkT H
Axk XTkBxk
−xTkHxk
xTkATxTkBTXk H
Axk XTkBxk
−xTkHxk xTk
ATHA−H
ATHXTkBBTXkHABTXkHXTkB xk xTk
ATHA−H
2BTXkHABTXkHXTkB xk
1.8
sinceATHXTkBT BTXkHA.
Since ρA < 1, for arbitrary positive-definite symmetric matrix C, the matrix Lyapunov equation
ATHA−H−C 1.9
has a unique solutionH—a positive-definite symmetric matrixe.g.,4, Theorem 4.30, page 216. We use such matrixHto estimate the stability domain. Then, as follows from1.8,
ΔVxk≤ −
λminC−2B · HA · xk − B2· H · xk2
· xk2. 1.10
Analysing1.10, we deduce that the first differenceΔVxkwill be negative definite if B2· H · xk22B · HA · xk ≤λminC, 1.11 that is, it will be negative definite in a neighborhoodUδ {x∈Rn:x< δ}of the steady- statexk ≡ 0, k 0,1, . . . ,ifδ is sufficiently small. In the case considered, the domain of stability can be described by means of two inequalities. The first inequality 1.11defines a part of the spaceRn, where the first differenceΔVxkis negative definite. The second inequality
Vx≤r, x∈Rn, r >0, 1.12 describes points inside a level surface. The guaranteed domain of stability is given by inequality1.12ifr is taken so small that the domain described by1.12is embedded in the domain described by inequality1.11.
Considering the investigated critical case, we will deal with a different structure of the right-hand side of the inequality from that seen in 1.10. Namely, we will show that, unlike the right-hand side of the inequality forΔVxkthat is multiplied byxk2with dimxk nin1.10, in the critical case considered, the right-hand side of the inequality or equalityforΔVxkwill be multiplied only by a termxn−1k2 with dimxn−1k n−1< nsee2.21in the casen2 and2.69in the general case below.
2. Main Results
In this section we derive the classes of the stable systems1.6in a critical case when the matrixAhas one simple eigenvalueλ1.
2.1. Instability in One-Dimensional Case
We start by discussing a simple scalar equation with the eigenvalue of matrixAequaling one, that is,a11 1. Then1.6takes the form
xk1 xk bx2k, k0, . . . , 2.1
and it is easy to see that the trivial solution is unstable for an arbitraryb /0to show this, we can apply, e.g., Theorem 1.154, page 29.
This elementary example shows that stability in the case of system 1.6 has an extraordinary significance and the results on stabilityforn /1lose their meaning forn1 when we deal with instability. We show that, ifn /1 andBsatisfies certain assumptions, the zero solution is stable. Moreover, the shape of the guaranteed domain of stability will be given.
We divide our forthcoming analysis into two parts. In the first one we give an explicit coefficient criterion in the subcase ofn2. Then we consider the generaln-dimensional case.
2.2. Stability in the General Two-Dimensional Case
Letn 2. Then system1.6with the matrixAhaving a simple eigenvalueλ 1 reduces after linearly transforming the dependent variables if necessaryto
x1k1 ax1k
b111x21k 2b112x1kx2k b122x22k , x2k1 x2k
b211x21k 2b212x1kx2k b222 x22k .
2.2
We will assume that|a|<1. Define auxiliary numbers as follows:
αh 1−a2
, β1hab111, β22hab112b112 , γ1 h
b1112
b2112
, γ2 4h b121 2
, δ12hb111b121 ,
2.3
wherehis a positive number.
Theorem 2.1. Lethandr be positive numbers. Assume that|a|<1 andb212 b122 b222 0. Then the zero solution of system2.2is stable in the Lyapunov sense and a guaranteed domain of stability is given by the inequality
hx21x22≤r2 2.4
ifris taken so small that the domain described by2.4is embedded in the domain
γ1x212δ1x1x2γ2x222β1x12β2x2 ≤α. 2.5
If, moreover,b211·b112/0, then a guaranteed domain of stability can be described using inequality
hx21x22≤r∗2 2.6
with
r∗ min
x1,x2
hx21x22, 2.7
wherex1, x2runs over all real solutions of the nonlinear system with unknownsx1andx2:
γ1x212δ1x1x2γ2x222β1x12β2x2α, hx1
δ1x1γ2x2β2
−x2
γ1x1δ1x2β1 0.
2.8
Proof. Define
B1
b111 b112 b121 b122
, B2
b211 b212 b212 b222
, x2y, x x1
y
. 2.9
We rewrite system2.2as
x1k1 ax1k xTkB1xk,
yk1 yk xTkB2xk. 2.10
To investigate the stability of the zero solution, we use, in accordance with the direct Lyapunov method, an appropriate Lyapunov functionV. Let a matrixH, defined as
H
h h12
h12 h22
, 2.11
where instead of the entryh11we put the numberh, be positive definite. We set
Vxk V
x1k, yk
:xTkHxk
x1k, yk h h12 h12 h22
x1k yk
hx21k 2h12x1kyk h22y2k.
2.12
The first difference of the functionV along the trajectories of system2.10equals
ΔVxk Vxk1−Vxk
hx21k1 2h12x1k1yk1 h22y2k1
−hx21k−2h12x1kyk−h22y2k h
ax1k xTkB1xk2 2h12
ax1k xTkB1xk
yk xTkB2xk h22
yk xTkB2xk2
−hx21k−2h12x1kyk−h22y2k h
a2−1
x21k 2h12a−1x1kyk h
2ax1k
xTkB1xk
xTkB1xk2 2h12
ax1k
xTkB2xk
yk
xTB1xk
xTkB1xk
xTkB2xk h22
2yk
xTkB2xk
xTkB2xk2 .
2.13
It is easy to see thatΔV does not preserve the sign ifh12/0. Therefore, we puth12 0 and ΔVreduces to
ΔVxk h a2−1
x21k h
2ax1k
xTkB1xk
xTkB1xk2 h22
2yk
xTkB2xk
xTkB2xk2 .
2.14
In the polynomial ΔV, with respect to x1 and y, we will put together the third-degree terms the expression F3x1k, yk below and the fourth-degree terms the expression F4x1k, ykbelow. In the computations we use the formulas
xTkBixk bi11x12k 2bi12x1kyk bi22y2k, i1,2, xTkBixk2
bi112
x41k 4 bi122
x21ky2k bi222
y4k
4bi11bi12x13kyk 2bi11bi22x21ky2k 4bi12bi22x1ky3k, i1,2.
2.15
We get
ΔVxk h a2−1
x21k F3
x1k, yk F4
x1k, yk
, 2.16
where
F3
x1k, yk
2hab111x13k 2
2hab112h22b112
x21kyk 2
hab1222h22b212
x1ky2k 2h22b222 y3k, F4
x1k, yk
h
b1112 h22
b2112
x41k 4
hb111b112h22b211b212
x31kyk 2
2h
b1122
2h22 b2122
hb111 b122h22b211b222
x21ky2k 4
hb112b221 h22b212b222
x1ky3k
h b1222
h22 b2222
y4k.
2.17
Analysing the increment of V, we see that, if |a| < 1, ΔV will be nonpositive in a small neighborhood of the zero solution if the multipliers of the termsx1y2,y3andx1y3are equal to zero and the multiplier of the termy4is nonpositive, that is, if
hab1222h22b122 0, h22b2220, hb121 b122h22b122 b222 0, h
b221 2 h22
b2222
≤0.
2.18
As long as the Lyapunov function is positive definite,h >0 andh22 >0. Therefore, conditions 2.18hold if and only if
b1220, b2120, b222 0. 2.19
Then, system2.2turns into
x1k1 ax1k
b111x21k 2b121 x1kx2k , x2k1 x2k b211x21k
2.20
andΔVwithout loss of generality, we puth22 1, i.e.,Vx1, y hx21y2into ΔVxk −
h 1−a2
−2hab111x1k−2
2hab112b112
yk− h
b1112
b2112
x12k
−4hb111b112x1kyk−4h b121 2
y2k
x21k −
α−2β1x1k−2β2yk−γ1x21k−2δ1x1kyk−γ2y2k x21k.
2.21 The first difference of the Lyapunov function is nonpositive in a sufficiently small neighborhood of the originthis is becauseh > 0, |a|< 1, andα h1−a2 > 0. In other words, the zero solution is stable in the Lyapunov sense.
Now we will discuss the shape of the guaranteed domain of stability. It can be defined by the inequalities
γ1x212δ1x1yγ2y22β1x12β2y≤α, hx21y2≤r2,
2.22
where r > 0. This means that inequalities 2.4 and 2.5 are correct. Both inequalities geometrically express closed ellipses ifb211·b112/0. For the second inequality, this is obvious.
For the first one, this follows from the following inequalities:γ1>0, γ2>0 and
γ1γ2−δ21
h b1112
b2112
·
4h b1122
−4
hb111b121 2 4h
b211b121 2
>0. 2.23
Moreover, forr → 0, the ellipse2.4
hx21y2≤r2 2.24
is contained because it shrinks to the originin the ellipse 2.5, that is, there exists such rr∗that, forr ∈0, r∗, the ellipse2.4lies inside the ellipse2.5without any intersection points and, forr r∗, there exists at least one common boundary point of both ellipses. Let us find the valuer∗. It is characterized by the requirement that the slope coefficientsk1and k2of both ellipses are the same at the point of contact. Therefore
k1−γ1x1δ1yβ1
δ1x1γ2yβ2, k2−hx1
y , 2.25
where we assumewithout loss of generalitythat the denominators are nonzero. Thus, we get a quadratic system of two equations to find the contact pointsx1, y:
γ1x212δ1x1yγ2y22β1x12β2yα, hx1
δ1x1γ2yβ2
−y
γ1x1δ1yβ1
0. 2.26
For the corresponding values ofr, we have
hx21y2r2. 2.27
In accordance with the geometrical meaning of the above quadratic system, we take such a solutionx1, yas a defintion of the minimal positive value ofrand setr∗r.
Example 2.2. Consider a system
xk1 0.5xk x2k−4xkyk,
yk1 yk x2k. 2.28
In our case,n 2,a0.5 <1, andb122 b212 b222 0. Therefore, byTheorem 2.1, the zero solution of system2.28is stable in the Lyapunov sense. We will find the guaranteed domain of stability. We have
b1111, b112 −2, b2111, 2.29
andb211·b112−2/0. Seth2. Then
αh 1−a2
21−0.25 1.5, β1hab111 2·0.5·11,
β22hab121 b2112·2·0.5·−2 1−3, γ1h
b111 2
b112 2
2·12123, γ24h
b1122
4·2·−2232, δ12hb111b1122·2·1·−2 −8.
2.30
That is, the guaranteed domain of stability is given by the inequalities
3x2−16xy32y22x−6y≤1.5, 2.31
2x2y2≤r2 2.32
if r is so small that the domain described by inequality2.32is embedded in the domain described by inequality2.31. We consider the case when the ellipse2.32is embedded in
−1
−0.5 0 0.5 1
−1.5 −1 −0.5 0 0.5 1 1.5 Figure 1: Graphical solution of system2.33and2.34.
the ellipse2.31and the boundaries of both ellipses have only one intersection point. Solving the system2.8, that is, the system
3x2−16xy32y22x−6y1.5, 2.33
2x
−8x32y−3
−y
3x−8y1
0, 2.34
with Mathematica software, we get the solutionsseeFigure 1where thex-axis is identified with the horizontal line and the y-axis is identified with the vertical line, the blue ellipse graphically depicts equation 2.33, and the red hyperbola graphically depicts equation 2.34:
x, y
x1, y1
−1.60766,−0.31220, x, y
x2, y2
−0.03568,−0.32187, x, y
x3, y3
0.01952,−0.13664, x, y
x4, y4
1.10728,0.37750.
2.35
Then, in accordance with2.7,
r∗ min
i1,2,3,4
2xi2yi2
2x32y32 . 0.1369, 2.36
−0.4
−0.2 0 0.2 0.4
−1.5 −1 −0.5 0 0.5 1 Figure 2: The guaranteed domain of stability.
and the guaranteed domain of stability
2x2y2≤r∗2 0.13692 2.37
obtained from2.31,2.32is depicted inFigure 2as an ellipsoidal domain shaded in red and bounded by the thick red ellipse, with the identification ofx-axis andy-axis being the same as before. Here, the domain2.31is bounded by the blue ellipse2.33.
2.3. Stability in the Generaln-Dimensional Case
Consider system1.6inRn. Assume that the matrixAhas a simple eigenvalue that is equal to unity with the others lying inside the unit circle. After linearly transforming the dependent variables if necessary, we can assume, without loss of generality, that the matrix Aof the linear terms in a block form, that is,
A A0 θ
θT 1
, A0 aij
, i, j1,2, . . . , n−1, 2.38
whereθ 0,0, . . . ,0T, is then−1-dimensional zero vector and all the eigenvalues of the matrix A0 lie inside the unit circle. In order to formulate the next result and its proof, we
have to introduce some new definitionsthey copy the ones used inSection 1.1, but we use dimension or sizen−1 instead ofnand note this change as a subscript if necessary:
xn−1 x1, x2, . . . , xn−1T, yxn,
B0i
⎛
⎜⎜
⎜⎜
⎜⎝
b11i b12i · · · bi1,n−1 b21i b22i · · · bi2,n−1
· · · · bin−1,1 bin−1,2 · · · bin−1,n−1
⎞
⎟⎟
⎟⎟
⎟⎠, i1,2, . . . , n, B
⎛
⎜⎜
⎝
b11n · · · bn−1,n1
· · · · bn−11n · · · bn−1,nn−1
⎞
⎟⎟
⎠,
BT
⎛
⎜⎜
⎝
b111 · · · b11,n−1 b211 · · · b12,n−1 · · · b1n−1,1 · · · b1n−1,n−1
· · · · b11n−1 · · · bn−11,n−1 bn−121 · · · bn−12,n−1 · · · bn−1n−1,1 · · · bn−1n−1,n−1
⎞
⎟⎟
⎠.
2.39
MatricesB0i, i1,2, . . . , n, are symmetric sincebisq biqs,i, s, q1,2, . . . , nseeSection 1.1.
Moreover, we assume that there exists a symmetric positive definiten−1×n−1matrix Hsuch that the symmetric matrix
CH−AT0HA0 2.40
is positive definite. Leth >0 be a positive number and
αλminC, β1 1
2
AT0HBT , β2 1
2
BHA03AT0HBT2h B0nT
, γ1
BHBT
hB0n2,
γ24BHBT, δ12B· H ·B.
2.41
Theorem 2.3. Lethandrbe positive numbers. Assume that
bnn1 bnn2 · · ·bnnn0, b1nn b2nn · · ·bnn−1,n0. 2.42
Then the zero solution of system1.6is stable by Lyapunov and the guaranteed domain of stability is described by the inequalities
γ1x22δ1xyγ2y22β1x2β2y≤α, 2.43
xTn−1Hxn−1hy2≤r2 2.44
ifris so small that the domain described by inequality2.44is embedded into the domain described by inequality2.43.
Proof. We will perform auxiliary matrix computations. With this in mind, we have defined ann−12×n−1matrixXn−1as
XTn−1
X1n−1T , XT2n−1, . . . , Xn−1n−1T
, 2.45
where all the elements of then−1×n−1matricesXTin−1, i1,2, . . . , n−1 are equal to zero except the rowi, which equalsxTn−1, that is,
XTin−1
⎛
⎜⎜
⎜⎜
⎜⎝
0 0 · · · 0
· · · · x1 x2 · · · xn−1
· · · · 0 0 · · · 0
⎞
⎟⎟
⎟⎟
⎟⎠. 2.46
Moreover, we define
avectorsYi, i 1,2, . . . , n−1, as a rown−1-dimensional vector with coordinates equal to zero except theith element, which equalsxn, that is,
Yi 0,0, . . . ,0, xn,0, . . . ,0, 2.47
b n−1×n−1zero matrixΘ,
cvectorsbi bi1n, bi2n, . . . , bin−1,nT,i1,2, . . . , n, dvectorb b1nn, b2nn, . . . , bnnn−1T.
It is easy to see that
XTk
⎛
⎝X1n−1T k Y1Tk · · · XTn−1n−1k Yn−1T k Θ θ
θT 0 · · · θT 0 xTn−1k yk
⎞
⎠,
B
⎛
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎝ B01 b1 bT1 b1nn
· · · · Bn0 bn bTn bnnn
⎞
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎠ .
2.48
Now we are able to rewrite system1.6in an equivalent form xn−1k1
yk1
A0 θ
θT 1
xn−1k yk
X1n−1T k Y1Tk · · · XTn−1n−1k Yn−1T k Θ θ
θT 0 · · · θT 0 xTn−1k yk
×
⎛
⎜⎜
⎜⎜
⎜⎝ B10 b1 bT1 b1nn
· · · · Bn0 bn bTn bnnn
⎞
⎟⎟
⎟⎟
⎟⎠
xn−1k yk
A0r11 r12
r21 1r22
xn−1k yk
,
2.49
where
r11 r11
xn−1k, yk n−1
j1
Xjn−1T kB0jYjTkbjT
BTXn−1k Byk,
r12r12
xn−1k, yk n−1
j1
XTjn−1kbjYjTkbjnn
BTxn−1k byk,
r21r21
xn−1k, yk
xTn−1kB0nykbTn, r22 r22
xn−1k, yk
xTn−1kbnykbnnn .
2.50
Before the following computations, for the reader’s convenience, we recall that for then− 1×n−1matricesA,A1, 1×n−1vectors , 1,n−1×1 vectorsC,C1and 1×1 “matrices”
m,m1, we have A C
m
×
A1 C1 1 m1
A × A1C × 1 A× C1C ×m1
× A1m× 1 × C1m×m1
. 2.51
To investigate the stability of system1.6, we use the Lyapunov function
V
xn−1k, yk
xn−1T kHxn−1k hy2k
xTn−1k, ykH θ θT h
xn−1k yk
,
2.52
whereHHn−1is ann−1×n−1constant real symmetric and positive definite matrix.
Let us find the first difference of the Lyapunov function2.52along the solutions of2.49.
We get
ΔV
xn−1k, yk
xTn−1k1, yk1H θ θT h
xn−1k1 yk1
−
xn−1T k, ykH θ θT h
xn−1k yk
xTn−1k, ykA0r11 r12
r21 1r22 T
H θ θT h
A0r11 r12
r21 1r22
xn−1k yk
−
xn−1T k, ykH θ θT h
xn−1k yk
xTn−1k, ykAT0 r11T r21T r12T 1r22
H θ θT h
×
A0r11 r12 r21 1r22
− H θ θT h
×
xn−1k yk
2.53
or, using formula2.51, ΔV
xn−1k, yk
xn−1T k, ykAT0Hr11TH hr21T r12TH h1r22
×
A0r11 r12 r21 1r22
− H θ
θT h
×
xn−1k yk
xn−1T k, yk
×
⎛
⎝
AT0r11T
HA0r11 r21Thr21−H
AT0 r11T
Hr12r21Th1r22 r12THA0r11 1r22hr21 r12THr12 1r22h1r22−h
⎞
⎠
×
xn−1k yk
.
2.54 Using formulas2.50, we have
ΔV
xn−1k, yk
xn−1T k, yk
× c11 c12 c21 c22
xn−1k yk
, 2.55
where
c11 c11
xn−1k, yk
A0BTXn−1k Byk T
H
A0BTXn−1k Byk h
Bn0T
xn−1k bnyk, xTn−1kBn0ykbTn
−H,
c12 c12
xn−1k, yk
A0BTXn−1k Byk T
H
BTxn−1k byk h
1xTn−1kbnbnnnyk Bn0T
xn−1k bnyk , c21 c21
xn−1k, yk cT12,
c22 c22
xn−1k, yk
BTxn−1k byk T
H
BTxn−1k byk h
1xTn−1kbnbnnn yk2
−h.
2.56
Then ΔV
xn−1k, yk
xTn−1kc11xn−1k 2xTn−1kc12yk ykc22yk xTn−1
A0BTXn−1k Byk T
H
A0BTXn−1k Byk h
B0nT
xn−1k bnyk
xTn−1kBn0ykbTn
−H
xn−1k 2xn−1T k
A0BTXn−1k Byk T
H
BTxn−1k byk h
1xn−1T kbnbnnnyk Bn0T
xn−1k bnyk
yk yk
BTxn−1k bykT H
BTxn−1k byk h
1xn−1T kbnbnnnyk2
−h
yk.
2.57
After further computation, we get ΔV
xn−1k, yk xTn−1k
AT0HA0−HAT0HBTXn−1k Xn−1T BHA0 Xn−1T kBHBTXn−1k
Xn−1T kBHBBTHBTXn−1k yk
BTHA0AT0HBT
yk BTHBy 2k h Bn0T
xn−1kxTn−1kB0n 2h
Bn0T
xn−1kbTnyk hbny2kbTn
xn−1k 2xTn−1k
AT0HBTxn−1k AT0Hbyk Xn−1T kBHBTxn−1k XTn−1kBHbyk BH BTxn−1kyk BH by 2k h
Bn0T
xn−1k bnyk
hxn−1T kbn
Bn0T
xn−1k bnyk
hbnnnyk B0nT
xn−1k bnyk
yk
xTn−1kBH
BTxn−1k byk bTH
BTxn−1k byk
yk h
xTn−1kbn
2
h
bnnn yk2
2hxTn−1kbn2hbnnn yk 2hxTn−1kbnbnnnyk y2k.
2.58