6/28/04

Factor A in LU = O(n^3)
Gaussian Elimination = O(n^3)
Back substitution = O(n^2)
Forward substitution = O(n^2)

Assume we have m systems of linear equations
Ax[0] = b[0]
Ax[1] = b[1]
...
Ax[m] = b[m]

Using Gaussian elimination costs O(mn^3)
Using LU factorization costs O(n^3) + mO(n^2)
If m = 1, use Gaussian elimination
If m > 1, use LU factorization

------------------------------------------------------------------------------

Iterative methods for systems of linear equations

Better than exact methods when
* n is extremely large
* A is a sparse matrix

Good example of use is solving partial differential equations numerically
(for which solving via Gaussian elimination or LU would take a long time)

Google uses an iterative approach due to large sparse relation matrices

Jacobi iteration
----------------

Calculates increasingly more accurate values for 2 variables.
Example:
3x + y = 5 and x + 3y = 7
Rewrite as x = (5-y)/3 and y = (7-x)/3
x[k+1] = (5-y[k])/3
y[k+1] = (7-x[k])/3
Start at x[0] = 0 and y[0] = 0
x[1] = (5-0)/3 ≈ 1.667
y[1] = (7-0)/3 ≈ 2.333
x[2] = .889
y[2] = 1.778
x[3] = 1.074
y[3] = 2.037
x[4] = .987
y[4] = 1.975
x[∞] = 1
y[∞] = 2
3(1) + 2 = 5
1 + 3(2) = 7

Gauss-Seidel iteration
----------------------

Similar to Jacobi iteration, but doesn't include waiting until the end of iteration to use the values of x and y
(it instead uses the new approximations of x and y for computation as soon as they are available), causing it to converge faster.

Example:
3x + y = 5 and x + 3y = 7
x[k+1] = (5-y[k])/3
y[k+1] = (7-x[k+1])/3, as x[k+1] is already available for use
x[0] = 0
y[0] = 0
x[1] = (5-0)/3 = 1.667
y[1] = (7-1.667)/3 = 1.778
x[2] = 1.074
y[2] = 1.975

Jacobi and Gauss-Seidel

May not converge in some cases

If the equations had been rearranged as
   3x + y = 5
   x + 3y = 7
then we would get the iteration equations
   x = 7 - 3y
   y = 5 - 3x
which would cause divergence

6/29/04

How can we ensure convergence?

Strictly diagonal dominant matrix
|a[k][k]| > sum[j=1;j!=k;j->n] |a[k][j]|, where k goes from 1 to n
The absolute value of each element in the diagonal has to be larger than the summations of the absolute values
of all the other elements in the same row

Both Jacobi and Gauss-Seidel will necessarily converge if only strictly diagonal dominant matrices are used

Systems of non-linear equations - we extend Newton's method for equations in 2 variables

6/30/04

Writing approximation in matrix form, we have:

[x[1] = [x[0] - JJ^(-1)[f[1](x[0],y[0])
y[1]] y[0]] f[2](x[0],y[0])]

[f[1](x,y) = [f[1](x[0],y[0]) + [df[1](x[0],y[0])/dx df[1](x[1],y[1])/dy
f[2](x,y)] f[2](x[0],y[0])] df[2](x[0],y[0])/dx df[2](x[1],y[1])/dy]

The last matrix is called the Jacobian, represented by double-J
To obtain the Newton method, we approximate f[1](x,y) and f[2](x,y) by using two terms and eliminating the error.
Then we set f[1](x,y) = 0 and f[2](x,y) = 0 for the next iteration.

Example
-------
f[1](x,y) = x^2 - y - .2 = 0
f[2](x,y) = y^2 - x - .3 = 0
df[1](x,y)/dx = 2x
df[1](x,y)/dy = -1
df[2](x,y)/dx = -1
df[2](x,y)/dy = 2y
JJ = [2x -1
      -1 2y]
JJ^(-1) = [2y 1    (1/(2x2y - 1))
           1 2x]

To compute inverse, use determinant and multiply by (-1)^(ij), where i is the row number
and j is the column number (or put side-by-side with identity matrix, performing the same functions on both)

i      x      y      f[1]      f[2]      JJ^(-1)
---------------------------------------------------
0      0      0      -.2       -.3       [0 1    (1/(4(0)-1))
                                          1 0]

1 -.3 -.2 .09 .04 [.526 -1.3158
-1.3158 .7895]

2 .2947 -.113158 .000006090 .0075

We want f[1] and f[2] to approach 0 as i approaches infinity

------------------------------------------------------------------------------

Taylor approximation -> a polynomial P[N] approximates a continuous function x such that
f(x) ≈ P[N] = sum[k=0;N] (f(x[0])/k!)(x-x[0])^k

Examples:
---------

