- Bookmarks
-- 6/28/2004
-- 6/29/2004
-- 6/30/2004
-- 7/1/2004
-- 7/2/2004

--* 6/28/2004

- Complexity of Gauss Elimination and LU Factorization
-- 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]
- If we use Gaussian Elimination n times
cost: m*O(n^3) = O(m*n^3)
- If we use LU factorization
cost: O(n^3)+m*O(n^2) = O(n^3)+O(m*n^2)
- LU will be faster if m is more than 1.
- Gauss will be faster if m is 1.


- Iterative Methods for Systems of Linear Equations
-- Iterative methods are better when:
1. n is very large.
2. A is a sparse matrix (most elements are 0).
- 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 will take a long time.
n = 100,000; O((10^5)^3) = O(10^15)

-- Jacobi Iteration
- Assume System
- 3x + y = 5
x +3y = 7
- We can write them as:
x = (5-y)/3 or
y = (7-x)/3
- This suggests:
Xk+1 = (5-Yk)/3 or
Yk+1 = (7-Xk)/3
- Let's Assume the starting point x0 = 0, and 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 -.889)/3 = 2.037
x4 = (5-2.037)/3 = 0.987 y4 = (7-1.074)/3 = 1.975
... ...
x(infinity) = 1 y(infinity) = 2

3(1) + (2) = 5 and
(1) + 3(2) = 7

-- Gauss-Siedel Iteration
- Similar to Jacobi, but we do not wait until end of iteration to use the values
of x,y,z...etc.
- Uses the new approximations as soon as they are available:
- Example:
3x + y = 5
x +3y = 7

x = (5 - y)/3 y = (7 - x)/3
- note: So far, same as Gauss-Siedel.

x0 = 0 y0 = 0
x1 = ( 5 - 0 )/3 = 1.667 y1 = ( 7 - 0 )/3 = 1.77
x2 = ( 5-1.77)/3 = 1.074 y2 = (7-1.074)/3 = 1.975
x3 = (5-1.975)/3 = 1.008 y3 = (7-1.008)/3 = 1.99
- note: Much faster than Jacobi.
- Gauss-Siedel and Jacobi may not converge in some situations.
- Example: Rearrange...
3x + y = 5 ----> x +3y = 7
x +3y + 7 ----> 3x + y = 5
- Now the iteration equations would be...
x = 7-3y
y = 5-3x
- Gauss-Siedel...
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.

--* 6/29/2004

- Testing for Convergence
-- Strictly Diagonal Dominant Matrices
- A matrix A is strictly diagonal matrix if abs(akk) > sum(j=1;n)[abs(akj)]; k=1,2,...n
( j!=k)
- This means that the absolute value of each element in the diagonal has to be larger
than the summations of the absolute value of all the other elements in the same row.
- Condition for the convergence of Jacobi and Gauss-Siedel...
- The Jacobi and Gauss-Siedel iteration will converge if A is strictly diagonal,
- Example:
3x + y = 5 ______\ [ 3 1 ] abs(3) > 1 (check)
x +3y = 7 / [ 1 3 ] abs(3) > 1 (check)

x(k+1) = (5-yk)/3 --> A is a diagonal dominant matrix so
y(k+1) = (7-xk)/3 Jacobi converges.

x +3y = 7 ______\ A = [ 1 3 ] abs(1) !> 3
3x + y = 5 / [ 3 1 ] abs(1) !> 3

x(k+1) = 7 - 3yk --> A is not a diagonal dominant matrix,
y(k+1) = 5 - 3yk so Jacobi will not converge.


- Systems of Non-Linear Equations
-- Newton's Method
- Assume:
f1(x,y) = (x^2)-y-0.2=0
f2(x,y) = (y^2)-x-0.3=0
- How to solve:
- We extend the Newton 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)+(df1(x0,y0)/dx)(x-x0)+(df1(x0,y0)/dy)(y-y0)+...
f2(x,y) = f2... f2... f2...

--* 6/30/2004

