Monday, June 28, 2004

 

-If we want to solve multiple systems of equations with the same A, then it is better to use LU factorization instead of Gaussian elimination:

            Assume we have systems of linear equations

                        Ax0  = b0  

                        Ax1  = b1 

                        Axn  = bn  

 

-If we use Gaussian Elimination n times

            Cost:  ln O(n 3) = O(ln n 3)           < --- Good for a single system

-If we use LU factorization

            Cost:  O(n 3) + ln O(n 2) = O(n 3) + O(ln n 2)        < --- Good for multiple systems

-LU will be faster

 

 

Iterative Methods for Systems of Linear Equations

 

-Iterative methods are better than the exact methods (Gaussian Elimination and LU factorization) when N is very large and 2A is a sparse matrix (sparse meaning that most of the elements in the matrix are zeros)

-For example, to solve partial differential equations numerically, very often you need to solve large (N > 100,000) and sparse matrices. Solving this kind of system using Gaussian Elimination or LU Factorization will take a long time

            N = 100,000        O((105)3) = O(1015)

 

Jacobi Iteration

-Assume the system:

            3x + y = 5

            x + 3y = 7

-We can write them as:

            x = (5 - y) / 3               y = (7 - x) / 3

-This suggests the following method:

            xk+1 = (5 - yk) / 3

            yk+1 = (7 - xk) / 3

-Assume the starting point:

            x0 = 0                                                   y0 = 0

            x1 = (5 - 0) / 3 = 1.667                         y1 = (7 – 0) / 3 = 2.33 

x2 = (5 – 2.33) / 3 = 0.889                   y2 = (7 – 1.667) / 3 = 1.777    

x3 = (5 – 1.777) / 3 = 1.074                 y3 = (7 – 0.889) / 3 = 2.037    

x4 = (5 – 2.037) / 3 = 0.987                 y4 = (7 – 1.074) / 3 = 1.475    

            xinfinity = 1                                             yinfinity = 2

 

Gauss-Seidel Iteration

-It is similar to Jacobi, but we do not wait until the end of the iteration to use the values of x and y etc…

-Gauss-Seidel uses the new approximations as soon as they are available

-Example:

            3x + y = 5

            x + 3y = 7

   We can write them as:

            x = (5 – y) / 3               y = (7 – x) / 3

   Now for the method:

            xk+1 = (5 - yk) / 3

            yk+1 = (7 - xk+1) / 3

   Assume the starting point:

            x0 = 0                                                   y0 = 0

            x1 = (5 - 0) / 3 = 1.667                         y1 = (7 – 1.667) / 3 = 1.778    

x2 = (5 – 1.778) / 3 = 1.074                 y2 = (7 – 1.074) / 3 = 1.975    

x3 = (5 – 1.975) / 3 = 1.008                 y3 = (7 – 1.008) / 3 = 1.99                  

            xinfinity = 1                                             yinfinity = 2

   Gauss-Seidel converges faster than Jacobi

 

-Jacobi and Gauss-Seidel may not converge in some cases

-For example, if we rearrange the previous equations:

            3x + y = 5                    x + 3y = 7

            x + 3y = 7        to         3x + y = 5

   Then if we create the iteration equations:

            x = 7 – 3y

            y = 5 – 3y

   Now for the method:

            xk+1 = 7 - 3yk

            yk+1 = 5 - 3xk+1

   Assume the starting point:

            x0 = 0                                                   y0 = 0

            x1 = 7 – 3(0) = 7                                  y1 = 5 – 3(7) = -16      

x2 = 7 – 3(-16) = 55                             y2 = 5 – 3(55) = -160               Does not converge!   

 

 

Tuesday, June 29, 2004

 

-How can we ensure convergence?

 

Strictly Diagonal Dominant Matrices

-A matrix A is strictly Diagonal Dominant if:

            | akk | > Σ | akj | where sigma is from j = 0 to N, where j ≠ k and where k = 1,2,…,N

-This means that the absolute value of each element in the diagonal has to be larger than the summation of the absolute value of all the other elements in the same row

 

Condition for convergence of Jacobi and Gauss-Seidel

-The Jacobi and Gauss-Seidel iteration will converge if A is strictly diagonal