sin(x) = x - x^3/3! + x^5/5! + ....
e^x = 1 + x + x^2/2! + x^3/3! + ....
Assume ?(x) = 2x^4 + 4x^3 + 3x^2 + 2x + 1
Straight-forward: 10 multiplications and 4 additions
This is too expensive, but we can factor f(x) using "Horner's method"
f(x) = (((2x+4)x+3)x+2)x+1
Horner's method: 4 multiplications and 4 additions

7/1/04

Polynomial approximation using n+1 points

Suppose that a function y = f(x) is known to have n+1 points (x[0],y[0]), (x[1],y[1]), ... (x[n],y[n]) which fit it

We want to build a polynomial function p(x) of degree n that will pass through the n+1 points

Once we have p(x), we can use it to approximate f(x)
* If x[0] < x < x[n], approximating f(x) at a given point via p(x) is referred to as interpolation
* If x < x[0] or x > x[n], the approximation is called extrapolation

Lagrange Approximation

Developed long before computers (used to make calculation easier/possible)

Given two points (x[0],y[0]) and (x[1],y[1]), the polynomial that passes through these two points is a straight line
(m = (y-y[0])/(x-x[0]), resulting in y = y[0] + ((y[1]-y[0])/(x[1]-x[0]))(x-x[0]))
     * The Lagrange Approximation is not dissimilar to this line-drawing method
     * y = p[1](x) = y[0]((x-x[0])/(x[0]-x[1]))
     * The aim is to have p[1](x[0]) = y[0], p(x[1]) = y[1] as with the lines

Degree 2 approximation by curve
p[2](x) = y[0]((x-x[1])(x-x[2]))/((x[0]-x[1])(x[0]-x[2]))
+ y[1]((x-x[0])(x-x[2]))/((x[1]-x[0])(x[1]-x[2]))
+ y[2]((x-x[0])(x-x[1]))/((x[2]-x[0])(x[2]-x[1]))

Degree 3 approximation by polynomial
p[3](x) = y[0]((x-x[1])(x-x[2])(x-x[3]))/((x[0]-x[1])(x[0]-x[2])(x[0]-x[3]))
        + y[1]((x-x[0])(x-x[2])(x-x[3]))/((x[1]-x[0])(x[1]-x[2])(x[1]-x[3]))
        + y[2]((x-x[0])(x-x[1])(x-x[3]))/((x[2]-x[0])(x[2]-x[1])(x[2]-x[3]))
        + y[3]((x-x[0])(x-x[1])(x-x[2]))/((x[3]-x[0])(x[3]-x[1])(x[3]-x[2]))

General form
p[n](x) = ∑[k=0;n] y[k]*(π[j=0;j∫k;n](x-x[j]))/(π[j=0;j∫k;n](x[k]-x[j]))

Example
-------

Approximating sin(x) on the interval [0,π/2] with degree 3 polynomial
x sin(x)
-----------
0 0
π/6 .5
π/3 .866
π/2 1
p[3](x) = 0((x-π/6)(x-π/3)(x-π/2))/((0-π/6)(0-π/3)(0-π/2))
+ .5((x-0)(x-π/3)(x-π/2))/((π/6-0)(π/6-π/3)(π/6-π/2))
+ .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))
p[3](x) = 1.74(x)(x-1.047)(x-1.5708)
- 3.0164(x)(x-.5236)(x-1.5708)
- 1.1611(x)(x-.5236)(x-1.047)

sin(π/4) = .7071
p[3](π/4) = .7054

7/2/04

Exam

Located at approximately the mid-point of the course

Will be during the evening

Can bring half a page of notes to the exam, will turn in at end

Newton Method for Polynomial Approximation

P[n](x) can be built from P[n-1](x)

P[1](x) = a[0] + a[1](x-x[0])

P[2](x) = a[0] + a[1](x-x[0]) + a[2](x-x[0])(x-x[1])

P[n](x) = P[n-1](x) + a[n](π[i=0;n](x-x[i]))

P[1](x) is built from points (x[0],f(x[0])), (x[1],f(x[1]))

a[0] = f(x[0])

a[1] = (f(x[1]) - f(x[0]))/(x[1]-x[0])

a[2] = (((f(x[2])-f(x[1]))/(x[2]-x[1]))

((f(x[1])-f(x[0]))/(x[1]-x[0]))) / (x[2]-x[0])

Each coefficient is built up from "divided differences", f[x[n]]
f[x[k]] = f(x[k])
f[x[k-1],x[k]] = (f[x[k]] - f[x[k-1]])/(x[k] - x[k-1])
f[x[k-2],x[k-1],x[k]] = (f[x[k-1],x[k]] - f[x[k-2],x[k-1]])/(x[k] - x[k-2])

Example
-------

Build polynomials of degrees 1, 2, and 3 to approximate f(x) = sin(x) over the interval [0,π/2]

x[k]   f[x[k]]       f[x[k-1],x[k]]   f[x[k-2],x[k-1]]
----------------------------------------------------------
0       sin(0) = 0
π /6   sin(π/6) = .5         .9549
π /3   sin(π/3) = .8699   .699          -.244
π /2   sin(π/2) = 1          .836          -.423