- Writing approximation in matrix form we have:
[f1(x,y)] = [f1(x0,y0)]+[(d/dx)f1(x0,y0)*(d/dy)f1(x0,y0)]*[x-x0]+...
[f2(x,y)] [f2(x0,y0)] [(d/dx)f2(x0,y0)*(d/dy)f2(x0,y0)] [x-x0] ...
\_________ Jacobian (J) ________/

[f1(x,y)] = [f1(x0,y0)]+J[x-x0] + ..error..
[f1(x,y)] [f2(x0,y0)] [x-x0] ..error..

- To obtain the Newton Method we approximate f1(x,y) and f2(x,y) by using the 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)] + J(x0,y0)[x1-x0]
[0] [f2(x0,y0)] [y1-y0]

(J(x0,y0)^-1){-[f1(x0,y0)] = J(x0,y0)[x1-x0]}
{ [f2(x0,y0)] [y1-y0]}

-(J(x0,y0)^-1)[f1(x0,y0)] = [x1-x0]
[f2(x0,y0)] [y1-y0]

[x0]-(J(x0,y0)^-1)[f1(x0,y0)] = [x1]
[y0] [f2(x0,y0)] [y1]

- Newton Method for systems of nonlinear equations
- Notice similarity with Newton Method for only one var...
x1 = x0 - (f(x0)/f`(x0))

[x1] = [x0] - (J(x0,y0)^-1)[f1(x0,y0)]
[y1] [y0] [f2(x0,y0)]

- Example
f1(x,y)=(x^2)-y-.2=0 __\ J=[2x -1]
f2(x,y)=(y^2)-x-.3=0 / [-1 2y]
- Eliminate row and column & compute determinant of remaining matrix times
(-1)^(i+j) ___\ (J^-1) = [2y 1]* ___1___
/ [1 2x] 2x*2y-1
- You can also do \ / \ one /
Gaussian Elimination cofactor over
to obtain inverse. matrix determinant

- Gaussian Inverse
[100][123] > Gauss Elimination in both matrices and make
[010][456] A an identity matrix to obtain inverse.
[001][789] > ie, perform same operations to both matrices.
I A

- continued example...
i x y f1((x^2)-y-.2) f2((y^2)-x-.3) (j^-1)
----------------------------------------------------------------------------
0 0 0 0-0-.2=-.2 0-0-.3=-.3 [2y=0 1](__1__) = [01]
[1 2x=0](4xy-1) [10]
1 -.3 -.2 .09 .04 [.526 -1.3]
[ -1.3 .79]
2 -.2947 -.1132 6.09*10^-7 .0075
... approaching 0

- 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) = sum(k=1 to n){f^(`k)(x0)/k!}(x-x0)^k

