Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 51
Текст из файла (страница 51)
This is exactly the type of equations encountered when solving fluidflow problems.In contrast, iterative methods are more appealing for these problems since thesolution of the linearized system becomes part of the iterative solution process. Addto that the low computer storage and low computational cost requirements of thisapproach relative to the direct method.There are many families of iterative methods and for a thorough review of thisapproach the reader is directed to dedicated books on the topic [11–14].
In thischapter a brief examination of basic iterative methods is provided along with anappraisal of multigrid algorithms that are generally used to address their deficiency.The Gauss elimination and LU decomposition direct methods were introduced forthe sole purpose of clarifying some fundamental numerical processes needed forunderstanding iterative methods.To unify the presentation of these methods, the coefficient matrix will be writtenin the following form:A¼DþLþUð10:46Þwhere D, L, and U refers to a diagonal, strictly lower, and strictly upper matrix,respectively.Iterative methods for solving a linear system of the type A/ ¼ b, compute aseries of solutions /ðnÞ that, if certain conditions are satisfied, converge to the exactsolution /.
Thus, for the solution, a starting point is chosen (i.e., /ð0Þ is selected as32010Solving the System of Algebraic Equationsthe initial condition or initial guess) and an iterative procedure that computes /ðnÞfrom the previously computed /ðn1Þ field is developed.A “fixed-point” iteration can always be associated to the above system bydecomposing matrix A asA¼MNð10:47ÞUsing this decomposition, Eq. (10.1) is rewritten asðM NÞ/ ¼ bð10:48ÞApplying a fixed point iteration solution procedure, Eq.
(10.48) becomesM/ðnÞ ¼ N/ðn1Þ þ bð10:49Þwhich can be rewritten in the following form:/ðnÞ ¼ B/ðn1Þ þ Cbn ¼ 1; 2; . . .ð10:50Þwhere B ¼ M1 N and C ¼ M1 . Different choices of these matrices define different iterative methods.Before embarking on the description of the various iterative methods, a minimalset of characteristics that an iterative method should possess to guarantee convergence is first presented.A.
The iterative equation can be written at convergence as/ ¼ B/ þ Cbð10:51Þwhich, after rearranging, becomesC1 ðI BÞ/ ¼ bð10:52ÞComparing Eq. (10.52) with Eq. (10.1), the coefficient matrix is obtained asA ¼ C1 ðI BÞð10:53ÞB þ CA ¼ Ið10:54Þor, alternatively, asThis relation between the various matrices ensures that once the exact solution isreached all consecutive iterations will not modify it.B. Starting from some guess /ð0Þ 6¼ /, the method should guarantee that /ðnÞ willconverge to / as n increases.
Since /ðnÞ can be expressed in terms of /ð0Þ as10.3Iterative Methods321/ðnÞ ¼ Bn /ð0Þ þn1XBi Cbð10:55Þi¼0then, for the above to be true, B should satisfyB B ffl{zfflfflfflfflfflfflfflfflfflfflfflfflB Blim Bn ¼ lim |fflfflfflfflfflfflfflfflfflfflfflfflffl} ¼ 0n!1n!1ð10:56Þn timesEquation (10.56) implies that the spectral radius of B should be less than 1, i.e.,qðBÞ\1ð10:57ÞThis condition guaranties that the iterative method is self corrective, i.e., it isrobust to any error adversely inserted into the solution vector /.More insight into the above condition can be obtained by defining the error e(n)in the solution as the difference between the exact value and the value at anyiteration (n), theneðnÞ ¼ /ðnÞ /andeðn1Þ ¼ /ðn1Þ /ð10:58ÞSubtracting Eq. (10.51) from Eq. (10.50) and using the definitions inEq.
(10.58), a relation between the error at iteration n and n − 1 is obtained aseðnÞ ¼ Beðn1Þð10:59ÞThus, for the method to converge, the following should be satisfied:lim eðnÞ ¼ 0n!1ð10:60ÞTo translate Eq. (10.60) into something meaningful, the eigenvectors of B areassumed to be complete and to form a full set, meaning that they form a basis forRN. This being the case then e can be expressed as a linear combination of theN eigenvectors v of B.
That ise¼NXai við10:61Þi¼1with each of the eigenvectors satisfyingBvi ¼ ki við10:62Þwhere ki is the eigenvalue corresponding to the eigenvector vi. Starting with thefirst iteration, Eq. (10.59) gives32210eð1Þ ¼ Beð0Þ ¼ BNXai vi ¼Solving the System of Algebraic EquationsNXi¼1NXai ðBvi Þ ¼i¼1ai ki við10:63Þi¼1For the second iteration, the error is obtained aseð2Þ ¼ Beð1Þ ¼ BNXai ki vi ¼i¼1NXai ki ðBvi Þ ¼i¼1NXai k2i við10:64Þi¼1This procedure can be continued and it is easily shown by induction thateðnÞ ¼NXai kni við10:65Þi¼1Therefore for the iterative procedure to converge as n approaches infinity, alleigenvalues should be less than 1.
If any of them is greater than 1 then the errorwill tend to infinity. This explains the importance of the spectral radius ρ of thematrix B defined asNqðBÞ ¼ maxðki Þð10:66Þi¼1that was mentioned above. The convergence of iterative methods is acceleratedby reducing the spectral radius of the iterative matrix. This is at the heart ofiterative techniques.C. Some type of a stopping criterion is needed with iterative methods.
Many usedcriteria are based on a variation of the norm of the residual error defined asrðnÞ ¼ A/ðnÞ bð10:67ÞOne criterion is to find the maximum residual in the domain and to require itsvalue to become less than some threshold ε to declare a solution converged, i.e.,NXðnÞMaxbi aij /j ei¼1Nð10:68Þj¼1or that the root mean square residual be smaller than ε, i.e.,PN i¼1bi PNðnÞj¼1 aij /jN2eð10:69ÞAnother possible criterion is for the maximum normalized difference betweentwo consecutive iterations to drop below e. This condition can be written as10.3Iterative Methods323/ðnÞ /ðn1Þ iiMax 100 eðnÞi¼1 /Nð10:70Þi10.3.1 Jacobi MethodPerhaps the simplest of the iterative methods for solving a linear system of equations is the Jacobi method, which is graphically presented in Fig. 10.2.φ (n)=φ (n−1)-Fig.
10.2 A graphical representation of the Jacobi methodConsidering the system of equations described by Eq. (10.1), if the diagonalelements are nonzero, then the first equation can be used to solve for /1 , the secondequation to solve for /2 , and so on. The solution process starts by assigningguessed values to the unknown vector /. These guessed values are used to calculatenew estimates starting with /1 , then /2 , and computations proceed until a newestimate for /N is computed. This represents one iteration. Results obtained aretreated as a new guess for the next iteration and the solution process is repeated.Iterations continue until the changes in the predictions between two consecutiveiterations drop below a vanishing value or until a preset convergence criterion issatisfied.
Once this happens the final solution is reached. In this method, given somecurrent estimate /ðn1Þ , an update is obtained using the following relation:10ðnÞ/jCBNX1Bðn1Þ CC i ¼ 1; 2; 3; . . .; NB¼ B bi aij /jCaii @Aj¼1j6¼ið10:71Þ32410Solving the System of Algebraic EquationsEquation (10.71) indicates that values obtained during an iteration are not used inthe subsequent calculations during the same iteration but rather retained for the nextiteration. Using matrices, the expanded form of Eq.
(10.71) is given by2a1166 066 ..6 .664 0020a22............00......00...aN1;N100a21...6666þ6664 aN1;1aN1a120...00.........32/1376 / 776 2 776776 .. 776 . 77677670 54 /N1 5aNN/NaN1;2......a1N1a2N1...0aN2aN3...3b176 / 7 6 b2 776 2 7 6 77 6 77676 .. 7 6 .. 776 . 7 ¼ 6 . 77 6 7767 6 . 776aN1;N 54 /N1 5 4 .. 50/NbNa1Na2N...32/132ð10:72ÞSolving for /ðnÞ , Eq. (10.72) yields2331a11 0 . . .006 ðnÞ 7 676 / 7 6 0 a22 .
. .00 76 2 7 677 6 .6.. 7....6 .. 7 ¼ 6 .. 7. ....6 . 7 6 .77 6676 /ðnÞ 7 4 0 . . . 0 aN1;N1054 N1 5ðnÞ00...0aNN/N02 3 2b10a12. . . a1N1B6 7 6B6 b2 7 6 a210. . . a2N1B6 7 6B6 . 7 6 ....B6 .. 7 6 .......B6 7 6 .B6 . 7 6B4 .. 5 4 aN1;1 aN1;2 . . .0@aN1aN2aN3...bNðnÞ/1232 ðn1Þ 31/17C76C6a2N 76 /ð2n1Þ 77C767.. 76 . 7C.. 7C. 7C767C76ðn1Þ76aN1;N 54 /N1 5CAðn1Þ0/a1NNð10:73ÞUsing Eqs.
(10.46) and (10.73) can be more concisely written as/ðnÞ ¼ D1 ðL þ UÞ/ðn1Þ þ D1 bð10:74ÞThe Jacobi method converges as long as q D1 ðL þ UÞ \1. This condition issatisfied for a large class of matrices including diagonally dominant ones wheretheir coefficients satisfy10.3Iterative Methods325N Xaij jaii ji ¼ 1; 2; 3; . . .; Nð10:75Þj¼1j6¼i10.3.2 Gauss-Seidel MethodA more popular take on the Jacobi is the Gauss-Seidel method, which has betterconvergence characteristics. It is somewhat less expensive memory-wise since itdoes not require storing the new estimates in a separate array.
Rather, it uses thelatest available estimate of / in its calculations. The iterative formula in the GaussSeidel method, schematically displayed in Fig. 10.3, is given as!i1NXX1ðnÞðnÞðn1Þ/i ¼bi aij /j aij /ji ¼ 1; 2; 3; . . .; Nð10:76Þaiij¼1j¼iþ1In matrix form Eq. (10.76) is written as/ðnÞ ¼ ðD þ LÞ1 U/ðn1Þ þ ðD þ LÞ1 bð10:77ÞIn effect the Gauss-Seidel method uses the most recent values in its iteration,ðnÞspecifically all /j values for j < i since by the time /i is to be calculated, thevalues of /1 ; /2 ; /3 ; .
. .; /i1 at the current iteration are already calculated. Thisapproach also saves memory since the newer value is always overwriting theprevious one. The Gauss-Seidel iterations converge as long asq ðD þ LÞ1 U \1ð10:78ÞAlthough in some cases the Jacobi method converges faster, Gauss-Seidel is thepreferred method.φ (n)=-Fig. 10.3 A graphical representation of the Gauss-Seidel methodφ (n−1)32610Solving the System of Algebraic EquationsExample 2Apply 5 iterations of the Gauss-Seidel and Jacobi methods to the system ofequations in Example 1 and compute the errors at each iteration using theexact solution.435 408 38219:/¼299 299 299299As an initial guess start with the field / ¼ ½ 000 0SolutionDenoting with a superscript ð Þ values from the previous iteration, theequations to be solved in the Jacobi method are as follows:1 / þ33 21/2 ¼ 2/1 þ /3 þ 461/3 ¼ 2/2 þ /4 þ 561/4 ¼ 2/3 37with the error given as e ¼ /exact /computed .
The solution for the firstiteration is obtained as/1 ¼1/1 ¼ ð0 þ 3Þ ¼ 1 ) e1 ¼ j1:4548 1j ¼ 0:454831/2 ¼ ð0 þ 0 þ 4Þ ¼ 0:6667 ) e2 ¼ j1:3645 0:6667j ¼ 0:697861/3 ¼ ð0 þ 0 þ 5Þ ¼ 0:8333 ) e3 ¼ j1:2776 0:8333j ¼ 0:444361/4 ¼ ð0 3Þ ¼ 0:4286 ) e4 ¼ j0:06354 þ 0:4286j ¼ 0:36507Computations proceed in the same manner with solution obtained treated asthe new guess. The results for the first five iterations are given in Table 10.1.Table 10.1 Summary of results obtained using the Jacobi iterative methodIter #/1e1/2e2/3e3/4001.454801.364501.27760e40.06354110.45480.66670.69780.83330.4443−0.42860.365021.22220.23261.13890.22570.98410.2935−0.19050.126931.37967.52E−021.23810.12651.18129.64E−02−0.14748.38E−0241.41274.22E−021.32344.11E−021.22155.61E−02−9.11E−022.75E−0251.44111.37E−021.34112.34E−021.25931.83E−02−7.96E−021.6E−0210.3Iterative Methods327Denoting with a superscript ð Þ values from the previous iteration, theequations to be solved in the Gauss-Seidel method are as follows:1 /2 þ 331/2 ¼ 2/1 þ /3 þ 461/3 ¼ 2/2 þ /4 þ 561/4 ¼ ð2/3 3Þ7with the error given as e ¼ /exact /computed .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
