Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 50
Текст из файла (страница 50)
. .; Nð10:24Þk¼1and the (i + 1)th through Nth rows of L are multiplied by the ith column of U,giving‘ki ¼aki Pi1j¼1 ‘kj ujiuiik ¼ i þ 1; i þ 2; . . .; Nð10:25ÞFor the Nth row of L, its coefficients are multiplied by the coefficients of the Nthcolumn of U from which uNN is obtained asuNN ¼ aNN N1X‘Nk ukNk¼1A summary of the LU factorization is shown algorithmically below.10.2.8 LU Decomposition Algorithmu1 j = a1 j j = 1 to Nai1i = 2 to Nu11For i = 2 to N 1i1={uij = aiji 1uik kjj = i,i + 1,..., Nk=1i 1akiki=kju jik = i + 1,i + 2,..., Nj=1uii}u NN = aNNN 1uNi iNi=1ð10:26Þ31210Solving the System of Algebraic Equations10.2.9 The Substitution StepHaving decomposed the original matrix A into L and U, the system of equationscan be solved in a two step procedure via Eqs.
(10.18) and (10.14). Note that thetwo-step procedure is equivalent to solving two linear systems of equations but nowsimplified by the fact that L and U are of lower and upper triangular form,respectively.In the first step the vector c is obtained from Eq. (10.18) by forward substitution. The process can be described asc 1 ¼ b1c i ¼ bi i1X‘ij cji ¼ 2; 3; .
. .; Nð10:27Þj¼1In the second step, the / values are found from Eq. (10.14) by back substitution.The process is described by/N ¼/i ¼cNuNNPci Nj¼iþ1 uij /juiið10:28Þi ¼ N 1; N 2; . . .; 3; 2; 1The elements of L and U can be directly stored in the original matrix A if it is nolonger needed. This is because the elements of A are only needed when the corresponding elements of either L or U are calculated. The number of operationsrequired to perform the LU factorization of a square matrix of size N × N is 2N3/3,which is double the number of operations required to solve the same system ofequations by Gauss elimination. Again the advantage of using LU factorization iswhen the same matrix A applies to many systems with different b vectors.Nevertheless, the main reason for introducing the LU factorization is because itforms the basis for developing some of the more efficient iterative solvers of linearalgebraic systems of equations, which will be introduced in the next section.10.2.10 LU Decomposition and Gauss EliminationIt may not be apparent, but Gauss elimination can be used to perform LUdecomposition.
It was shown that the forward elimination step results in an uppertriangular matrix U. In the process however, L is actually produced. The elementsof L are the factors (denoted by ratio in the Gauss elimination algorithm) by whichthe rows are multiplied during the various elimination steps.
The below algorithm,which assumes a unit lower triangular matrix L, performs the LU decomposition ofA by Gauss elimination.10.2Direct or Gauss Elimination Method31310.2.11 LU Decomposition Algorithm by Gauss Eliminationu1 j = a1 j j = 1 to NFor k = 1 to N 1{For i = k + 1 to N{ik=aikakk{For j = k + 1 to Nuij = aij}}ik* akj}Example 1Solve the following system of linear algebraic equations using the LUdecomposition technique:23 16 2 664 0 2000162332 3 230/176 7 60 776 /2 7 ¼ 6 4 71 54 /3 5 4 5 537/4SolutionThe system is of the form A/ ¼ b. The L and U should satisfy216‘6 21LU ¼ 64 ‘31‘410100‘32‘421‘43320u116 00776760 54 010u12u22u13u2300u3303 23u146 2u24 77 67¼6u34 5 4 0u44016012062300 7771 5731410Solving the System of Algebraic EquationsFollowing the procedure described above the elements are calculated asfollows:8u11 ¼ 3>><u12 ¼ 1u1j ¼ a1j j ¼ 1; 2; 3; 4 )u ¼0>>: 13u14 ¼ 08a212>>¼> ‘21 ¼>3u>11><ai1a31‘i1 ¼i ¼ 2; 3; 4 ) ‘31 ¼¼0>u11u11>>>a>>: ‘41 ¼ 41 ¼ 0u11816>< u22 ¼ a22 ‘21 u12 ¼3u2j ¼ a2j ‘21 u1j j ¼ 2; 3; 4 ) u ¼ a ‘ u ¼ 12321 13>: 23u24 ¼ a24 ‘21 u14 ¼ 08< ‘ ¼ a32 ‘31 u12 ¼ 3ai2 ‘i1 u1232u228‘i2 ¼i ¼ 3; 4 ): ‘ ¼ a42 ‘41 u12 ¼ 0u2242uij ¼ aij i1X‘ik ukj i ¼ 3; j ¼ 3; 4 )k¼1‘ki ¼aki 8<:u34 ¼ a34 ‘31 u14 ‘32 u24i ¼ 3; k ¼ 4 ) ‘43 ¼uiiuNN ¼ aNN N1Xa43 ‘41 u13 ‘42 u2316¼45u33‘Nk ukN ) u44 ¼ a44 ‘41 u14 ‘42 u24 ‘43 u34 ¼k¼1Therefore the L and U matrices are given by216 266 36L¼66 06400138000370777710775161450458¼ 1u33 ¼ a33 ‘31 u13 ‘32 u23 ¼Pi1j¼1 ‘kj ujiu222366066U¼6606401163000145800370 77771 77299 5452994510.2Direct or Gauss Elimination Method315The c vector should satisfy216 266 36Lc ¼ b ) 66 06400003372 3 23c10776 7 676 c2 7 6 4 77¼710 74 c3 5 4 5 5735 c41614513800with its solution obtained asc1 ¼ 32 c1 þ c2 ¼ 4 ) c2 ¼ 63329 c2 þ c3 ¼ 5 ) c3 ¼841619 c3 þ c4 ¼ 3 ) c4 ¼ 4545The solution to the original equation is obtained by solving2366066U/ ¼ c ) 6606401163000145803272 3/60 776 1 7 676 /2 7 674 5 ¼ 661 7 /374299 5 /4045336 729 7774 751945with the solution found as9324352991919 >>/4 ¼ ) /4 ¼ >76>4545299 >6 299 7>>>76>4529382 >6 408 7>/ /4 ¼) /3 ¼=768 342996 299 7)/¼676 382 716408 >>>76/2 /3 ¼ 6 ) /2 ¼>>763299 >>6 299 7>>54435 >19>;3/1 /2 ¼ 3 ) /1 ¼29929910.2.12 Direct Methods for Banded Sparse MatricesThe Gauss elimination and LU decomposition methods are applicable to any systemof equations.
In specific it can be used for solving the system of equations resultingfrom the discretization of the conservation equations of interest in this book on31610Solving the System of Algebraic Equationsstructured or unstructured grid networks. When a structured grid method is used thediscretization process results in a system of equations with the non-zero elements ofits matrix of coefficients aligning along few diagonals. Depending on the discretization stencil used and the dimension of the problem being solved, tridiagonal orpentadiagonal matrices may arise, for which efficient algorithms have been developed as described next.10.2.13 TriDiagonal Matrix Algorithm (TDMA)The TriDiagonal Matrix Algorithm (TDMA), also known as Thomas algorithm[5, 6], solves the system of algebraic equations with a tridiagonal coefficient matrixwritten asai /i þ bi /iþ1 þ ci /i1 ¼ dii ¼ 1; 2; 3; .
. .; Nc 1 ¼ bN ¼ 0ð10:29ÞFor the grid arrangement adopted in this book, i refers to the grid point locationshown in Fig. 10.1.i=1 234N2N 1 NFig. 10.1 One dimensional grid arrangementFor i = 1 the equation can be used to solve for /1 in term of /2 asi ¼ 1 ) a1 /1 ¼ b1 /2 þ d1 ) /1 ¼ b1d1/ þa1 2 a1ð10:30ÞSimilarly for i = 2 Eq. (10.29) with the help of Eq.
(10.30) allows expressing /2solely in term of /3 asi ¼ 2 ) a2 /2 ¼ b2 /3 c2 /1 þ d2 ) /2 ¼ a1 b2d2 a1 c 2 d1/3 þa1 a2 c 2 b1a1 a2 c 2 b1ð10:31ÞThe same can be repeated for /3 through /N suggesting that in general /i can beexpressed as function of /iþ1 according to/i ¼ Pi /iþ1 þ Qii ¼ 1; 2; 3; . . .; NEquation (10.32) for i − 1 when combined with Eq. (10.29) results inð10:32Þ10.2Direct or Gauss Elimination Method/i1 ¼ Pi1 /i þ Qi1ai /i þ bi /iþ1 þ ci /i1 ¼ di317) /i ¼ bidi ci Qi1/ þð10:33Þai þ ci Pi1 iþ1 ai þ ci Pi1Comparing Eq. (10.32) with Eq. (10.33) the following recurrence relations for Piand Qi are found:Pi ¼ biai þ ci Pi1Qi ¼di ci Qi1ai þ ci Pi1i ¼ 1; 2; . .
.; Nð10:34ÞFor i = 1 the values for P1 and Q1 are computed from Eq. (10.30) asP1 ¼ b1a1Q1 ¼d1a1ð10:35ÞFor i = N, since bN = 0 the following is deduced:bN ¼ 0 ) PN ¼ 0 ) /N ¼ QNð10:36ÞThe TDMA solution algorithm can be summarized as follows:1. Compute the values for P1 and Q1 using Eq. (10.35)i = 2, 3,..., N use forward recursion to compute the values of Pi2. Forand Qi from Eq. (10.34)3. Set N = QN as given by Eq. (10.36)4. For i = N 1, N 2,...3,2,1 use backward recursion to compute the valuesofifrom Eq.
(10.32)10.2.14 PentaDiagonal Matrix Algorithm (PDMA)The PentaDiagonal Matrix Algorithm (PDMA) [7–10], solves the system of algebraicequations with a pentadiagonal coefficient matrix arising from discretization schemesthat relate the value of /i at grid point i to the values at its two upstream (i − 1 andi − 2) and two downstream (i + 1 and i + 2) neighboring node values.
For the notationillustrated schematically in Fig. 10.1, the general algebraic equation is written asai /i þ bi /iþ2 þ ci /iþ1 þ di /i1 þ ei /i2 ¼ fii ¼ 1; 2; 3; . . .; Nð10:37ÞSubject tod1 ¼ e 1 ¼ e 2 ¼ 0bN1 ¼ bN ¼ cN ¼ 0ð10:38Þ31810Solving the System of Algebraic EquationsFor i = 1, Eq. (10.37) gives/1 ¼ b1c1f1/ / þa1 3 a1 2 a1ð10:39Þwhile for i = 2 the value of /2 is found to be/2 ¼ a1 b2a1 c 2 b1 d2a1 f2 d2 f1/ / þa1 a2 d2 c 1 4 a1 a2 d2 c 1 3 a 1 a2 d2 c 1ð10:40ÞThe process can be continued for other values of i and in general /i can beexpressed as/i ¼ Pi /iþ2 þ Qi /iþ1 þ Rii ¼ 1; 2; 3; .
. .; Nð10:41ÞComputing /i1 and /i2 using Eq. (10.41) and substituting their values inEq. (10.37), an equation for /i is derived asbi/ai þ ei Pi2 þ ðdi þ ei Qi2 ÞQi1 iþ2ci þ ðdi þ ei Qi2 ÞPi1/ai þ ei Pi2 þ ðdi þ ei Qi2 ÞQi1 iþ1fi ei Ri2 ðdi þ ei Qi2 ÞRi1þai þ ei Pi2 þ ðdi þ ei Qi2 ÞQi1/i ¼ ð10:42ÞComparing Eqs. (10.41) and (10.42) Pi, Qi, and Ri are found asbiai þ ei Pi2 þ ðdi þ ei Qi2 ÞQi1ci þ ðdi þ ei Qi2 ÞPi1Qi ¼ ai þ ei Pi2 þ ðdi þ ei Qi2 ÞQi1fi ei Ri2 ðdi þ ei Qi2 ÞRi1Ri ¼ai þ ei Pi2 þ ðdi þ ei Qi2 ÞQi1Pi ¼ ð10:43Þwith their values for i = 1 and 2 given byP1 ¼ b1a1P2 ¼ b2a2 þ d2 Q 1Q1 ¼ c1a1f1a1c2 þ d2 P1Q2 ¼ a2 þ d2 Q 1R1 ¼R2 ¼f2 d2 R1a2 þ d2 Q 1ð10:44ÞSince bN1 ¼ bN ¼ cN ¼ 0 then PN1 ¼ PN ¼ QN ¼ 0. Thus, the equations for/N1 and /N are found from10.2Direct or Gauss Elimination Method/N ¼ RN/N1 ¼ QN1 /N þ RN1319ð10:45ÞThe PDMA solution algorithm can be summarized as follows:R1 , P2 , Q2 , and R2 using Eq.1.
Compute the values for P1 , Q1 ,(10.44).2. For i = 3, 4,..., N use forward recursion to compute the values of Pi ,Qi , and Ri from Eq. (10.43).3. Compute4. For i = NiNandN 1from Eq. (10.45).2,...3,2,1 use backward recursion to compute the values offrom Eq. (10.41).10.3 Iterative MethodsDirect methods are generally not appropriate for solving large systems of equationsparticularly when the coefficient matrix is sparse, i.e., when most of the matrixelements are zero. This is more so when the linearized system of equations isnonlinear with solution dependent coefficients, or when dealing with time dependent problems.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
