By Patrick Piemonte
Index:
Announcements
Newthon's method for solving systems of Non-linear Equations
Example:
ƒ1(x, y) = x2 - y - .2 = 0
ƒ2(x, y) = y2 - x - .3 = 0
ƒ1(x, y) = ƒ1(x0, y0) + d/dx ƒ1(x0, y0)(x - x0) + d/dy ƒ1(x0, y0)(y - y0)
ƒ2(x, y) = ƒ2(x0, y0) + d/dx ƒ2(x0, y0)(x - x0) + d/dy ƒ2(x0, y0)(y - y0)
Expressing these two functions in matrix form:
|
= |
|
|
+ |
|
Jacobain Matrix (J):
| | d/dx ƒ1(x0, y0) d/dy ƒ1(x0, y0) | |
| | d/dx ƒ2(x0, y0) d/dy ƒ2(x0, y0) | |
Iteration Equation:
|
= |
|
+ |
|
-J |
|
= |
|
|
= | -J-1 |
|
|
= |
|
- | J-1 |
|
See similarity with Newthon-Rapson of one variable.
x1 = x0 - (ƒ(x0)/ƒ'(x0))
Example:
ƒ1(x, y) = x2 - y - .2 = 0
ƒ2(x, y) = y2 - x - .3 = 0
| Jx0, y0 | = |
|
| J-1 | = |
|
1/(4xy - 1) |
Start with x0 = 0, y0 = 0
x0 = 0
y0 = 0
ƒ1(x0, y0) = 0 - 0 - .2 = -.2
ƒ2(x0, y0) = 0 - 0 - .3 = -.3
| J-1x0, y0 | = |
|
1/(4(0)(0) - 1) |
=
| | 0 | -1 | |
| | -1 | 0 | |
|
= |
|
- |
|
|
=
| | -.2| |
| | -.3| |
x1 = -.3, y1 = -.2
ƒ1= (-.3)2 - (-.2) - .2 = .09
ƒ2= (-.2)2 - (-.3) - .3 = .04
| J-1x1, y1 | = |
|
1/(4(-.3)(-.2) - 1) |
= -.76
| J-1x1, y1 | = |
|
| | -.2947 | |
| | -.113158 | |
x2 = -.2947, y2 = -11.3158
ƒ1= (-.2947)2 - (-11.3158) - .2 = .000006090
ƒ2= (11.3158)2 - (.2947) - .3 = .6075
Why linear equations?
They are easy to solve. When you have a hammer everything is a nail.
What's the problem with that?
Not everything in the world is linear.
Interpolation and Polinomial Approximation
Taylor Expansion
The taylor expansion expresses a function with a polinomial of the form
ƒ(x) = (∑ from K=0 to Infinity) (ƒK(x0)/K!) (x - x0)K
If we want to approximate the function using a polinomial of degree N.
ƒ(x) = PN(x) + EN(x)
EN(x) = ƒN + F(c)/(N + 1)!
for some c such that c lies between x and x0
ex = ? approximate using a polynomial of degree 3. Assume x0= 0.
ƒ(x) = (∑ from K = 0 to 3) (ƒ(x)(x0)/K!)(x - x0)K
ƒ(x) = ex, ƒ(0) = e0= 1
ƒ'(x) = ex, ƒ'(0) = e0= 1
ƒ''(x) = ex, ƒ''(0) = e0= 1
ƒ'''(x) = ex, ƒ'''(0) = e0= 1
ƒ(x) = (ƒ0(x0)/0!)(x - x0)0 + (ƒ1(x0)/1!)(x - x0)1 + (ƒ'(x0)/2!)(x - x0)2 + (ƒ'''(x0)/3!)(x - x0)3
= (1/1)(1) + (1/1)(x)1+ (1/2)(x2) + (1/6)x3 ...
sin(x) = x - (x3/3!) + (x5/5!) - (x7/7!) + ...
Disadvantages of Taylor Expansion
The problem is that higher order derivatives must be known and maybe hard to compute.
Methods for evaluating polynomials:
Horner's Method (also called nested multiplication)
ƒ(x) = x4 + 2x3 + 3x2+ 2x + 1
(9 Multiplications, 4 sums)
ƒ(x) = (((x + 2) x + 3) x + 2)x + 1
Announcements
Suppose that a function ƒ(x) is known to have N + 1 points
(x0, y0), (x1, y1), ..., (xN, yN)
such that
a < x0 < x1 < x2 ... xN < b
and yK = ƒ(xK)
We want to construct a polynomial p(x) of degree N that passes through these N + 1 points.
The polynomial can be used for interpolation:
If x0 < x < xN
the approximation is PN(x)
Extrapolation:
If x < x0 or xN < x
PN(x) is the extrapolate value
Lagrange Approximation
Assume 2 points (x0, y0), (x1, y1)
Lets build the line that passes through these two points.
m = (y1 - y0)/(x1 - x0) = (y - y0)/(x - x0)
y - y0 = (x - x0)((y1 - y0)/(x1 - x0))
y = y0 + (y1 - y0)/(x1 - x0)
Lagrange used a different method to obtain this polynomial.
The line equation can be written as:
p1(x) = y = y0(x - x1)/(x0 - x1) + y1(x - x0)/(x1 - x0)
1 x = x0, 0 x = x0
0 x = x1, 1 x = x1
p1(x) = y0(x - x1)/(x0 - x1) + y1(x - x0)/(x1 - x0)
= y0
P2(x1) = y = y0(x1 - x1)/(x0 - x1) + y1(x1 - x0)/(x1 - x0)
= y1
The terms
L1, 0 = (x - x1)/(x0 - x1) and L1, 1 = (x - x1)/(x1 - x0)
are called "lagrange polynomials"
The same idea can be used for polynomials of degree 2, 3, ..., N
For P2(x) -> (x0, y0), (x1, y1), (x2, y2)
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))
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 - x3))
In general
PN(x) = (∑ of K=0 to N) yKLN, K(x)
LN, K(x) = ((∏ from J=0 to N) (x - xJ)/(∏ from J=0 to N) (xK - xJ))
Example:
We want to approximate sin(x) in the interval 0, pi/2 with a polynomial of degree 3.
| x | sin(x) |
| 0 | 0 |
| pi/6 | .5 |
| pi/3 | .866 |
| pi/2 | 1 |
P3(x) = 0((x - (pi/6))(x - (pi/3))(x - (pi/2))/(0 - (pi/6))(0 - (pi/3))(0 - (pi/2)))
+ .5((x - 0)(x - (pi/3))(x - (pi/2))/((pi/6) - 0)((pi/6) - (pi/3))((pi/6) - (pi/2)))
+ .866((x - 0)(x - (pi/6))(x - (pi/2))/((pi/3) - 0)((pi/3) - (pi/6))((pi/3) - (pi/2)))
+ 1((x - 0)(x - (pi/6))(x - (pi/3))/((pi/2) - 0)((pi/2) - (pi/6))((pi/2) - (pi/2)))
(.5/(pi/6)(-pi/6)(-pi/3))x(x - (pi/3))(x - (pi/2))
+ (.866/(pi/3)(pi/6)(-pi/6))x(x - (pi/6))(x - (pi/2))
+ (1/(pi/2)(pi/3)(pi/6))x(x - (pi/6))(x - (pi/3))
Example:
sin(pi/4) = .7071
P3(pi/4) = 1.74 (.7854)(.7854 - 1.047)(.7854 - 1.5708)
- 3.0164 (.7854)(.7854 - .5236)(.7854 - 1.5708)
+ 1.1611(.7854)(.7854)(.7854 - .5236)(.7854 - 1.047)
= .28078 + .4871 - .06246
= .7054
Error = |((.7071 - .7054)/1.7071)| = .0017/.7071
= .2% (.0024)
Polynomial Approximation
Newthon Polynomials
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)
PN(x) = a0 + a1(x - x0) + an(x - x0)...(x - xN-1)
How to compute a0, a1, ..., aN
Assume we want to build P1(x) that passes through the points (x0, ƒ(x0)) and (x1, ƒ(x1))
p0(x0) = ƒ(x0)
p1(x1) = ƒ(x1)
ƒ(x1) = P1(x0) = a0 + a1(x0- x0) = a0 = ƒ(x0)
P1(x1) = ƒ(x1) = a0 + a1(x1- x0)
ƒ(x0) = a0 + a1(x1- x0)
a1 = (ƒ(x1) - ƒ(x0)) / (x1- x0)
Example:
Build polynomial P1(x), P2(x), P3(x) that approximate ƒ(x) = sin(x) in the interval 0,pi/2
xk ƒ(xk) ƒ[xk-1,xk]
x0 0 sin(0) = 0 (.5-0)/(0.5236-0) = .9544 (a1)
x1 pi/6 = 0.5236 sin(pi/6) = .5
x2 pi/3 = 1.047 .8659 (0.8659-0.5)/(1.047-0.5236) = .699
x3 pi/2 = 1.5708 1 (1-0.869)/(1.5708-1.047) = 0.256
ƒ[xk-2, xk-1, xk] ƒ[xk-3, xk-2, xk-1, xk]
(.699-0.09549)/(1.047-0) = -.2404 (-.4230+0.244)/(1.5708-0) = -.1137
(.256-0.699)/(1.5708-0.5236) = -0.4230
P1(x) = 0 + 0.9549(x - 0) = 0.9549x
P2(x) = 0.9549x + (-.2444)(x - x0)(x - x0)
P3(x) = 0. .9549x + (-.2444)(x - x0)(x - x0) - 0.1137(x - x0)(x - x1)(x - x2)
sin(pi/4) = 0.7071
P1(.7854) = 0.9549(.7854) = .74998
P2(.7874) = .6997
P3(.7874) = .70582
Announcements
Padé Approximations
There are qutients of polynomial that are used to approximate functions:
RN, M(x) = PN(x)/QM(x) for a < x < b
Where
1) PN(x) = P0 + P1x + P2 x2 + PN xN
2) QM(x) = 1 + Q1x + Q2x2 + … + QMxM
The polynomials are constructed so ƒ(x) and RN, M(x) agree at x = 0 and also the derivatives agree at x = 0 (up to M ≠ N)
We want
ƒ(x) = PN, M(x) = PN(x)/QM(x)
3) ƒ(x) Qm(x) = Pn(x) = z(x)
Assume that ƒ(x) has the
4) ƒ(x) = a0 + a1x + a2x^2 +… Akx^k + …
Substitute ƒ(x) in 3 and also function 1 and 2
ƒ(x) QM(x) - PN(x) = C0x0 + C1 x + … CN+MxN+M = 0
Performing multiplication of polynomial and factoring
(a0 - p0) + x(a1 + q1(a0 - p1)) + x^2(a2 + q1a1 + q2(a0 - p2)) + … = c0x^0 + c1x + c2x^2 + …
To make the left side as small as possible the initial coefficeient can be 0.
a0 - p0 = 0
q1a0 + a1 - p1 = 0
q2a0 + q1a1 + a2 - p2 = 0
qM+N + qM-1qN+1….qN+M = 0
we get n+m+1 equations and we have n+m+1 variables to compute.
Example:
Find the Padé approximation R2, 2(x) for ƒ(x) = ln(1+x)/x start with the Maclurin expansion
ƒ(x) = 1 x/2 + x^2 /3 -…
Solution:
RN, M(x) = PN(x)/QM(x)
R2, 2(x) = P2(x)/Q2(x) = (P0 + P1x + P2x^2)/(1 + Q1x + Q2x^2)
We want ƒ(x) QN(x)
- PN(x) = 0
(1 x/2 + x^2/3 -…)(1 + Q1x + Q2x^2) - (P0+P1x + P2x^2) = 0
1-x/2 + x^2 /3 - x^3 /4 +…
Qx -Q/2 x^2 + Q/3x^3 - …
-P0 - P1x P2 x^2 = 0
1 - P0 = 0
-1/2 + Q1 - P1 = 0
1/3 - Q1/2 + Q2 - P2= 0
-1/4 + Q1/3 - Q2/2 = 0
1/5 - Q1/4 + Q2/3 = 0
P2 = 1/30
P1 = 7/10
P0 = 1
Q1 = 6/5
Q2 = 3/10
R2,2(x) = (1 + 7/10x + 1/30x^2)/(1 + 6/5x + 3/10x^2)
For x = 1 (ln(1+x))/x = .6931
R2, 2 = .6933
Given the approximation for ln(1+x)/x obtain an approximation for ln(1+x)
(ln(1+x))/x = R2,2(x)
ln(1+x) = x R2,2(x) = x.((1+7/10x+1/30x^2)/(1+6/5x+3/10x^2))
= (30x+21x^2+x^3)/(36+36x+9x^2)
Announcements
Curve Fitting
Given a set of points (x0, y0)(x1, y1)….(xN,yN) come up with the curve that best fit these points.
Polynomial Approximation Polynomial passes through all points.
Curve fitting - Curve does not necessarily passes through the points but it passes cloes to them.
Example:
An experimental procedure gives a set of points and you would like to obtain a function y = ƒ(s) that relates the points.
Least Squares Line
Assue that you have recorded the points (x0, y0)(x1, y1)….(xN,yN).
And we want to find the line y = Ax + B that best fits the recorded data. We have that ƒ(x) = y = Ax + B
We can express the approximation error as ek=ƒ(xk)-yk
There are several norms that can be used to tell us how far the curve ƒ(x) is form the data
Max error Einf (ƒ)=max 1<=k<=n {|ƒ(xk)-yk|}
Average error E1(ƒ)=1/N ∑(From K=1 to N){|ƒ(xk)-yk|}
Rott mean Square Error E2(ƒ)=1/N
∑(From K=1 to N)(ƒ(xk)-yk)^2
Given the points (x0, y0)(x1, y1)….(xN,yN).
The least square liney = ƒ(x) = Ax + B is the line that minimized the root mean square error E2(ƒ).
∂E(A,B)/∂A=1/N ∑(From K=1 to N) 2(Axk + B - Yk)Xk
1) =A ∑(From K=1 to N)xk^2 + B∑(From K=1 to N) xk-∑(From K=1 to N) xkYk= 0
∂E(A,B)/∂B=1/N ∑(From K=1 to N) 2(Axk+B-Yk)
2) = A ∑(From K=1 to N)Xk + BN-∑(From K=1 to N)Yk = 0
Solve 1 and 2 to obtain A and B for the line Y = Ax + B that best fits the data.
Example:
1 6 7 36 42
2 9 6 81 54
3 14 3 196 42
4 17 1 289 17
5 21 0 441 0
67 17 1043 185
Obtain the line Y = Ax + B that best fits the data
From (1) A 1043+67B-135=0
From (2) 67A+5B-17=0
ƒ(x) = - .5613 x + 10.11
Sometimes instead of fitting to a function ƒ(x) = Ax+B, we would like to fit the data points to the curve ƒ(x)
= AX^M where M is known and A is unknown.
We also want to minimize the square error
E(A)= A ∑(From K=1 to N) (Axk^M-Yk)^2
∂E(A,B)/∂A= ∑(From K=1 to N) 2(Axk^M-Yk)Xk^M=0
A= ∑(From K=1 to N) (Yk Xk^M/∑k=1 N Xk^(2M))
Exc: y=Ax^3
2 5.9 47.2 64
2.3 8.3 100.86 14.86
2.6 10.7 1888.06 308.92
2.9 13.7 334.13 599.32
3.2 17 557.66 1075
1227.44 2056.2
A = 1227.44/2056.28 = .5969
Y=.5969X^3