Solving Equations by Iteration

Solving Equations by Iteration Purpose These notes provide an introduction to a set of files that illustrate the details of iteration calculations. Iteration equations We have seen that the finite-volume equations resulting from terms considering convection, diffusion and a source term result in a finite difference equation of the following form. AS φ S + AW φW + AP φ P + AE φ E + AN φ N = bP [1] The finite volume grid on which this equation is derived is shown at the right. The points at the center of the finite volumes, indicated by the N NW symbol •.and marked as the center point, P, nw ne n and it’s four nearest neighbors to the North, South, East, and West are the nodes whose φ W E values are related by equation [1]. Equation P w e [1] applies to a single central node, P. We will have one such equation for the central node of each control volume. This set of equations s has to be solved simultaneously to find the sw se values of φ at each node. In order to define the entire set of nodes we will first examine a S SE simple grid where the geometry is rectangular with uniform grid spacing, Δx and Δy. In this case we can define the grid coordinates xi and yj by the following equations. xi = x0 + i∆x, i = 0,, N x y j = y0 + j∆y, i = 0,, N y [2] With this grid notation, the value of φ at a given location, φ(xi,yj), is denoted as φij. If φij is the center node, P, in equation [1], the North node, located at (xi, yj+1), is denoted as φij+1; the South node, located at (xi, yj-1), is denoted as φij-1; the East node, located at (xi+1, yj), is denoted as φi+1j; finally, the West node, located at (xi-1, yj), is denoted as φi-1j. The various A coefficients may be different at each grid point. To denote the different values of the A coefficients at different points, the following notation is used. Each of the five coefficients in equation [1] is written with the P, N, S, E, W, as a superscript and an ij subscript. This shows that the particular set of coefficients applies to the equation for the central P node located at (xi,yj). With this notation we can rewrite equation [1] as follows. AijS φij −1 + AijW φi −1 j + AijP φij + AijE φi +1 j + AijN φij +1 = bij [3] We have an equation like this for each node in our two-dimensional grid. The set of equations represents a set of linear algebraic equations that we have to solve simultaneously. If we have a problem in which the values of u are specified at all the boundaries (called the Dirichlet boundary conditions) we will have (Nx – 1)(Ny – 1) nodes at which we will have to solve an equation like [3]. We can illustrate a typical set of equations by considering a very small grid with Nx = 6 and Ny = 5. This sample grid is shown in Figure 1. In this case we have (6 – 1)(5 – 1) = 20 nodes in the grid (not counting boundary nodes.) Figure 1. Example of small finite-difference grid. i=0 i=1 i=2 i=3 i=4 i=6 i=6 j=5 φ44 φ33 φ12 φ43 j=4 φ53 φ42 φ11 φ21 j=2 φ51 φ10 j=3 Boundary nodes Computational Molecule j=1 j=0 The nodes related by the typical finite-difference equation are said to form a computational molecule as shown in Figure 1. Although we can write this set of equations in a matrix form, we usually solve the equations by iteration. In doing this we denote the value of φij at a particular iteration, n, as φij(n). In the iteration process, we first have to make an initial guess for each unknown value of φ. This initial guess is iteration zero and we denote the initial guesses as φij(o). A typical initial guess is to take all values of φij(0). = 0. The basic iteration equation is found by solving equation [3] for φij. φij = bij − AijS φij −1 − AijW φi −1 j − AijE φi +1 j − AijN φij +1 AijP [4] We can now show the simplest iteration approaches. All iteration approached find the new values at iteration n in terms of old values at iteration n – 1. The first approach, Jacobi iteration, finds the new value at iteration n in terms of old values at iteration n – 1.This iteration process is written by taking all the φ values on the right side of equation [4] at iteration n – 1. Here we use the superscripts (n) and (n-1) to denote the values of φ at the two iteration levels. φij( n) = bij − AijS φij( n−−11) − AijW φi(−n1−1j ) − AijE φi(+n1−1j ) − AijN φij( n+−11) AijP [5] This iteration approach is known to be slow to converge. A faster alternative, known as GaussSeidel iteration, uses the order in which the variables are found to get newer values for some of the φ values on the right side of equation [4]. We start each set of iterations by finding a new value of φ11 in the lower left corner of the grid. When finding φ11, we will have two known boundary values at φ10 and φ01. Since these boundary values do not change, they represent the values for the new iteration. Once we compute φ11 we can use that value in the computation of φ12. We continue solving for all the values in the first row, φi1, from left to right. For each φi1, we know a “new” φi0 from the lower row that is a boundary condition and a new value of φ11, from the previous calculation. Once we complete the solution of the first row, we then repeat this process for all other rows from bottom to top. In this fashion, when we solve for φ12, we know the boundary value, φ02, and the value φ11 from the new iteration that was found in the solution for the first row. Similarly, when we solve for φ22, we know the value, φ12, that we just found and the value φ21 from the new iteration that was found in the solution for the first row. In general, we will always know φi-1j and φij-1, from boundary values or previous results from the new iteration, when we are solving for φij. If we use the values for φi-1j and φij-1, from the new iteration, in our solution for φij, our iteration form for equation [4] can be written as follows. φij( n ) = bij − AijS φij( n−)1 − AijW φi(−n1) j − AijE φi(+n1−1j ) − AijN φij( n+−11) AijP [6] Gauss-Seidel iterations converge about twice as fast as Jacobi iterations. A significant improvement to the convergence rate for iterations is obtained by an approach known as successive over relaxation or SOR. This approach is based on visualizing the GaussSeidel iterations as providing an update, ∆φij to the old φij value. This is illustrated below where the symbol φij(n,GS) denotes the result from Gauss-Seidel iteration as found from equation [6]. [ φij( n) = φij( n −1) + ∆φij = φij( n −1) + φij( n,GS ) − φij( n −1) ] [7] In SOR one tries to find a relaxation factor, ω, that multiplies the correction, ∆φij so that the iterations will converge faster. This is like saying that is some correction is good, more of the correction is better. Thus SOR uses the following modification of equation [7] where the correction term ∆φij is multiplied by the relaxation factor, ω. [ φij( n) = φij( n −1) + ω∆φij = φij( n −1) + ω φij( n,GS ) − φij( n −1) ] [8] Typically w is a value between 1 and 2 that increases the correction to the old φij value. This value of ω is said to provide overrelaxation that makes convergence faster. Alternatively, the value of ω may be between 0 and 1. The use of such a value of ω is called underrelaxation. It slows down the convergence in cases where application of the full correction factor will lead to divergence of the iterations. Various approaches are used to estimate the best value of ω to improve convergence. For the finite difference version of Laplace’s equation (shown below) on a square grid with uniform spacing of N = Nx = Ny nodes in each coordinate direction, the optimum relaxation factor is given by the following approximation for large values of N. ωoptimum ≈ 2 − 2π N [9] Regardless of how the relaxation factor is found, the iteration equation for SOR can be found by combining equations [8] and [6]. [ ] φij( n) = φij( n −1) + ω φij( n,GS ) − φij( n −1) = ωφij( n,GS ) + (1 − ω)φij( n −1) =ω bij − AijS φij( n−)1 − AijW φi(−n1) j − AijP AijE φi(+n1−1j ) − AijN φij( n+−11) + (1 − ω)φij( n −1) [10] Example problem The operation of the various iteration processes can be illustrated by using a finite-difference approximation of Laplace’s equation in two dimensions. ∂ 2u ∂x 2 + ∂ 2u ∂y 2 =0 0 ≤ x ≤ L, 0 ≤ y ≤ H [11] We can use the second-order central-difference expression for the second derivative. ∂ 2T ∂x 2 = ij ui +1 j + ui −1 j − 2uij (∆x) 2 ∂ 2T + O[(∆x) 2 ] ∂y 2 = uij +1 + uij −1 − 2uij ij (∆y ) 2 + O[(∆y ) 2 ] [12] Substituting these finite-difference representations into Laplace’s equation and ignoring the truncation error, which is second order in each independent variable, gives ui +1 j + ui −1 j − 2uij (∆x) 2 + uij +1 + uij −1 − 2uij (∆y ) 2 =0 [13] Multiplying through by (∆x)2 and defining β = ∆x/∆y gives the following finite-difference equation. ( ) β2uij −1 + ui −1 j + − 2 − 2β2 uij + ui +1 j + β2uij +1 = 0 [14] For this equation, the coefficients have the same values at all the ij points on the grid. AijN = AijS = β2 AijE = AijW = 1 AijP = −2 − 2β2 [15] AijP = −4 [16] For ∆x = ∆y, β = 1, and we have AijN = AijS = AijE = AijW = 1 The example problem uses a square grid with sides of length L, and ∆x/L = ∆y/L = 0.25 (Nx = Ny = 4). The boundary values have a different value at each node of the grid. The boundary values and the initial guesses for the nodes inside the grid, which are all zero, are shown below. j j j j j = = = = = 4 3 2 1 0 i = 0 4.00000 3.00000 2.00000 1.00000 0.00000 i = 1 6.00000 0.00000 0.00000 0.00000 2.00000 i = 2 8.00000 0.00000 0.00000 0.00000 4.00000 i = 3 10.00000 0.00000 0.00000 0.00000 6.00000 i = 4 12.00000 11.00000 10.00000 9.00000 8.00000 Files with iteration results Details of nine iteration steps, using Jacobi, Gauss-Seidel, and SOR are shown on the files jacobi.txt, gs.txt, and sor.txt. For each node, at each iteration, these files show the iteration equations, the values used in each equation, and the new iteration result for the node. The values for the current iteration are designated as u and the values for the old iteration are designated as uOld. For example, the first iteration at the first node for Gauss-Seidel iterations are shown as follows: Calculating u(1,1) = -[ AS(1,1) + AE(1,1) u(1,1) = [ -(1.00000)(2.00000) - (1.00000)(0.00000) - * u(1,0) + AW(1,1) * u(0,1) * uOld(2,1) + AN(1,1) * uOld(1,2) ] / AP(1,1) (1.00000)(1.00000) (1.00000)(0.00000) ] / (-4.00000) = 0.75000 For Gauss-Seidel and SOR, these equations are followed by a new table of nodal values. In the example above the following table follows the equations. This table shows the new value of u(1,1) just computed in the table. This value will be used in the subsequent computations of u(1,2) and u(2,1). j j j j j = = = = = 4 3 2 1 0 i = 0 4.00000 3.00000 2.00000 1.00000 0.00000 i = 1 6.00000 0.00000 0.00000 0.75000 2.00000 i = 2 8.00000 0.00000 0.00000 0.00000 4.00000 i = 3 10.00000 0.00000 0.00000 0.00000 6.00000 i = 4 12.00000 11.00000 10.00000 9.00000 8.00000 The detailed output tables are not shown for the Jacobi iteration, since that iteration process uses only the old values. For Jacobi iteration, the tables of iteration results are presented only at the end of an iteration step. At the end of each iteration step a table of residuals and relative changes are presented. The residuals, in general, are found from the following equation. Rij( n) = φij( n) + AijS φij( n−)1 + AijW φi( −n1) j + AijE φi( +n1−1j ) + AijN φij( n+−11) − bij AijP [17] For the problem solved here, bij = 0. The relative differences are found from the following equation. δij( n) = φij( n) − φij( n −1) φij( n) [18] The files with the iteration details should be available on the same web site where you found this document.