The Transpose-Free QMR (TFQMR) [7] for solving general non-Hermitian linear systems is derived from CGS algorithm. By incorporating the QMR approach, residual norm of TFQMR produces smoother convergence behavior than CGS. In the first, we derive TFQMR from the CGS algorithm. We double all subscripts in the CGS algorithm, that is to say
α2j = (r2j, r0∗)/(Ap2j, r0∗) (4.2)
q2j = u2j− α2jAp2j (4.3)
x2j+2 = x2j+ α2j(u2j+ q2j) (4.4) r2j+2 = r2j − α2jA(u2j+ q2j) (4.5) β2j = (r2j+2, r0∗)/(r2j, r0∗) (4.6)
u2j+2 = r2j+2+ β2jq2j (4.7)
p2j+2 = u2j+2+ β2j(q2j+ β2jp2j) (4.8) Observe that approximate solution x2j+2 in (4.4) can be split into the following two half-steps:
x2j+1 = x2j + α2ju2j (4.9)
x2j+2 = x2j+1+ α2jq2j. (4.10)
When m is odd, set ½
um = qm−1
αm = αm−1 . (4.11)
Then (4.9) and (4.10) can be simplified by (4.11). Whether m is even or odd, the single equation is
xm = xm−1+ αm−1um−1.
It must be noted that the intermediate iterates xm does not exist in the original CGS algorithm with m is odd. And we define the n × m matrix,
Um = [u0, . . . , um−1]
and the vector
zm = (α0, . . . , αm−1)T. General iterate xm and residual vector rm are
xm = xm−1+ αm−1um−1 (4.12)
= x0+ Umzm, (4.13)
rm = rm−1− αm−1Aum−1 (4.14)
= r0− AUmzm. (4.15)
From a result of (4.14), we have
Aum−1 = 1
αm−1(rm−1− rm).
Judging from the above relation, we can translate into matrix form, this relation becomes
AUm = Rm+1B¯m (4.16)
where Rm+1 is the n × (m + 1) matrix,
Rm+1 = [r0, r1, . . . , rm]
and where ¯Bm is the (m + 1) × m matrix with the following form,
B¯m =
1 0 0 · · · 0
−1 1 0 · · · 0 0 −1 1 · · · 0 ... ... ... ... ...
0 · · · 0 −1 1 0 0 · · · 0 −1
× diag
½ 1 α0, 1
α1, · · · , 1 αm−1
¾ .
Let the (m + 1) × (m + 1) scaling matrix is
4m+1 = diag{δ0, δ1, . . . , δm},
in addition, each inverse of diagonal element can rescale the corresponding column of Rm+1 equal to one. Then, the relation of (4.16) becomes
AUm = Rm+14−1m+14m+1B¯m
= ˜Vm+1H¯m, (4.17)
if we define the n × (m + 1) matrix ˜Vm+1 = Rm+14−1m+1 and (m + 1) × m matrix ¯Hm = 4m+1B¯m. With this result, equation (4.15) becomes
rm = r0− AUmzm = Rm+1£
e1− ¯Bmzm¤
= Rm+14−1m+1£
δ0e1− 4m+1B¯mzm¤
. (4.18)
This will lead us further into a consideration of whether we can exploit above relations by QMR. But before we come on to that, let us pause here to look at the following lemma.
Lemma 1. Let ˜Rm be the m × m upper part of the matrix Qm−1H¯m in FOM and let Rm be the m × m upper part of the matrix QmH¯m in GMRES. Similarly, let ˜gm be the vector of the first m components of Qm−1(βe1) and let gm be the vector of the first m components of Qm(βe1). Define
˜
ym = ˜R−1m ˜gm, ym = R−1m gm
the y vectors obtained for an m-dimensional FOM and GMRES methods, respectively.
Then
ym−
µym−1 0
¶
= c2m µ
˜ ym−
µym−1 0
¶¶
in which cm is the cosine used in the m-th rotation Ωm, as defined in DQGMRES.
We need mention here only the result of lemma, this proof of lemma can be seen in [9, p.180]. Let us return to our main subject to derive TFQMR. By above lemma, we know it is also valid for the CGS/TFQMR pair. This relation provides a starting point to derive the TFQMR algorithm. Thus, the TFQMR iterates satisfy the following relation
xm− xm−1 = c2m(xCGSm − xm−1). (4.19) Setting the vector
dm = 1
αm−1(xCGSm − xm−1) (4.20)
and the scalar
ηm = c2mαm−1. (4.21)
Then (4.19) can be expressed as
xm = xm−1+ ηmdm. (4.22)
Moreover, by (4.12), (4.19), (4.20), and (4.21), we have dm = 1
αm−1
£(xCGSm − xCGSm−1) + (xCGSm−1 − xm−2) − (xm−1 − xm−2)¤
= um−1+ 1
αm−1(xCGSm−1 − xm−2) − 1
αm−1(xm−1− xm−2)
= um−1+1 − c2m−1
αm−1 (xCGSm−1 − xm−2)
= um−1+1 − c2m−1 c2m−1
ηm−1 αm−1dm−1.
Define the m-st tangent by θm = sm/cm, then we have the new relation as below dm+1 = um+θm2ηm
αm dm. (4.23)
By the structure of ¯Hm in (4.17), the angle used in the (j + 1)-th rotation can obtain by (4.1), that is
sj+1= q−δj+1
τj2 + δj+12 , cj+1= q τj
τj2+ δ2j+1, θj+1 = −δj+1
τj , (4.24)
where τj be (j + 1)-th diagonal element of Ωj· · · Ω1H¯m. Moreover, after j + 1 rotations, next diagonal element δj+1 becomes τj+1, which is
τj+1 = δj+1× cj+1 = qτjδj+1
τj2+ δj+12
= −τj × sj+1 = −τjθj+1cj+1. (4.25) Since only the square of scalar θj+1 is invoked in the update of the direction dm+1. We can ignore the negative symbol in (4.24) and (4.25). So far, we have seen the following relations:
• dm+1 = um+ (θ2m/αm)ηmdm
• θm+1 = δm+1/τm
• cm+1 = (1 + θ2m+1)−12
• τm+1 = τmθm+1cm+1
• ηm+1 = c2m+1αm.
The question which we must consider next is to derive the remain terms. We note a little earlier that the vectors rm in (4.14) are no longer the actual residuals, we rewrite the new notation by wm. Then we have
wm = wm−1− αm−1Aum−1.
We may note, in passing, that the scalar δm equal to kwmk2. Setting the new vector v2j = Ap2j and multiplying (4.8) by matrix A,
v2j = Au2j+ β2j−2(Aq2j−2 + β2j−2Ap2j−2)
= Au2j+ β2j−2(Au2j−1+ β2j−2v2j−2). (by (4.16)) The same observation applies to (4.3) and (4.7), we have
u2j+1 = u2j − α2jv2j u2j+2 = w2j+2+ β2ju2j+1.
To sum up all relations above, we have TFQMR algorithm as follow:
ALGORITHM 17. TFQMR
1. w0 = u0 = r0 = b − Ax0, v0 = Au0, d0 = 0 2. τ0 = kr0k2, θ0 = η0 = 0.
3. Choose r0∗ such that ρ0 = (r∗0, r0) 6= 0.
4. For m = 0, 1,. . . , until convergence Do:
5. If m is even then
6. αm+1 = αm = ρm/(vm, r∗0) 7. um+1 = um− αmvm
8. EndIf
9. wm+1 = wm− αmAum 10. dm+1 = um+ (θ2m/αm)ηmdm
11. θm+1 = kwm+1k2/τm; cm+1 = (1 + θ2m+1)−12 12. τm+1 = τmθm+1cm+1; ηm+1 = c2m+1αm 13. xm+1 = xm+ ηm+1dm+1
14. If m is odd then
15. ρm+1 = (wm+1, r0∗); βm−1 = ρm+1/ρm−1 16. um+1 = wm+1+ βm−1um
17. vm+1 = Aum+1+ βm−1(Aum+ βm−1vm−1)
18. EndIf
19. EndDo
5 Algorithm and Numerical Results
5.1 General Solution Algorithm
In the first, we refer to the finest grid (with 41 × 129 nodes) as level 1. Two additional coarser grids are constructed by successively discarding every other mesh line from one grid to the coarser one. These grids are numbered with increasing level number, we have
• level 1: 41 × 129 grid,
• level 2: 21 × 65 grid,
• level 3: 11 × 33 grid.
According to the characteristic of the results are given in section 1.2, we have a nonlinear system of the discretized equations in residual form,
F (Φ) = 0. (5.1)
Assume that nonlinear system start from an initial guess Φ0 and have an accurate solution Φ∗. The vector Φ has NcNxNy three components for two dimensional problems, where Nc is the number of the solution components; Nx and Ny are the number of grid points in the r and z direction, respectively. In order to stabilize convergence, we use a damped Newton method [3] instead of Newton method. Suppose an initial solution Φ0 is close to the solution Φ∗ enough, the equations (5.1) can be solved by using a damped Newton method, that is
J(Φn)(Φn+1− Φn) = −λnF (Φn), n = 0, 1, . . . (5.2) where J(Φn) is the Jacobian matrix at Φn with the form
J(Φn) = ∂F
∂Φ(Φn),
and the damping parameter λn is in the range of (0 < λn ≤ 1). We denote the update vector as 4Φn = Φn+1− Φn and (5.2) with convergence tolerance
k 4 Φnk2 < 10−5.
In our numerical applications, we use likely modified Newton method in (5.2). This means that the Jacobian matrix is re-evaluated only twice Newton iterative step. If the rate of convergence is too slow, we form a new damping parameter.
We know that the governing equations of a combustion process are difficult to solve and a good initial solution estimate for an iteration process is hard to determine. For eliminating these difficulties, a time relaxation process is used. Due to the nonlinearity of the governing equations of a combustion process, we can use a pseudo transient process to produce a parabolic in time problem and bring the starting estimate into the convergence domain of the steady Newton method. Thus, the unsteady form of the governing equations can be obtained by adding the unsteady term to the steady-state equations. We have
F(Φn+1) = F (Φn+1) + DΦn+1− Φn
4tn+1 = 0, (5.3)
where D is a diagonal scaling matrix with nonnegative entries and 4tn+1 = tn+1−tnis n-st time step. In our calculations of program, in time dependent part, we use modified Newton method to solve the nonlinear system (5.3) which is similar to the system of equations in (5.2). In steady-state part, we use one way multigrid method and damped Newton method to solve the nonlinear system (5.1). One way multigrid method means that the coarser meshes are used only to initialize the next finer ones. Our goal is to obtain a converged numerical solution on level 1. For that purpose, we solve the nonlinear problem (5.1) and (5.3) starting at the coarsest level and ending at the finest. To summarize above informations, we write down these phases as following sample processes:
1. Time stepping on level 3 (in coarsest level).
2. Steady Newton iterations on level 3 and interpolation of the numerical solution from level 3 onto level 2.
3. Steady Newton iterations on level 2 and interpolation of the numerical solution from level 2 onto level 1.
4. Steady Newton iterations on level 1.