-Example:

            3x + y = 5                    A = | 3  1 |      |3| > |1|         

            x + 3y = 7                           | 1  3 |      |3| > |1|              A is a strictly diagonal dominant matrix, so the Jacobi method converges

-Example:

            x + 3y = 7                    A = | 1  3 |      |1| !> |3|

            3x + y = 5                           | 3  1 |      |1| !> |3|             A is not a strictly diagonal dominant matrix, so the Jacobi method may not converge

 

Systems of non-linear equations

-Newton’s Method for systems of non-linear equations

 Assume:

            f1(x,y) = x² – y – 0.2 = 0

            f2(x,y) = y² – x – 0.3 = 0

-How do we solve this system?

-We extend Newton’s Method for equations of two variables

-Using Taylor’s expansion for functions of two variables around x0,y0 we have:

            f1(x,y) = f1(x0,y0) + (∂f1(x0,y0) / ∂x)(x - x0) + (∂f1(x0,y0) / ∂y)(y - y0) + …

            f2(x,y) = f2(x0,y0) + (∂f2(x0,y0) / ∂x)(x - x0) + (∂f2(x0,y0) / ∂y)(y - y0) + …

 

 

Wednesday, June 30, 2004

 

-Writing the approximation in matrix form we have:

            | f1(x,y) |    | f1(x0,y0) |    | (∂f1(x0,y0) / ∂x)(∂f1(x0,y0) / ∂y) | | x - x0 |

            | f2(x,y) | = | f2(x0,y0) | + | (∂f2(x0,y0) / ∂x)(∂f2(x0,y0) / ∂y) | | y - y0 | + …

-The matrix containing the partial derivatives is known as the Jacobian Matrix J

            | f1(x,y) |    | f1(x0,y0) |              | x - x0 |

| f2(x,y) | = | f2(x0,y0) | +  Jx0,y0 | y - y0 | + …<error>…

-To obtain the Newton Method we approximate f1(x,y) and f2(x,y) by using two terms and eliminating the error

-Then we set f1(x,y) = 0 and f2(x,y) = 0 and x = x1 and y = y1 for the next iteration

            | 0 |    | f1(x0,y0) |              | x1 - x0 |

| 0 | = | f2(x0,y0) | +  Jx0,y0 | y1 - y0 |

 

                       (  | f1(x0,y0) |              | x1 - x0 |  )

Jx0,y0-1 (- | f2(x0,y0) | =  Jx0,y0 | y1 - y0 |  )

 

                         | f1(x0,y0) |     | x1 - x0 | 

-Jx0,y0-1  | f2(x0,y0) | =  | y1 - y0 | 

 

            | x0 |              | f1(x0,y0) |     | x1 | 

| y0 | -Jx0,y0-1  | f2(x0,y0) | =  | y1 |

 

-Newton Method for systems of non-linear equations

-Notice the similarity with Newton Method for only one variable

