Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 49
Текст из файла (страница 49)
This procedure is also known asLU factorization. In addition, direct methods benefiting from the banded structureof A, applicable to structured grid methods, will be discussed.10.2.1 Gauss EliminationThe best way to describe the Gauss elimination technique is to start with a simpleexample. For that purpose a linear system of equations in the two unknowns /1 and/2 is considered.
The equations are given bya11 /1 þ a12 /2 ¼ b1ð10:3Þa21 /1 þ a22 /2 ¼ b2ð10:4ÞThe system can be solved by eliminating one of the variables from one of theequations [say /1 from Eq. (10.4)]. This can be done by multiplying Eq. (10.3) bya21 =a11 and subtracting the resulting equation from Eq. (10.4). This yieldsa22 a21a21a12 /2 ¼ b2 b1a11a11ð10:5ÞThe obtained equation involves only one unknown.
As such it can be used to solvefor /2 , which is obtained asa21b1a11/2 ¼a12 a21a22 a11b2 ð10:6ÞKnowing /2 , its value can be substituted back into Eq. (10.3) to find /1 .Performing this step, /1 is found to be30610Solving the System of Algebraic Equationsa21b bb1 a12 2 a11 1/1 ¼a11 a11 a22 a12 a21a11ð10:7ÞThe above procedure is composed of two steps.
In the first step the equations aremanipulated in order to eliminate one of the unknowns. The end result of this step isan equation with one unknown. In the second step, this equation is solved directlyand the result back-substituted into one of the equations to solve for the remainingunknown. The same procedure can be generalized to a system of N equationsdescribed by Eq. (10.1) or (10.2), as detailed next.10.2.2 Forward EliminationIn the derivations to follow, the first row of A refers to the discretized equation for/1 , the second row represents the equation for /2 , and in general the ith row refersto the equation for /i . The procedure starts by eliminating /1 from all equationsbelow row 1 in A.
To eliminate /1 from the ith row (i = 2, 3, …, N), the coefficientsof the first row are multiplied by ai1 =a11 and the resulting equation is subtractedfrom the ith row. The system of equations by the end of this step becomes2 332 /1 32b1a11 a12 . . . a1N1 a1N6 0 7767 / 7 6b 7626 276 0 a022 . . . a02N1 a02N 7676 . 76.. 77¼6ð10:8Þ.76 .... 76..... 7 6766 . 76 .......7654 ..
5 6 .. 744 . 5.0 a0N2 . . . a0NN1 a0NNb0N/NThen /2 is eliminated from all equations below row 2 in the modified A. Toeliminate /2 from the ith row (i = 3, 4, …, N), the coefficients of the second row aremultiplied by a0i2 =a022 and the resulting equation is subtracted from the ith row. Then/3 is eliminated from all rows below the third row in the modified coefficient matrixand the process is continued until /N1 is eliminated from the Nth row leading tothe following equivalent system of equations with its matrix A transformed into anupper triangular matrix:323 232a11 a12 a13 . .
. a1N1 a1Nb1/10000767 6766 0 a22 a23 . . . a2N1 a2N 76 /2 7 6 b02 77776660000000 a33 . . . a3N1 a3N 76 /3 7 ¼ 6 b003 76 0ð10:9Þ776 .66.. 76 . 7 6 . 7......76 ...4 .. 54 . 5 4 . 5. .....N1/NbNN1000 ...0aNNThe resulting algorithm is as described below.10.2Direct or Gauss Elimination Method30710.2.3 Forward Elimination AlgorithmFor k = 1 to N 1{For i = k + 1 to N{Ratio =aikakk{For j = k + 1 to Naij = aijRatio * akj}bi = bi}Ratio *bk}10.2.4 Backward SubstitutionThe modified system of equations given by Eq.
(10.9) indicates that the onlyunknown in the Nth equation is /N . Therefore this equation can be used to obtain/N as/N ¼bNN1N1aNNð10:10ÞThe (N − 1)th equation is function of /N1 and /N . Having found /N then thisequation can be used to find /N1 as/N1 ¼N2N2bN1 aN1N/NN2aN1N1ð10:11ÞThe process continues moving backward and by the time the ith equation isreached, the values /iþ1 ; /iþ2 ; /iþ3 ; . . .; /N1 ; /N would have become availableand as such /i can be computed using30810/i ¼bi1iSolving the System of Algebraic EquationsPNj¼iþ1ai1iiai1ij /jð10:12Þthe process is continued until /1 is calculated. Algorithmically, this is representedas shown below.10.2.5 Back Substitution AlgorithmN=bNaNNFor i = N 1 to 1{Term = 0{For j = i + 1 to NTerm = Term + aij *j}i=iTermaii}Techniques to improve the performance of the method to avoid division by zerothrough pivoting (interchanging rows in order to select the largest pivoting element)and to reduce roundoff errors in large systems are available but are not discussedhere.
Interested readers may consult specialized textbooks on the subject [1–4]. Thepresented algorithm shows that the method is expensive and the number of operations required to solve a linear system of N equations is proportional to N3/3 ofwhich only N2/2 arithmetic operations are required for back substitution. This highcomputational cost has enticed researchers to look for more efficient specializedsolvers for systems with sparse matrices.10.2.6 LU DecompositionAnother direct method for solving linear algebraic systems of equations is the LU ormore generally the PLU (where P refers to the pivoting process mentioned above),in this book reference is made only to LU factorization method, which are variants10.2Direct or Gauss Elimination Method309of the Gauss elimination method.
The advantage of these methods over the Gausselimination method is that once the (P)LU factorization is performed the linearsystem can be solved as many times as needed for different values of the right handside vector b without performing any additional elimination, which would still berequired with the Gauss method.Based on the elimination performed in the previous section, Eq. (10.1) wastransformed into an upper triangular matrix as given by Eq. (10.9), which can bewritten as2u116 066 066 ..4 .0u12u220...u13u23u33...00.
. . u1N1. . . u2N1. . . u3N1.........u1Nu2Nu3N...0323 2 3/1c176 /2 7 6 c2 7767 6 776 /3 7 6 c3 7767¼6 776 .. 7 6 .. 754 . 5 4 . 5uNN/Nð10:13ÞcNUsing compact matrix notation Eq. (10.13) can be simplified toU/ c ¼ 0ð10:14ÞLet L be a unit lower triangular matrix (diagonal elements are set to 1 in order tomake the factorization unique) given by216 ‘2166L ¼ 6 ‘316 ..4 .‘N101‘32...‘N2001...‘N3...............000...‘NN130077077.. 7.5ð10:15Þ1such that if Eq. (10.14) is multiplied by L, Eq. (10.1) is recovered.
If this holds,then the following can be written:LðU/ cÞ ¼ LU/ Lc ¼ A/ bð10:16ÞBased on matrix properties it follows thatLU ¼ Að10:17ÞLc ¼ bð10:18ÞandEquation (10.17) indicates that A is written as the product of a proper lower and anupper triangular matrix, known as LU factorization.31010Solving the System of Algebraic Equations10.2.7 The Decomposition StepThe efficient procedure to find the L and U coefficients described next is denoted bythe Crout decomposition [1–4]. In the original Crout algorithm a unit upper triangular matrix is used whereas here a unit lower triangular matrix is assumed. Theprocedure is based on multiplying L and U to obtain A such that2166 ‘2166 ‘3166 .6 .4 .‘N1232u110760 76 07660776 07..
766 .. 54 ..0100......00‘32...1.........0...‘N2‘N3...‘NN11a11a12a13...a1N1a22a23...a2N1a32...aN2a33...aN3...66 a2166¼ 6 a316 .6 .4 .aN10a3N1....... . . aNN1a1Nu12u22u13u23......u1N1u2N10...u33.........u3N1...030...03u1N7u2N 77u3N 77.. 77. 5uNN7a2N 77a3N 77..
77. 5aNNð10:19ÞThe calculation of the coefficients starts by multiplying the first row of L by allcolumns of U, and equating with the corresponding coefficients of A to yieldu1j ¼ a1jj ¼ 1; 2; 3; . . .; Nð10:20ÞThen the second through Nth rows of L are multiplied by the first column ofU leading to‘i1 u11 ¼ ai1 ) ‘i1 ¼ai1u11i ¼ 2; 3; . . .; Nð10:21ÞThe process is repeated by multiplying the second row of L by the second throughNth columns of U to giveu2j ¼ a2j ‘21 u1jj ¼ 2; 3; . .
.; Nð10:22Þafter that the third through Nth rows of L are multiplied by the second column ofU to give‘i2 u22 þ ‘i1 u12 ¼ ai2 ) ‘i2 ¼ai2 ‘i1 u12u22i ¼ 3; 4; . . .; Nð10:23Þ10.2Direct or Gauss Elimination Method311In general, the ith row of L is multiplied by the ith through Nth columns of U,resulting inuij ¼ aij i1X‘ik ukjj ¼ i; i þ 1; .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
