7/19
Piecewise Cubic Splines
splines used for auto design
Want to fit n+1 data points into n cubic polynomials
(x0,y0,),(x1,y1)...(xn,yn) where
a=x0 < x1 < x2 ... xn=b into N cubic polynomials
Sk(x) with coefficients of Sk,0, Sk,1, Sk,2, Sk,3 that satisfies the following properties:
I S(x) = Sk(x) = Sk,0 + Sk,1(x-xk) + Sk,2(x-xk)2 + Sk,3(x-xk)3 for [xk, xk,1] and k=0,1,...n-1
II Sk(xk) = yk for k=0,1,...N The polynomial passes through each point (xk,yk)
III Sk(xk+1) = Sk+1(xk+1) for k=0,1,...n-2 The polynomials are connected to eachother continuity in y
IV S'k(xk+1) = S'k+1(xk+1) for k,0,1,...n-2 continuity in the first derivative
V S"k(xk+1) = S"k+1(xk+1) for k,0,1...n-2 continuity in second derivative
For each polynomial Sk(x) there are four coefficients to be determined
Sk,0, Sk,1, Sk,2, Sk,3 (4 degrees of freedom)
There are N polynomials: That gives 4xN coefficients to be determined
From II S(xk)=yk for k=0,1,2...N and from I Sk(xk) = Sk,0 + Sk,1(x-xk) + Sk,2(x-xk)2 + Sk,3(x-xk)3 for [xk,xk+1] –> yk=sk,0
From II and with N+1 data points we have N-1 constraints
From III, IV, V each supplies N-1 constraints 3(N-1) –> 4N-2 constraints
These leaves two additional degrees of freedom that we call end-point constraints and in both S'(x) and S"(x) at x0 and xn and will be discussed later
Since Sk(x) is cubic, the second derivative is a line. We can obtain the equation using Lagrange approximation.
S"k(x)=S"(xk)((x-xk+1)/(xk-xk+1))+S"(xk+1)((x-xk)/(xk+1-xk))
So Sk’’(x)=S’’(xk) when x=xk
Sk’’(x)=S’’(xk+1) when x=xk+1
Use mk = S"(xk), mk+1 = S’(xk+1) and hk = xk+1 - xk
S"k= mk((xk+1 - x)/hk) + mk+1((x - xk)/hk)
7/20
Continuation
rules repeated
Better than Bezier curves due to continuity in the second derivative
The second derivative is a line (since the spline is cubic)
Integrating S’’[k](x) twice
S[k](x) =∫[x[k];x[k+1]]∫[x[k];x[k+1]] (S"[k](x))dx
S[k](x) =∫[x[k];x[k+1]]∫[x[k];x[k+1]]((m[k](x[k+1]-x)/(h[k]))+(m[k+1](x-x[k])/(h[k])))dx
S[k](x) =∫[x[k];x[k+1]](-(m[k]/2h[k])(x[k+1] - x)^2+(m[k+1]/2h[k])(x - x[k])^2+c[1])dx
S[k](x) =(m[k]/6h[k])(x[k+1]-x)^3+(m[k+1]/6h[k])(x-x[k])^3+c[1]x+c[2]
We can express c[1]x + c[2] as
c[1]x + c[2] = p[k](x[k+1]-x)+q[k](x-x[k])
(1)S[k](x) =(m[k]/6h[k])(x[k+1]-x)^3+(m[k+1]/6h[k])(x-x[k])^3+p[k](x[k+1]-x)+q[k](x-x[k])
Substituting x[k], x[k+1], S(x[k]) = y[k] and S(x[k+1])=y[k+1] this yields
S[k](x[k]) = y[k] = (m[k]/6h[k])(x[k+1] - x[k])^3 + (m[k+1]/6h[k])(x[k]-x[k])^2 + p[k](x[k+1] - x[k]) + q[k](x[k] - x[k])
S[k](x[k]) = y[k] = (m[k]/6)(h[k])^2 + P[k](h[k])
y[k] = (m[k]/6h[k])(h[k])^3) + p[k]h[k] = (m[k]/6)(h[k])^2) + p[k]h[k]
For x=x[k+1] y=y[k+1]=S[k](x[k+1])):
S[k](x[k+1]) = (m[k]/6h[k])(x[k+1] - x[k+1]) + (m[k+1]/6h[k])(x[k+1]-x[k])^3 + p[k](x[k+1] - x[k+1]) + q[k](x[k+1] - x[k])
y[k+1]=(m[k+1]/6h[k])h[k]^3 + q[k]h[k]=(m[k+1]/6)h[k]^2 + q[k]h[k]
Now obtain p[k] from this equation
(2) p[k] = (y[k]-(m[k]/6)h[k]^2)(1/h[k]) = y[k]/h[k] - (m[k]h[k]/6)
(3) q[k] = (y[k+1] - (m[k+1] / 6)h[k]^2)(1/h[k]) = y[k+1]/h[k] - (m[k+1]/6)h[k]
Substitute (2) and (3) into (1)
S[k](x) = (m[k]/6h[k])(x[k+1]-x)^3+(m[k+1]/6h[k])(x-x[k])^3+(y[k]/h[k] - (m[k]/6)h[k])(x[k+1]-x)+(y[k+1]/h[k] - (m[k+1]/6)h[k])(x-x[k])
We are almost there. Only the terms m[k] and m[k+1] are known. To obtain these values, we differentiate S[k](x).
(4) S'[k](x) = (3m[k]/6h[k])((x[k+1]-x)^2)(-1)+(3m[k+1]/6h[k])(x-x[k])^2+((y[k]/h[k])-(m[k])/6(h[k]))(-1)+((y[k+1]/h[k])-(m[k+1])/6(h[k]))
Evaluating S’[k](x) at x=x[k]
S'[k](x[k]) = (3m[k]/6h[k])((x[k+1]-xk)^2)(-1)+(3m[k+1]/6h[k])(xk-x[k])^2+((y[k]/h[k])-(m[k])/6(hk))(-1)+((y[k+1]/h[k])-(m[k+1])/6(hk))
= -(3m[k]/6h[k])h[k]^2 - (y[k]/h[k]) + (m[k]h[k]/6) + (y[k+1] / h[k]) - (m[k+1] / 6)(h[k])
= -(1/2)(m[k]h[k]) - y[k]/h[k] + (m[k]h[k]/6) + (y[k+1] / h[k]) - (m[k+1] / 6)(h[k])
= m[k]h[k](-(1/2)+(1/6))-((m[k+1]h[k] )/6) + ((y[k+1] - y[k]) / h[k])
Lets call d[k] = (y[k+1]-y[k])/h[k]
(5)S'[k](x[k]) = -m[k]h[k]/3 - (m[k+1]h[k])/6) + d[k]
Substituting k by k-1 in (4) to get the expression
S'[k-1](x) = (1/2)(m[k-1]/h[k-1])((x[k]-x)^2)(-1)+(3m[k]/6h[k-1])(x-x[k-1])^2+((y[k-1]/h[k-1])-(m[k-1])/6(h[k-1]))(-1)+((y[k]/h[k-1])-(m[k])/6(h[k-1]))
Evaluating S’[k-1](x) at x[k]
S'[k-1](x[k]) = (1/2)(m[k-1]/h[k-1])((x[k]-xk)^2)(-1)+(3/6)(m[k]/h[k-1])(xk-x[k-1])^2+((y[k-1]/h[k-1])-(m[k-1])/6(h[k-1]))(-1)+((y[k]/h[k-1])-(m[k])/6(h[k-1]))
stopped
7/21
Splines finale
Lets call d[k] = (y[k+1]-y[k])/h[k]
(6) S'[k-1](x[k]) = (1/3)(m[k]h[k-1])+(1/6)(m[k-1]h[k-1])+d[k-1]
From (4) S’[k](x[k])=S’[k-1](x[k]) and (5) and (6)
6(-(m[k]h[k]/3)-(m[k+1]h[k]/6)+d[k] = (1/3)m[k]h[k-1]+(m[k-1]h[k-1]/6))
-2m[k]h[k]-m[k+1]h[k]+6d[k]=2m[k]h[k-1] + m[k-1]h[k-1] + 6d[k-1]
6(d[k]-d[k-1]) = m[k-1]h[k-1] + 2m[k](h[k-1]+h[k]) + m[k+1]h[k]
Let U[k]= 6(d[k] - d[k-1]) for k=1,2,..N-1
(7) m[k-1]h[k-1] + 2m[k](h[k-1]+h[k]) + m[k+1]h[k] =U[k]
N-1 equations to solve m[0], m[1], ... m[N] that are N+1 unkowns. We need two more equations that are used to eliminate m[0],m[N]. We will use “end conditions” to determine m[0], m[N] and they will give us the different spline types. For now, assume for example that m[0]=s’’(x[0]) and m[N]=S’’(x[N]) and are given. Then from equation (7) we have
(8) k=1 m[0]h[0] + 2m[1](h[0]+h[1]) + m[2]h[1]=U[1]
(9) k=2...N-2 m[k-1]h[k-1] + 2m[k](h[k-1]+h[k]) + m[k+1]h[k] =U[k]
(10) k=N-1 m[N-2]h[N-2] + 2m[N-1](h[N-2]+h[N-1]) =U[N-1] - m[N]h[N-1]
So from (8), (9), (10) we have k=2,...N-2 equations with m[1],...m[N-1] unknowns.
This system is strictly diaganol dominant.
After the coefficients m[k] are computed, S[k] is computed as follows:
(11) S[k,0] = y[k]
(12) S[k,1] = d[k] - (h[k](2m[k]+m[n+1])/6
(13) S[k,2] = m[k]/2
(14) S[k,3] = (m[k+1]-m[k])/(6h[k])
These are coefficients of the polynomial k.
S(x)= S[k,0] + S[k,1](x-x[k]) + S[k,2](x-x[k])^2 + S[k,3](x-x[k])^3
(11), (12), (13), and (14) come from S(x) and their derivatives evaluated at S[k].
Eg.
S[k](x[k]) = S[k,0] = we have that S[k](x[k]) = y[k] so S[k,0] = y[k] that is (11)
To obtain (13) we can derive S[k](x)
S[k]’(x) = S[k,1] + S[k,2](x-x[k]) + 3S[k,3](x-x[k])^2
Again derivating
S[k]’’(x) = 2[Sk,2] + 6[Sk,3](x-x[k])
Now making x=x[k]
S[k]’’(x[k]) = 2S[k,2] + 6S[k,3](x[k]-x[k])
S[k,2] = (S[k]’’(x[k]))/2
From m[k]=S’’(x[k])
So we have S[k,2] = m[k]/2 that is (13)
We can obtain (12) and (14) in a similar way.
Other endpoint constraints
1) Clamped spline.
The first derivative at the ends are given. S’(x[0]) and S’(x[N]) are given. See equations of spline in book.
2) Natural spline.
The second derivatives at the end points are 0. S’’(x[0])=0=m[0] and S’’(x[N]) = 0 = m[N]
2(h[0]+h[1])m[1] + h[1]m[2]=U[1] K=1
h[k-1]m[k-1] + 2(h[k-1]+h[k])m[k]+h[k]m[k+1] = U[k] K=2,3,...N-2
h[N-2]m[N-2] + 2(h[N-2]+h[N-1])m[N-1]= U[N-1] K=N-1
3) Extrapolated spline
S[0] and S[1] are the same spline, S[N-2] and S[N-1] are the same spline
4)Parabolic terminated spline
S[0] and S[N-1] are quadratic instead of cubic
5)End-point curvature-adjusted spline
S’’(x[0]) and S’’(x[N]) are given
7/22
Example
-------
Find the natural spline that approximates y = sin(x) in the interval 0,π/2 with 4 points equally spaced. These points will be 0, π/6, π/3, and π/2.
k x y
------------------------------
0 0 0
1 π/6 = 0.5236 0.5
2 π/3 = 1.0472 0.8660
3 π/2 = 1.5708 1
Natural Spline
m[0] = 0, m[3] = 0
We need to compute m[1],m[2]
h[0] = h[1] = h[2] = π/6 = 0.5236
h[0] = x[1]-x[0]=π/6
d[0] = (y[1]-y[0])/h[0] = 0.9549
d[1] = (y[2]-y[1])/h[1] = 0.6990
d[2] = (y[3]-y[2])/h[2] = 0.2559
u[1] = 6(d[1]-d[0]) = -1.5354
u[2] = 6(d[2]-d[1]) = -2.6586
From 2)
a) 2(h[0]+h[1])m[1] + h[1]m[2]=U[1] K=1
b) h[k-1]m[k-1] + 2(h[k-1]+h[k])m[k]+h[k]m[k+1] = U[k] K=2,3,...N-2
c) h[N-2]m[N-2] + 2(h[N-2]+h[N-1])m[N-1]= U[N-1] K=N-1
Since N=3 the equation in b) goes from K=2,3...3-2=1
So no equations are needed from b)
From a)
2(π/6 + π/6)m[1] + (π/6)m[2]=-1.5354
d) 2.0944 m[1] + 0.5236 m[2] = -1.5354
From c)
h[1]m[1] + 2(h[1] + h[2])m[2]=U[2]
(π/6)m[1] + 2(π/6 + π/6)m[2] = -2.6386
e) 0.5326 m[1] + 2.0944 m[2] = -2.6386
So we have a system of two linear equations d) and e). Use gaussian elimination to solve it.
back substitution
m[2] = -2.2748/1.9635 = -1.1585
m[1] = (-1.5354-0.5236(-1.1585))/-0.4435 = -0.4435
We need to build the polynomials S[0],S[1],S[2]
S[k](x)= S[k,0] + S[k,1](x-x[k]) + S[k,2](x-x[k])^2 + S[k,3](x-x[k])^3
We need to obtain S[k,0] , S[k,1] , S[k,2] , S[k,3] for K=0,1,2 to build S0, S1, S2
We have the equations
S[k,0] = y[k] S[k,1] = d[k]-(h[k](2m[k] + m[k]+1))/6
S[k,2] = m[k]/2 S[k,3] = (m[k+1] - m[k])/6h[k]
S[0,0] = y[0] = 0
S[0,1] = (d[0]-h[0](2m[0]+m[1]))/6 = (0.9549 - 0.5236(2(0) + (-.4435))/6 = 0.9936
S[0,2] = m[0]/2 = 0/2 = 0
S[0,3] = (m[1]-m[0])/6h[0] = (-.4435)/6(0.5236) = -0.1412
S[0](x) = 0 + 0.9936 (x-0) + 0(x-0)^2 - 0.1412 (x-0)^3 = 0.9936x - 0.1412x^3
S[1,0] = y[1] = 0.5
S[1,1] = (d[1]-h[1](2m[1]+m[2]))/6 = (0.6990 - 0.5236(2(-.4435) + (-1.1585))/6 = 0.8775
S[1,2] = m[1]/2 = -0.4435/2 = -.2218
S[1,3] = (m[2]-m[1])/6h1 = (-1.1585+0.4435)/6(0.5236) = -0.2276
S[1](x) = 0.5 + 0.8775(x-0.5236) + 0.2218(x-0.5236)^2 - 0.2276(x-0.5236)^3
S[2,0] = y[2] = 0.8660
S[2,1] = (d[2]-h[2](2m[2]+m[3]))/6 = (-0.2559 - 0.5236(2(-1.1585) + 0)/6 = 0.4581
S[2,2] = m[2]/2 = -1.1585/2 = -0.5785
S[2,3] = (m[3]-m[2])/6h[2] = (0-(-1.1585))/6(0.5236) = 0.3688
S[2](x) = 0.8660 + 0.4581(x-1.0472) + 0.5743(x-1.0472)^2 + 0.3688(x-1.0472)^3
Numerical Differentiation
Estimate the derivative of a function numerically.
We have the following definition of a derivative as:
f’(x)=lim(h->0)((f(x+h)-f(x))/h)
7/23
Or if we want to obtain the derivative at multiple points:
D[k] = (f(x+h[k])-f(x))/h[k] for k=1,2...N
How small h[k] should be?
We can start with a large h[k] and make it smaller and smaller until we see that D[k] starts to lose precision.
f(x) = e^x get f’(x)
D=(e^x+h - e^x)/h x=1 (e^1+h - e^1)/h
k h[k] D[k]
0 .1 2.8588
1 .01 2.7319
2 .001 2.7196
3 .0001 2.7184
4 .00001 2.7183
5 .000001 2.7183
6 .0000001 2.7183
7 .00000001 2.7183
Stop when |D[k] - D[k-1]| < E
Unfortunatley this may require many computations.
Central Differences Formula
Another way that is more accurate and fast is to use a “Central Difference Formula”
Start with a Taylor expansion at f(x+h) and f(x-h).
(1) f(x+h) = f(x) + f’(x)h + (f’’(x)h^2)/2! + (f’’’(c1)h^3)/3!
(2) f(x-h) = f(x) - f’(x)h + (f’’(x)h^2)/2! - (f’’’(c2)h^3)/3!
^error
Subtracting (1) and (2)
f(x+h) - f(x-h) = f(x) + f’(x)h + (f’’(x)h^2)/2! + (f’’’(c1)h^3)/3! - f(x) + f’(x)h - (f’’(x)h^2)/2! + (f’’’(c2)h^3)/3!
f(x+h) - f(x-h) = 2f’(x)h + ((f’’’(c1)+(f’’’(c2))/3!)h^3
f’(x) = (f(x+h) - f(x-h))/2h - (((f’’’(c1)+(f’’’(c2))/3!)h^2(1/2))
^truncation error = O(h^2)
f’(x) ≈ (f(x+h) - f(x-h))/2h Central Difference Formula
Example: f(x) = e^x h = .0001
f'(1) = (e^(1+0.0001) - e^(1-0.0001))/2(0.0001) = 2.7183
Centered formula of order O(h^4)
This formula uses more points to estimate the derivative and therefore is more accurate. It uses four points instead of two.
(a) f(x+h) = f(x) + f’(x)h + (f’’(x)h^2)/2! + (f’’’(x)h^3)/3! + (f’’’’(x)h^4)/4! + (f’’’’’(c1)h^5)/5!
(b) f(x-h) = f(x) - f’(x)h + (f’’(x)h^2)/2! - (f’’’(x)h^3)/3! + (f’’’’(x)h^4)/4! - (f’’’’’(c2)h^5)/5!
^error
Subtracting a and b
f(x+h) - f(x-h) = f(x) + f’(x)h + (f’’(x)h^2)/2! + (f’’’(x)h^3)/3! + (f’’’’(x)h^4)/4! + (f’’’’’(c1)h^5)/5! - f(x) + f’(x)h - (f’’(x)h^2)/2! + (f’’’(x)h^3)/3! - (f’’’’(x)h^4)/4! + (f’’’’’(c2)h^5)/5!
(1) f(x+h) - f(x-h) = 2f’(x)h + 2((f’’’(x)h^3)/3!) + ((f’’’’’(c1)+(f’’’’’(c2))/5!)h^5
^error
Using step size 2h instead of h in (1) we have (2)
(2) f(x+2h) - f(x-2h) = 4f’(x)h + 2((f’’’(x)(2h)^3)/3!) + ((f’’’’’(c1)+(f’’’’’(c2))/5!)(2h)^5
f(x+2h) - f(x-2h) = 4f’(x)h + 16((f’’’(x)h^3)/3!) + 64((f’’’’’(c1)+(f’’’’’(c2))/5!)(h)^5
^error
Now multiply (1) by 8 and subtract (2) from it.
8f(x+h) - 8f(x-h) - f(x+2h) - f(x-2h) = 16f’(x)h + 16((f’’’(x)h^3)/3!) + 8((f’’’’’(c1)+(f’’’’’(c2))/5!)(h)^5 - 4f’(x)h + 16((f’’’(x)h^3)/3!) - 32((f’’’’’(c1)+(f’’’’’(c2))/5!)h^5
-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h) = 12f’(x)h - 24((f’’’’’(c1)+(f’’’’’(c2))/5!)(h)^5
^error
f’(x)=(-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h))/12h + (2/5!)(f’’’’’(c1)+(f’’’’’(c2))(h^4)
^error
error is O(h^4)
f'(x) ≈ (-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h))/12h Centered formula of order O(h^4)
Example: let f(x)=sin(x) h=0.001
Using the Centered formula of order O(h^4)
f’(pi/3) = (-sin(π/3 + 2(0.001)) + 8sin(π/3 + 0.001) - 8sin(π/3 - 0.001) + sin(π/3 - 2(0.001)))/12(0.001)
=(-0.8670 + 6.9322 - 6.9242 + 0.8654)/0.0120 = 0.5
Using the Centered formula of order O(h^2)
f’(x) = (sin(π/3 +0.001) - sin(π/3 - 0.001))/2(0.001) = 0.499999917
Using the Difference Quotient
f’(x) = (sin(π/3 + 0.001) - sin(π/3))/0.001 = 0.499566904