f(x) = Pn(x) + En(x)
En(x) = (((f`(n+1)(c))/(n+1)!)(c-y0) -> For some error
(f prime n+1 of c) between x0 <<<<x

- Examples:
sin(x) = x - (x^3)/3! + (x^5)/5! + (x^7)/7! ...
e^x = 1 + x + (x^2)/2! + (x^3)/3! ...
- Methods to evaluate a polynomial

Assume: f(x) = 2(x^4) + 4(x^3) + 3(x^2) + 2x + 1
- 2xxxx + 4xxx + 3xx + 2x + 1
yields...
-> 10 multiplications
-> 4 sums
which is expensive.

Factor using Horner's: f(x) = 1+x(2+x(3+x(4+2x)))
which yields...
-> 4 multiplications
-> 4 sums
which saves time.



--* 7/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.

-- 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)
2. m = (y1-y0)/(x1-x0)

- combine 1 & 2...
(y-y0) = (y1-y0)
------ -------
(x-x0) (x1-x0)


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

- at x = x0, y = y0
- at x = x1, y = y1

- Lagrange uses similar approach...
- The function can be written as
- degree > 1

y = P1(x) = y0(x-x1)/(x0-x1) + y1(x-x0)/(x1-x0)
P1(x0)= y0(x0-x1)/(x0-x1) + y1(x0-x0)/(x1-x0) = y0
L__ equals 1 L__ equals 0
P1)x1_= y0(x1-x1)/(x0-x1) + y1(x1-x0)/(x1-x0) = y1
L__ equals 0 L__ equals 1

- degree > 2

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

- Example:

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

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

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

-- In general:

Pn(x) = sum(k=0 to n) Yk * L(n,k)(x) --> (L n k of x)

prod(j=0 to n while j!=k) (x-xj)
L(n,k)(x) = -----------------------------------
prod(j=0 to n while j!=k) (xk - xj)
L__ x sub k,j...

- Example:

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

x | sin(x)
------+---------
0 | 0
pi/6 | .5
pi/3 | .866
pi/2 | 1


(x-pi/6)(x-pi/3)(x-pi/2) (x-0)(x-pi/3)(x-pi/2)
P3(x) = 0*------------------------ +.5------------------------------>>
(0-pi/6)(0-pi/3)(0-pi/2) (pi/6-0)(pi/6-pi/3)(pi/6-pi/2)

(x-0)(x-pi/6)(x-pi/2) (x-0)(x-pi/6)(x-pi/3)
>> .866------------------------------ +1------------------------------
(pi/3-0)(pi/3-pi/6)(pi/3-pi/2) (pi/2-0)(pi/2-pi/6)(pi/2-pi/3)


= 1.74x(x-1.047)(x-1.5708)-3.0164x(x-.5236)(x-1.5708)+1.1611x(x-.5236)(x-1.047)

--> so for example, sin(pi/4) = .7071; pi/4 = .7854....
P3(pi/4) = .7054

--* 7/2/2004


- Newton Polynomials
- In Lagrange Polynomials P(n-1)(x) and Pn(x) are not related so we cannot use P(n-1)(x) to build Pn(x).
- Newton Polynomials can be built incrementally, so we can build Pn(x) using P(n-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)
... and so on

SO,
P3(x) = P2(x) + a3(x-x0)(x-x1)(x-x2)
P2(x) = P1(x) + a2(x-x0)(x-x1)

- in general...
Pn(x) = a0+a1(x-x0)...+an(x-x0)*...*(x-x(n-1))
Pn(x) = P(n-1)(x) + an(x-x0)*...*(x-x(n-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 degreen n=1)
- substitute 1. in 3.
P1(x0) = f(x0) = a0 + a1(x0-x0)---> equals 0
f(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)
we want to obtain a2...
 

-- Divided differences

> f[xk] = f(xk)
f[Xk]-f[Xk-1]
f[Xk-1,Xk] = -------------
(Xk) - (Xk-1)

f[Xk-1,Xk]-f[Xk-2,Xk-1]
f[Xk-2,Xk-1,Xk] = -----------------------
(Xk) - (Xk-2)
:
:
in general...


> f[Xk-j+1,....Xk] - f[Xk-j ... Xk-1]
< f[Xk-j,...Xk] = -----------------------------------
> (Xk) - (Xk-j)
< where ak = f[X0,X1,... Xk]


- Example:
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)

i | Xk | f[Xk] | f[Xk-1,Xk] | f[Xk-2,Xk-1,Xk]
--+--------+----------------+-----------------------+----------------------------
0 | 0 | sin(0) = 0=a0 | |
--+--------+----------------+-----------------------+----------------------------
1 | pi/6= | sin(pi/6)=.5 | (.5-0)/(.5236-0)= |
| .5236 | | .9549 = a1 |
--+--------+----------------+-----------------------+----------------------------
2 | pi/3= | sin(pi/3)= | (.8659-.5) | (.699 - .9549)
| 1.047 | .8659 | ---------- = .699 | -------------- = -.244 = a2
| | | (1.047-.5236) | (1.047 - 0)
--+--------+----------------+-----------------------+----------------------------
3 | pi/2= | sin(pi/2) = 1 | (1-.8659) | (.256-.699)
| 1.5078 | | ----------- = .256                           | -------------- = -.4230
| |              (1.5078-1.047)                               | ( 1.5078-.5263)
|--------+----------------+-----------------------+-----------------------------
| f[Xk-3,Xk-2,Xk-1,Xk]
|-------------------------------------------------------------------------------
| (-.4230-(-.244))
| ---------------- = -.1137 = a3
| (1.5078-0)
|_______________________________________________________________________________