07/21/2003
To obtain Mk, for k = 0, 1, 2, ... , N + 1, we differtiate Sk(x)
(4) Sk'(x) = (3mk/6hk)(xk+1 - x)^2(-1) + (3mk+1/6hk)(x - xk)^2 + ((yk/hk) - (mkhk/6)(-1) + ((yk+1/hk) - (mk+1/6)hk)
Evaluating (4) Sk'(x) at x = xk
Sk'(xk) = (3mk/6hk)(xk+1 - xk)^2(-1) + (3mk+1/6hk)(xk - xk)^2 + ((yk/hk) - (mkhk/6)(-1) + ((yk+1/hk) - (mk+1/6)hk)
= (-1/2)((mk/hk)(hk)^2) - ((yk/hk) - (mk*hk)/6) + ((yk+1/hk) - (mk+1/6)hk)
= mkhk((-1/2) + (1/6)) - (1/6)(mk+1)(hk) + ((yk+1 - yk)/hk)
Assume ((yk+1 - yk)/hk) => dk
(5) Sk'(xk) = -((mkhk)/3) - ((mk+1hk)/6) + dkUsing (4) we obtain S'k-1(x)
- Substitute k by k-1S'k-1(x) = (3mk-1/6hk-1)(xk - x)^2(-1) + (3mk/6hk-1)(x - xk-1)^2 + ((yk-1/hk-1) - (mk-1hk-1/6)(-1) + ((yk/hk-1) - (mk/6)hk-1)
We evaluate S'k-1(x) at x = xk
S'k-1(xk) = (1/2)(mk-1/hk-1)(xk - xk)^2(-1) + (3mk/6hk-1)(xk - xk-1)^2 + ((yk-1/hk-1) - (mk-1/6)(hk-1)(-1) + ((yk/hk-1) - (mk/6)hk-1)
= (3mk/6hk-1)(hk-1)^2 - ((yk-1/hk-1) + (mk-1/6)(hk-1) + ((yk/hk-1) - (mk/6)(hk-1))
= (mk)(hk-1)((1/2) - (1/6)) + (mk-1hk-1/6) + (yk - yk-1)/hk-1
Assume (yk - yk-1)/hk-1 => dk-1
(6) S'k-1(xk) = (1/3)(mkhk-1) + (mk-1hk-1/6) + dk-1From property (4) we can make S'k-1(xk) = Sk'(xk)
(1/3)(mkhk-1) + (mk-1hk-1/6) + dk-1 = -((mkhk)/3) - ((mk+1hk)/6) + dk (mult by 6)
2(mkhk-1) + (mk-1hk-1) + 6dk-1 = -(2mkhk) - (mk+1hk) + 6dk
(mk-1hk-1) + 2mk(hk-1 + hk) + mk+1hk = 6(dk - dk-1)
Assume 6(dk - dk-1) => uk [for k = 1, 2, 3, ... , N-2]
(7) (mk-1hk-1) + 2mk(hk-1 + hk) + mk+1hk = ukWe have N - 1 equations and we have m0, m1, m2, ... , mN unknowns. These equations are needed to eliminate m0 and mN (extremes);
and for that we use different criterias.
From (7) we can tell that the system of equations to determine mk is :
This system is strictly diagonal dominant.
How to obtain m0 and mN? There are several criterias.
(1) CLAMPED SPLINE : The first derivatives at the ends are given -> S'(a) and S'(b) are given
((3/2)(h0 + 2h1))m1 + h1m2 = u1 - 3(d0 - S'(x0))
(hk-1mk-1) + 2(hk-1 + hk)mk + hkmk+1 = uk
(hN-2mN-2) + 2(hN-2+ (3/2)hN-1)mN-1 = uN-1 - 3(S'(xN) - dn-1)(2) NATURAL SPLINE : 2nd derivatives at the ends are zer0 -> S''(a) = 0 and S''(b) = 0
Normally used when trying to fit points.
(2.1) 2(h0 + h1)m1 + h1m2 = u1
(2.2) (hk-1mk-1) + 2(hk-1 + hk)mk + hkmk+1 = uk
(2.3) (hN-2mN-2) + 2(hN-2+ hN-1)mN-1 = uN-1(3) EXTRAPOLATED SPLINE : S0 and SN-1 are part of the same cubic polynomial. Also SN-2 and SN-1 are the same cubic polynomial.
(3h0 + 2h1 + (h0^2/h1)m1 + (h1 - (h0^2/h1))m2 = u1
(hk-1mk-1) + 2(hk-1 + hk)mk + hkmk+1 = uk
(hN-2 - (hN-1^2/hN-2)mN-2 + 2(hN-2 + 3hN-1 + (hN-1^2/hN-2))mN-1 = uN-1(4) PARABOLIC TERMINATED SPLINE : S0 and SN-1 will be quadratic polynomials instead of cubic
(3h0 + 2h1) + h1m2 = u1
(hk-1mk-1) + 2(hk-1 + hk)mk + hkmk+1 = uk
(hN-2mN-2) + 2(hN-2+ 3hN-1)mN-1 = uN-1(5) END-POINT CURVATURE ADJUSTED SPLINE : S''(a) and S''(b) are given
2(h0 + h1)m1 + h1m2 = u1 - h0S''(x0)
(hk-1mk-1) + 2(hk-1 + hk)mk + hkmk+1 = uk
(hN-2mN-2) + 2(hN-2+ hN-1)mN-1 = uN-1 - hN-1S''(xN)
07/22/2003
Solving a spline:
- Input a set of points (x0, y0) ... (XN, YN)
- Choose the criteria for restraints at the end points
- Solve the system of linear equations to determine m0 ... mN
- Find coefficients Sk,0 = yk , Sk,1 = dk - (hk(2mk - mk+1)/6) , Sk,2 = (mk/2) , Sk,3 = (mk+1 - mk) /6hkSk(x) = Sk,0 + Sk,1(x - xk) + Sk,2(x - xk)^2 + Sk,3(x - xk)^3
Example : Find the natural spline that approximates y = sin(x) in the interval [0, pi/6] with four points equally spaced.k x y
0 0 0
1 (pi/6) = .5236 sin(pi/6) = 0.5
2 (pi/3) = 1.0472 sin(pi/3) = .866
3 (pi/2) = 1.5708 sin(pi/2) = 1
natural spline => S''(a) = 0 and S''(b) = 0
m0 = S''(0) , m3 = mN = S''(pi/2) = 0 , m1 = ? , m2 = ?
h0, h1, h2 = pi/6 -> the change of xd0 = (y1 - y0) / h0 = (0.5 - 0) / 0.5236 = 0.9549
d1 = (y2 - y1) / h1 = (0.866 - 0.5) / 0.644
d2 = (y3 - y2) / h2 = (1 - 0.866) / 0.5236 = 0.2559u1 = 6(d1 - d0) = -1.5354
u2 = 6(d2 - d1) = -2.6586From the equations for natural spline criteria :
From (2.1) :
(a) 2(0.5236 + 0.5236)m1 + 0.5236m2 = -1.5354
We don't have equations coming from (2.2) since k goes from k = 2, .... , N-2 -> 3 -2 = 1
From (2.3) :
(b) 0.5236m1 + 2(0.5236 + 0.5236)m2 = -2.6586
Solving (a) and (b) :
Back Substitution : m2 = -2.2748 / 1.9635 = -1.1585
m1 = -(1.5354 - 0.5236(1.9635)) / 2.0944 = -0.4435So now we have that : m0 = 0, m1 = -0.4435, m2 = -1.1585, m3 = 0
We need to obtain the spline coefficients
Sk(x) = Sk,0 + Sk,1(x - xk) + Sk,2(x - xk)^2 + Sk,3(x - xk)^3
Sk,0 = yk
Sk,1 = dk - (hk(2mk - mk+1)) / 6
Sk,2 = mk / 2
Sk,3 = (mk+1 - mk) / 6hkk = 0 :
S0,0 = y0 = 0
S0,1 = 0.9549 - (0.5236(2(0) - 0.4435) / 6 ) = 0.9936
S0,2 = (0 / 2) = 0
S0,3 = (-0.4435 - 0) / 6(0.5236) = -0.1412S0(x) = 0.9936x - 0.1412x^3 -> the first spline
This spline have to pass through (0,0) and (0.5236, 0.5)
S0(0) = 0, S0(0.5236) = 0.5k = 1 :
S1,0 = y1 = 0.5
S1,1 = 0.6440 - (0.5326(2(-0.4435) - 1.1535) / 6 ) = 0.8775
S1,2 = (-0.4435 / 2) = -0.8218
S1,3 = (-1.1585 + 0.4435) / 6(0.5236) = -0.2276S1(x) = 0.5 + 0.8775(x - 0.5236) - 0.8218(x - 0.5236)^2 - 0.2276(x - 0.5236)^3
This spline have to pass through (0.5236, 0.5) and (1.0478, 0.866)
S1(0.5236) = 0.5, S1(1.0472) = 0.866Similarly, S2(x) = 0.866 + 0.4581(x - 1.0472) - 0.5293(x - 1.0472)^2 + 0.3688(x - 1.0472)^3
Numerical Differentiation : In some methods we need to know the derivative of a function
Limit of the Difference Quotient : Approximate the limit to 0 of the change of x by using a small value
f'(x) = lim(n->0) (f(x + h) - f(x)) / h
How small should h be to compute the limit numerically?
The value of h depends on the function. We can start with a large h and make it smaller and smaller until
- the division cannot be made
- the numerator is equal to zero
- you lose precision
07/23/2003
Numerical Differentiation:
Example :
f(x) = e^x, when x = 1
f'(x) ~= Dk = (e^(x+h) - e^x) / h
f'(1) = e^1 = 2.7183
Central Difference Formulas :
Center Formula of Order O(h^2) -> error
Using Taylor Expansion :
(1) f(x+h) = f(x) + f'(x)h + f''(x)h^2/2! + f'''(x)h^3 / 3!
where f'''(x)h^3 / 3! = error term
(2) f(x-h) = f(x) - f'(x)h + f''(x)h^2/2! - f'''(x)h^3 / 3!
where f'''(x)h^3 / 3! = error termSubtracting (2) from (1) :
f(x + h) - f(x - h) = 2f'(x)h + (f'''(c1) + f'''(c2))h^3 / 3!
f'(x) = ((f(x+h) -f(x-h)) / 2h) - (f'''(c1) + f'''(c2))h^3 / 3!(2h)
where (f'''(c1) + f'''(c2))h^3 / 3!(2h) = Truncation Error O(h^2)
f'(x) ~= (f(x+h) - f(x-h)) / 2hExample : f(x) = e^x, h = 0.0001
f'(1) ~= (e^(1 + 0.0001) - e^(1 - 0.0001)) / 2(0.0001) = 2.7183Example : f(x) = e^x, h = 0.001
f'(1) ~= (e^(1 + 0.001) - e^(1 - 0.001)) / 2(0.001) = 2.7183
Center Formula of Order O(h^4) :Do the Taylor Expansion :
(3) f(x+h) = f(x) + f'(x)h + f''(x)h^2/2! + f'''(x)h^3 / 3! + f''''(x)h^4 / 4! + f'''''(x)h^5 / 5!
(4) f(x-h) = f(x) - f'(x)h + f''(x)h^2/2! - f'''(x)h^3 / 3! + f''''(x)h^4 / 4! - f'''''(x)h^5 / 5!
(5) f(x + h) - f(x - h) = 2f'(x)h + 2(f'''(x))h^3 / 3! + 2(f'''''(c1) + f'''''(c2))h^5 / 5!Using step size 2h instead of h in (3):
f(h+2h) = f(x) + f'(x)2h + f''(x)4h^2/2! + f'''(x)8h^3 / 3! + f''''(x)16h^4 / 4! + f'''''(x)32h^5 / 5!
Doing the same to (5) :
f(x + 2h) - f(x - 2h) = 2f'(x)2h + 2(f'''(x))8h^3 / 3! + 2(f'''''(c1) + f'''''(c2))32h^5 / 5!
(6) f(x + 2h) - f(x - 2h) = 4f'(x)h + 16(f'''(x))h^3 / 3! + 64(f'''''(c1) + f'''''(c2))h^5 / 5!Multiply (5) by 8 and subtract (6) from it :
8f(x + h) - 8f(x - h) - f(x + 2h) - f(x - 2h) = (16f'(x)h + 16(f'''(x))h^3 / 3! + 16(f'''''(c1) + f'''''(c2))h^5 / 5!)
- (4f'(x)h + 16(f'''(x))h^3 / 3! + 64(f'''''(c1) + f'''''(c2))h^5 / 5!)
- f(x + 2h) + 8f(x + h) - 8f(x - h) + f(x - 2h) = 12f'x)h + kf'''''(c3)h^5
f'(x) = (- f(x + 2h) + 8f(x + h) - 8f(x - h) + f(x - 2h)) / 12h + k2(f'''''(c3)h^4) / 12
where k2(f'''''(c3)h^4) / 12 = O(h^4) error and
f'(x) = (- f(x + 2h) + 8f(x + h) - 8f(x - h) + f(x - 2h)) / 12h is formula of order O(h^4)Example : Let f(x) = sin(x), f'(x) = cos(x), x = pi/3, f'(pi/3) = 0.5, h = 0.001
Centered Formula of O(h^4):
f'(pi/3) = (-sin(pi/3 + 0.002) + 8sin(pi/3 + 0.001) - 8sin(pi/3 - 0.001) + sin(pi/3 - 0.002)) / 12(0.001) = 0.5
Centered Formula of O(h^2):
f'(pi/3) = (-sin(pi/3 + 0.001) - sin(pi/3)) / 2(0.001) = 0.499566904Numerical Integration : Obtaining area under the curve using numerical methods
07/24/2003
Numerical Integration :
Trapezoidal Rule :
A1 = ((f0 + f1) / 2)h , A2 = ((f1 + f2) / 2)h , A3 = ((f2 + f3) / 2)h , A4 = ((f3 + f4) / 2)h
A = A1 + A2 + A3 + A4 = h[f0/2 + f1 + f2 + f3 + f4/2]In general A = h[f0/2 + f1 + f2 + .... + fm-1 + fm/2] -> Trapezoidal Rule
h = (b - a) / mExample :
[-1, 1] e^(-x^2) dx , M = 4 , h = (1+1)/4 = 0.5
Example :
Exact Solution =Simpson Rule : Approximate f(x) using polynomials every 3 points.
Find quadratic polynomial that pass through (xo, f0) , (x1, f1) , (x2, f2)
Using LaGrange Polynomials : P2(x) = f0((x-x1)(x-x2))/((x0-x1)(x0-x2)) + f1((x-x0)(x-x2))/((x1-x0)(x1-x2)) + f2((x-x0)(x-x1))/((x2-x0)(x2-x1))The area below this P2(x) is :
Change x->x0 + ht , dx-> hdt , t = (x-x0) / h , t0 = 0 , t2 = (x2-x0) / h = 2
A0,2 =
+ f2
= f0
= f0= (h/3) [f0 + 4f1 + f2] =>
A = A0,2 + A2,4 + A4,6 = (h/3) [f0 + 4f1 + f2] + (h/3) [f2 + 4f3 + f4] + (h/3) [f4 + 4f5 + f6]
= (h/3)[f0 + 4f1 + 2f2 + 4f3 + 2f4 + 4f5 + f6]In general with 0 to 2m points A = (h/3)[f0 + 4f1 + 2f2 + .... + 2fm-1 + fm]
Example :
With 2m = 2 , h = ((pi/2) - 0)/3 = pi/4
A = ((pi/4)/3)[sin(0) + 4sin(pi/4) + sin(pi/2)] = 1.002279878
07/25/2003
Numerical Optimization : Find a maximum or minimum of an equation.
- Global maximum/minimum : max/min point in a function for -inf < x < inf
- Local maxmimum/minimum
- Unimodal function [a,b] : only one minimum or maximum in the function in the range [a,b]
The maximization of f(x) can be reduced to minimizing -f(x)Minimization of an equation :
- Local maximum value at x = P in the interval I : if f(x) <= f(P) for all x in I
- Local minimum value at x = P in the interval I : if f(x) >= f(P) for all x in I
A function is increasing in I if x1 < x2 and f(x1) < f(x2) for all x1,x2 in I
A function is increasing in I if x1 < x2 and f(x1) >f(x2) for all x1,x2 in I
Suppose that f(x) is continuous in I = [a,b] and is differentiable.
if f'(x) > 0 for all x = (a,b) then f(x) is increasing in I
if f'(x) < 0 for all x = (a,b) then f(x) is increasing in I
if there is a maximum or minimum value at x=P then f'(p) = 0if f''(P) < 0 then f(P) is a local maximum
if f''(P) >0 then f(P) is a local minimum
Search Methods :
The Golden Ratio : Assume that f(x) is unimodal in the interval [a,b], such that c = a + (1-r)(b-a) , d = a+r(b-a) , 0.5 < r < 1
Two cases :
(1) if f(c) < f(d) then make new search interval [a,d]
(2) if f(c) > f(d) then make new search interval [c,b]
Stop when the length of interval is less than some epsilon.
We can choose r in such a way that the new c or d can be a previous d or cif f(c) < f(d) then new b = d, we just need to compute c and f(c)
if f(c) > f(d) then new a = cWhat is the value of r that can allow us to do this?
1-r = r^2 , r^2 + r - 1 , r = (sqrt(5) - 1) / 2 => Golden RatioExample : f(x) = x^2 + 1 , r = (sqrt(5) - 1) / 2 = 0.61803
i = 1 :
c = a+(1-r)(b-a) = -0.8541
d = a+r(b-a) = -0.14589f(a) = f(-2) = 5
f(b) = f(1) = 2
f(c) = f(-0.8541) = 1.72949
f(d) = f(-0.14589) = 1.02128f(c) > f(d) so new interval = [-0.8541, 1]
i = 2 :
a = -0.8541 , b = 1 , c = previous d = -0.14589
d = a+r(b-a) = 0.297189f(a) = 1.72949
f(b) = 2
f(c) = 1.02128
f(d) = 1.08514f(c) < f(d) so new interval = [-0.8541, 0.297189]
i = 3 :
a = -0.8541 , b = 0.297189 , c = a+(1-r)(b-a) = 0.4164 , d = previous c = -0.14589
