Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 53
Текст из файла (страница 53)
Through mathematical manipulations the gradient isfound as11Q0 ð/Þ ¼ AT / þ A/ b22If A is symmetric A ¼ AT , then Eq. (10.94) impliesQ0 ð/Þ ¼ A/ bð10:94Þð10:95ÞThe minimum is obtained when Q0 ð/Þ ¼ 0, leading toQ0 ð/Þ ¼ 0 ) A/ ¼ bð10:96ÞTherefore minimizing Qð/Þ is equivalent to solving Eq. (10.1) and the solution ofthe minimization problem yields the solution of the system of linear equations.Now for the function Qð/Þ to have a global minimum it is necessary for thecoefficient matrix A to be positive definite, i.e., it should satisfy the inequality/T A/ [ 0 for all / 6¼ 0. This requirement can be established by considering therelationship between the exact solution / and its current estimate /ðnÞ . If e ¼ /ðnÞ / denotes the difference between the exact solution and the current estimate, thenEq. (10.93) gives1Qð/ þ eÞ ¼ ð/ þ eÞT Að/ þ eÞ bT ð/ þ eÞ þ c21111¼ /T A/ þ eT A/ þ /T Ae þ eT Ae bT / bT e þ c222 02111B1C¼ /T A/ bT / þ c þ @ eT A/ þ /T AeA bT e þ eT Ae|fflffl{zfflffl}|fflffl{zfflffl}2|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 22TTTQ ð /Þe b¼b eb e1¼ Qð/Þ þ eT Ae2ð10:97Þ10.3Iterative Methods335indicating that if A is positive definite, the second term will be always positiveexcept when e = 0, in which case the required solution would have been obtained.Moreover, when A is positive definite, all its eigenvalues are positive and thefunction Qð/Þ has a unique minimum.Thus with a symmetric and positive definite matrix a converging series of /ðnÞcan be derived such that/ðnþ1Þ ¼ /ðnÞ þ aðnÞ d/ðnÞð10:98Þwhere aðnÞ is some relaxation factor, and d/ðnÞ is related to the correction needed tominimize the said function at each iteration.
This can be accomplished in a varietyof ways leading to different methods.10.3.12 The Method of Steepest DescentThe method of steepest descent for solving linear systems of equations of the formgiven by Eq. (10.1) is based on minimizing the quadratic form given by Eq. (10.93).If / is a one dimensional vector with its components given by the scalar /, thenQð/Þ will represent a parabola. Finding the minimum of a parabolic functioniteratively starting at some point /0 , involves moving down along the parabola untilhitting the minimum.The same idea is used in N dimensions.
In this case Qð/Þ may be depicted as aparaboloid and the solution is iteratively found starting from an initial position /ð0Þand moving down the paraboloid until the minimum is reached. For quick convergence the sequence of steps /ð0Þ ; /ð1Þ ; /ð2Þ ; . . . should be selected such that thefastest rate of descent occurs, i.e., in the direction of Q0 ð/Þ. According toEq. (10.95) this direction is also given byQ0 ð/Þ ¼ b A/ð10:99ÞThe exact solution being /, the error and residual at any step n, denoted respectively by eðnÞ and rðnÞ , are computed asrðnÞeðnÞ ¼ /ðnÞ / ¼ b A/ðnÞ ¼ Q0 /ðnÞ)) rðnÞ ¼ AeðnÞð10:100ÞMoving linearly in the direction of the steepest descent, the value of / at step n + 1can be expressed in terms of the value of / at step n according to33610Solving the System of Algebraic Equations/ðnþ1Þ ¼ /ðnÞ þ aðnÞ rðnÞThe value of aðnÞ that minimizes Qð/Þ should satisfydðnþ1ÞQ/¼0daðnÞð10:101Þð10:102ÞThis can be expanded into"# T d/ðnþ1ÞTddðnþ1Þðnþ1ÞQ /Q /¼ 0 ) rðnþ1Þ rðnÞ ¼ 0¼0)ðnÞðnþ1ÞðnÞdadad/ð10:103Þindicating that the new step should be in a direction normal to the old step.
Thevalue of aðnÞ is calculated by using Eq. (10.103) as follows:TTrðnþ1Þ rðnÞ ¼ 0 ) b A/ðnþ1Þ rðnÞ ¼ 0hiTrðnÞ ¼ 0) b A /ðnÞ þ aðnÞ rðnÞTT) b A/ðnÞ rðnÞ ¼ aðnÞ ArðnÞ rðnÞð10:104Þ T T) rðnÞ rðnÞ ¼ aðnÞ rðnÞ ArðnÞ ðnÞ T ðnÞrrðnÞ)a ¼TðnÞðr Þ ArðnÞThe steepest descent algorithm can be summarized as follows:r(0) = b A(0)(choose residual as starting direction)iterate starting atr(n)(n)=b A=( n+1)until convergence(Compute the residual vector)( r ) r (Compute the factor in the orthogonal direction)( r( ) ) Ar( )T(n)(n)(n)=n(n)T+(n)n(n) (n)r(Obtain new)As presented above, the algorithm necessitates performing two matrix-vectormultiplications per iteration.
One of them can be eliminated by multiplying bothsides of Eq. (10.101) by −A and adding b to obtain10.3Iterative Methods337/ðnþ1Þ ¼ /ðnÞ þ aðnÞ rðnÞ ) b A/ðnþ1Þ¼ b A /ðnÞ þ aðnÞ rðnÞ ) rðnþ1Þ ¼ rðnÞ aðnÞ ArðnÞð10:105ÞThe equation for rðnÞ in step 1 is needed only to calculate rð0Þ , while Eq. (10.105)can be used afterwards. With this formulation there will be no need to computeA/ðnÞ , as it was replaced by ArðnÞ . However a shortcoming of this approach, is thelack of feedback from the value of /ðnÞ into the residual, which may cause thesolution to converge to a value different from the exact one due to accumulation ofroundoff errors.
This deficiency can be resolved by periodically computing theresidual using the original equation.10.3.13 The Conjugate Gradient MethodWhile the steepest descent method guarantees convergence, its rate of convergenceis low. This slow convergence is caused by oscillations around local minimaforcing the method to search in the same direction repeatedly. To avoid thisundesirable behavior every new search should be in a direction different from thedirections of previous searches [19].
This can be accomplished by selecting a set ofsearch directions dð0Þ ; dð1Þ ; dð2Þ ; . . .; dðN1Þ that are A-orthogonal. Two vectors dðnÞand dðmÞ are said to be A-orthogonal if they satisfy the following condition: TdðnÞ AdðmÞ ¼ 0ð10:106ÞIf in each search direction the right step size is taken, the solution will be foundafter N steps. Step n + 1 is chosen such that/ðnþ1Þ ¼ /ðnÞ þ aðnÞ dðnÞð10:107ÞSubtracting ϕ from both sides of the above equation, an equation for the error isobtained aseðnþ1Þ ¼ eðnÞ þ aðnÞ dðnÞð10:108ÞCombining Eq. (10.100) with Eq.
(10.108), an equation for the residual is found asrðnþ1Þ ¼ Aeðnþ1Þ¼ A eðnÞ þ aðnÞ dðnÞ¼ rðnÞ aðnÞ AdðnÞð10:109Þ33810Solving the System of Algebraic EquationsEquation (10.109) shows that each new residual rðnþ1Þ is just a linear combinationof the previous residual and AdðnÞ .It is further required that eðnþ1Þ be A-orthogonal to dðnÞ . This new condition isequivalent to finding the minimum point along the search direction dðnÞ .
Using thisA-orthogonality condition between eðnþ1Þ and dðnÞ along with Eq. (10.108) anexpression for aðnÞ can be derived asdðnÞT T Aeðnþ1Þ ¼ 0 ) dðnÞ A eðnÞ þ aðnÞ dðnÞ ðnÞ T ðnÞdrðnÞ¼ 0 ) a ¼ TðnÞdAdðnÞð10:110ÞThe above requirement also implies that T TdðnÞ Aeðnþ1Þ ¼ 0 ) dðnÞ rðnþ1Þ ¼ 0ð10:111ÞIf the search directions are known then aðnÞ can be calculated.To derive the search direction, it is assumed to be governed by an equation of theformdðnþ1Þ ¼ rðnþ1Þ þ bðnÞ dðnÞð10:112ÞThe A-orthogonality requirement of the d vectors implies thatTdðnþ1Þ AdðnÞ ¼ 0ð10:113ÞSubstituting the value of dðnþ1Þ from Eq.
(10.112) in Eq. (10.113) yieldsbðnÞTrðnþ1Þ AdðnÞ¼ TdðnÞ AdðnÞð10:114ÞFrom Eq. (10.109) an expression for AdðnÞ is obtained asAdðnÞ ¼ 1 ðnþ1Þr rðnÞðnÞað10:115Þ10.3Iterative Methods339Combining Eqs. (10.110), (10.114), and (10.115) leads to ðnþ1Þ T ðnþ1Þrr rðnÞðnÞb ¼ ðnÞ TrðnÞd ðnþ1Þ T ðnþ1Þ ðnþ1Þ T ðnÞrr rr|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}¼0¼ ðnÞ TðnÞdr ðnþ1Þ T ðnþ1Þrr¼ TdðnÞ rðnÞThe denominator of the above equation can be further expressed asTTdðnÞ rðnÞ ¼ rðnÞ þ bðn1Þ dðn1Þ rðnÞ TT¼ rðnÞ rðnÞ þ bðn1Þ dðn1Þ rðnÞ|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}¼ rðnÞð10:116Þð10:117Þ¼0TrðnÞUsing Eqs. (10.116) and (10.117), the final expression for bðnÞ is obtained as ðnþ1Þ T ðnþ1ÞrrðnÞð10:118Þb ¼TðrðnÞ Þ rðnÞThe Conjugate Gradient algorithm becomesd(0) = r(0) = b A(0)(choose residual as starting direction)iterate starting at(n)until convergence( n )T ( n )d r(Choose factor in d direction)d( n )T Ad ( n )( n+1)= ( n ) + ( n )d( n ) (Obtain new )(n)=(n)r ( n+1) = r ( n )Ad ( n ) (calculate new residual)r ( n+1)T r ( n+1)(n)= ( n )T ( n ) (Calculate coefficient to conjugate residual)r rd( n+1) = r ( n+1) + ( n )d( n ) (obtain new conjugated search direction)The convergence rate of the CG method may be increased by preconditioning.This can be done by multiplying the original system of equations by the inverse of34010Solving the System of Algebraic Equationsthe preconditioned matrix P1 , where P is a symmetric positive-definite matrix, toyield Eq.
(10.79). The problem is that P1 A is not necessarily symmetric even if Pand A are symmetric. To circumvent this problem the Cholesky decomposition isused to write P in the formP ¼ LLTð10:119ÞTo guarantee symmetry, the system of equations is written asL1 ALT LT / ¼ L1 bð10:120Þwhere L1 ALT is symmetric and positive-definite. The CG method can be used tosolve for LT /, from which ϕ is found.
However, by variable substitutions, L can beeliminated from the equations without disturbing symmetry or affecting the validityof the method. Performing this step and adopting the terminology used with the CGmethod, the various steps in the preconditioned CG method are obtained.The preconditioned CG method can be summarized as follows:r(0) = b Ad( 0 ) = P 1r ( 0 ) (choose starting direction)iterate starting at ( n ) until convergenceTn( n+1)rand( r( ) ) P r( ) (Choose factor in=( d( ) ) Ad( )n(n)(0)( n+1)=(n)=r(n)T+n(n) (n)d (Obtain new(n)( r( ) )=( r( ) )n+1( n+1)n1nd direction))Ad ( n ) (calculate new residual)TP 1r ( n+1)TP 1r ( n )d( n+1) = P 1r ( n+1) +(Calculate coefficient to conjugate residual)( n+1) ( n )d(obtain new conjugated search direction)Many pre-conditioners have been developed with a wide spectrum of sophistication varying from a simple diagonal matrix whose elements are the diagonalelements of the original matrix A (Jacobi pre-conditioner) to more involved onesusing incomplete Cholesky factorization.
Nonetheless, the CG method shouldalways be used with a pre-conditioner when solving large systems of equations.10.3.14 The Bi-conjugate Gradient Method (BiCG)and Preconditioned BICGThe matrix of coefficients resulting from the discretization of the diffusion equationpresented in Chap. 8 and some other equations like the incompressible pressure or10.3Iterative Methods341pressure correction equation that will be presented in Chap.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