x1 = x0 – (f(x0) / f'(x0))

 

| x1 |    | x0 |              | f1(x0,y0) |     

| y1 | = | y0 | -Jx0,y0-1  | f2(x0,y0) | 

 

-Example:

            f1(x,y) = x² – y – 0.2 = 0

            f2(x,y) = y² – x – 0.3 = 0

 

                  | (∂f1(x0,y0) / ∂x)(∂f1(x0,y0) / ∂y) |

            J = | (∂f2(x0,y0) / ∂x)(∂f2(x0,y0) / ∂y) |

 

                  | 2x   -1 |                        | 2y   1 |

J = | -1   2y |                 J-1 = | 1   2x |  * (1 / (2x * 2y -1))  ί Where J-1 is the cofactor matrix divided by the determinant

 

-You can also Gaussian Elimination to obtain the inverse of a matrix.

                   | 3  6  1 |

            A = | 7  8  2 |

                   | 4  5  1 |

 

| 1  0  0 | 3  6  1 |

| 0  1  0 | 7  8  2 |   ί Do Gaussian Elimination on both matrices and make A an identity matrix to get the inverse of A, A-1

| 0  0  1 | 4  5  1 | 

 

-Example:

 i             x                     y                                  f1                                             f2                                       J-1                    

0          x0 = 0            y0 = 0                   02 - 0 - 0.2 = -0.2                  02 - 0 - 0.3 = -0.3                 | 2(0)0     1   |

                                                                                                                                                             J-1 = |  1      2(0)0 | (1/(4(0)(0) – 1))

            | x1 |    | 0 |        | 0  1 | | -0.2 |    | -0.3 |

            | y1 | = | 0 | - -1 | 1  0 | | -0.3 | = | -0.2 |

 

1       x1 = -0.3        y1 = -0.2            (x1)2 – (y1) - 0.2 = 0.09              (y1)2 – (x1) - 0.3 = 0.04          | 2(-0.2)     1   |

                                                                                                                                                             J-1 = |  1      2(-0.3) | (1/(4(-0.3)(-0.2) – 1))

            | x2 |    | -0.3 |     | 0.526     -1.3158  | | 0.09 |    | -0.2947     |

            | y2 | = | -0.2 | -  | -1.3158  0.78951 | | 0.04 | = | -0.113158 |

 

2     x2 = 0.2947   y2 = -0.113158   (x2)2 – (y2) - 0.2 = 0.00000609   (y2)2 – (x2) - 0.3 = 0.0075

 

Interpolation and Polynomial Approximation

-We want to approximate functions using polynomials

-Taylor Approximation

            A polynomial PN(x) can be used to approximate a continuous function f(x)

                        f(x) ≈ PN(x) = ∑ (fk(x0) / k!) * (x - x0)k

                        f(x) = PN(x) + EN(x)

                        EN(x) = (fN+1(c) / (N + 1)!) * (c - y0)           for some c between x0 < c < x

-Examples:

            sin(x) = x – (x3 / 3!) + (x5 / 5!) + (x7 / 7!) + …

            ex = 1 + x + (x2 / 2!) + (x3 / 3!) + …

 

Methods to evaluate a polynomial

-Assume:

            f(x) = 2x4 + 4x3 + 3x2 + 2x + 1

                      2*x*x*x*x + 4*x*x*x + 3*x*x + 2*x*x +

          Number of multiplications = 10

          Number of sums = 4

-This is too expensive, we can factor the x using “Horner’s Method” 

            f(x) = (((2x + 4)x + 3)x + 2)x + 1

                      Number of multiplications = 4

                      Number of sums = 4

 

 

Thursday, July 1, 2004

 

Polynomial Approximation using N+1 points

-Suppose that a function y = f(x) is known to have N+1 points (x0, y0), (x1, y1), … , (xn, yn)

-We want to build a polynomial P(x) of degree N that will pass through the N+1 points

-Once we have the polynomial, we can use it to approximate f(x)

-If x0 < x < xN the approximation of f(x) using p(x) is called interpolation

-If x < x0 or xN < x the approximation of f(x) using p(x) is called extrapolation

 

 

-We are going to see two methods to build the polynomial from the N+1 points

 

Lagrange Approximation

-Assume two points (x0, y0), (x1, y1). The polynomial that passes through these points is a straight line.

 

 

#1 -> m = (y - y0) / (x - x0)  and  #2 -> m = (y1 - y0) / (x1 - x0)

 

Combining #1 and #2

            (y - y0) / (x - x0) = (y1 - y0) / (x1 - x0)

            (y - y0) = ((y1 - y0) / (x1 - x0)) * (x - x0)

            y = y0 + ((y1 - y0) / (x1 - x0)) * (x - x0)

 

            at x = x0                                                           at x = x1

                y = y0 + ((y1 - y0) / (x1 - x0)) * (x0 - x0)             y = y0 + ((y1 - y0) / (x1 - x0)) * (x1 - x0)

                            y = y0                                                               y = y1

 

-Lagrange uses a similar approach

-The function can also be written as:

            y = P1(x) = y0((x - x1) / (x0 - x1)) + y1((x - x0) / (x1 - x0))

            y0 = P1(x0) = y0((x - x1) / (x0 - x1)) + y1((x - x0) / (x1 - x0))

            y1 = P1(x1) = y0((x - x1) / (x0 - x1)) + y1((x - x0) / (x1 - x0))

-For P2(x)

            P2(x) = y0(((x - x1)(x - x2)) / ((x0 - x1)(x0 - x2))) + y1(((x - x0)(x - x2)) / ((x1 - x0)(x1 - x2))) + y2(((x - x0)(x - x1)) / ((x2 - x0)(x2 - x1)))

 

-Example:

            P2(x2) = y0(0) + y1(0) + y2(1) = y2

     For P3(x)

            P3(x) = y0(((x - x1)(x - x2)(x – x3)) / ((x0 - x1)(x0 - x2)(x0 - x3))) + y1(((x - x0)(x - x2)(x - x3)) / ((x1 - x0)(x1 - x2)(x1 - x3)))

+ y2(((x - x0)(x - x1)(x - x3)) / ((x2 - x0)(x2 - x1)(x2 - x3))) + y3(((x - x0)(x - x1)(x - x2)) / ((x3 - x0)(x3 - x1)(x3 - x2)))

-In general

            PN(x) = ∑ yK LN,K(x)                                                     where K = 0 and goes to N

            Where LN,K(x) = ( ∏ (x – xJ) ) / ( ∏ (xK – xJ) )              where J = 0 and K ≠ J and they go to N

 

-Example

            We want to approximate sin(x) in the interval 0, π/2 with a polynomial of degree 3 and 4 values

                        x             sin(x)

                        0              0

                        π/6           0.5

                        π/3           0.866

                        π/2           1

           

            P3(x) = 0(((x - π/6)(x - π/3)(x – π/2)) / ((0 - π/6)(0 - π/3)(0 - π/2))) + 0.5(((x - 0)(x - π/3)(x - π/2)) / ((π/6 - 0)(π/6 - π/3)(π/6 - π/2)))

+ 0.866(((x - 0)(x - π/6)(x - π/2)) / ((π/3 - 0)(π/3 - π/6)(π/3 - π/2))) + 1(((x - 0)(x - π/6)(x - π/3)) / ((π/2 - 0)(π/2 - π/6)(π/2 - π/3)))

 

            P3(x) = (0.5 / ((π/6)(-π/6)(-π/3))) * (x)(x - π/3)(x - π/2) + (0.866 / ((π/3)(-π/6)(-π/6))) * (x)(x - π/6)(x - π/2)

                        + (1 / ((π/2)(π/3)(π/6))) * (x)(x - π/6)(x - π/3)

 

                     = 1.74(x)(x – 1.047)(x – 1.5708) – 3.0164(x)(x – 0.5236)(x – 1.5708) + 1.1611(x)(x – 0.5236)(x – 1.047)

 

            Example for sin(π/4) = 0.7071

            π/4 = 0.7854

            P3(π/4) = 1.74(0.7854)(0.7854 – 1.047)(0.7854 – 1.5708) – 3.0164(0.7854)(0.7854 – 0.5236)(0.7854 – 1.5708)

   + 1.1611(0.7854)(0.7854 – 0.5236)(0.7854 – 1.047)

= 0.7054          (close to the true value of sin(π/4))

 

 

Friday, July 2, 2004

 

Newton Polynomials

-In Lagrange polynomials PN-1(x) and PN(x) are not related so we cannot use PN-1(x) to build PN(x)

-Newton Polynomials can be built incrementally, so we can build PN(x) using PN-1(x)

-Assume:

            P1(x) = a0 + a1(x - x0)

            P2(x) = a0 + a1(x - x0) + a2(x - x0)(x – x1)

P3(x) = a0 + a1(x - x0) + a2(x - x0)(x – x1) + a3(x - x0)(x – x1)(x - x2)

                     = P2(x) + a3(x - x0)(x – x1)(x - x2)

PN(x) = a0 + a1(x - x0) + a2(x - x0)(x – x1) + … + aN(x - x0)…(x - xN-1)

 

-How to compute a0, a1, … aN

-Assume we want to build P1(x) with points (x0, f(x0)) and (x1, f(x1))

-We want

            #1 -> P1(x0) = f(x0)                   #2 -> P1(x1) = f(x1)

                        #3 -> P1(x) = a0 + a1(x - x0)     Newton Polynomial N = 1

-Substitute #1 in #3

            P1(x0) = f(x0) = a0 + a1(x0 - x0) = a0

-Substitute #2 in #3

            P1(x1) = f(x1) = a0 + a1(x1 - x0)

                             a1 = (f(x1) - a0) / (x1 - x0)

                             a1 = (f(x1) - f(x0)) / (x1 - x0)

-Now if we want to obtain P2(x) with an additional point (x2, f(x2))

            P2(x) = a0 + a1(x - x0) + a2(x - x0)(x – x1)

P2(x2) = f(x2)

P2(x2) = f(x2) = a0 + a1(x2 - x0) + a2(x2 - x0)(x2 - x1)

 

            -Substituting the values of a0 and a1

            f(x2) = f(x0) + ((f(x1) – f(x0)) / (x1 – x0)) * (x2 - x0) + a2(x2 - x0)(x2 - x1)

            a2 = (1 / (x2 - x0)(x2 - x1)) * [ f(x2) – f(x0) - ((f(x1) – f(x0)) / (x1 – x0)) * (x2 - x0) ]

            a2 = (1 / (x2 - x1)) * [ (f(x2) – f(x0)) / (x2 - x0) - (f(x1) – f(x0)) / (x1 – x0) ]

            a2 = (1 / (x2 - x1)) * [ ((f(x2) – f(x0))(x1 – x0) - (f(x1) – f(x0))(x2 - x0)) / ((x2 - x0)(x1 – x0)) ]

            a2 = (1 / (x2 - x1)) * [ (x1f(x2) – x1f(x0) - x0f(x2) – x0f(x0) - x2f(x1) + x2f(x0) + x0f(x1) – x0f(x0) + x1f(x1) - x1f(x1)) / ((x2 - x0)(x1 – x0)) ]

a2 = (f(x2)(x1 – x0) - f(x1)(x1 – x0) - f(x1)(x2 – x1) +  f(x0)(x2 – x1)) / ((x2 - x1)(x2 - x0)(x1 – x0))

a2 = ((f(x2) - f(x1))(x1 – x0) – (f(x1) - f(x0))(x2 – x1)) / ((x2 - x1)(x2 - x0)(x1 – x0))

a2 = (1 / (x2 - x0)) * [ ((f(x2) - f(x1)(x1 – x0)) / ((x2 - x1)(x1 – x0)) - (f(x1) - f(x0)(x2 – x1)) / ((x2 - x1)(x1 – x0)) ]

a2 = (1 / (x2 - x0)) * [ ((f(x2) - f(x1)) / (x2 - x1) - (f(x1) - f(x0)) / (x1 – x0) ]

where:

   m1 = (f(x1) - f(x0)) / (x1 – x0)

   m2 = (f(x2) - f(x1)) / (x2 - x1)

 

-Observe that m1 and m2 are slopes for the lines between x1 ΰ x2 and x0 ΰ x1 they are called “divided differences”

-The divided differences of a function f(x) are defined as follows:

            f[xK] = f[xK]

            f[xK-1, xK] = (f[xK] – f[xK-1]) / (xK - xK-1)

            f[xK-2, xK-1, xK] = (f[xK-1, xK] – f[xK-2, xK-1]) / (xK - xK-2)

            …       

            f[xK-J, …, xK] = (f[xK-J+1, …, xK] - f[xK-J, …, xK-1] ) / (xK - xK-J)

            where aK = f[x0, x1, …, xK]

 

-Example

  -Build a polynomial of degree 1, 2, 3, to approximate f(x) = sin(x) in the interval of 0, π/2

            x0 = 0

            x1 = π/6

            x2 = π/3

            x3 = π/2

            P1(x) = a0 + a1(x - x0)

            P2(x) = a0 + a1(x - x0) + a2(x - x0)(x – x1)

                     = P1(x) + a2(x - x0)(x – x1)

P3(x) = a0 + a1(x - x0) + a2(x - x0)(x – x1) + a3(x - x0)(x – x1)(x - x2)

         = P2(x) + a3(x - x0)(x – x1)(x - x2)

 

xK                        f[xK]                    f[xK-1, xK]                                           f[xK-2, xK-1, xK]                                         f[xK-3, xK-2, xK-1, xK]

0                      sin(0) = 0 = a0              

π/6 = 0.5236      sin(π/6) = 0.5          .5 / .5236 = 0.9549 = a1

π/3 = 1.047       sin(π/3) = 0.8659     (.8659 – .5) / (1.047 – .5236) = 0.699    (.699 - .9549) / 1.047 = -.244 = a2

π/2 = 1.5078      sin(π/2) = 1             (1 - .8659) / (1.5078 – 1.047) = 0.256    (.256 - .699) / (1.5078 - .5236) = -.4230     (-.4230 +.244) / 1.5078 = -0.1137 = a3