Transcript ENGRD 241 / CEE 241: Engineering Computation 1

Engineering Computation
Curve Fitting: Interpolation 1
Curve Fitting by
Interpolation
Engineering Computation
Curve Fitting: Interpolation 2
Curve Fitting:
1. We have discussed Least-Squares Regression where the
function is "best fit" to points but does not necessarily pass
through the points.
2. We now discuss Interpolation & Extrapolation
The function passes through all (or at least most) points.
interpolation
extrapolation
Engineering Computation
Curve Fitting: Interpolation 3
Curve Fitting by Interpolation:
C&C covers four approaches:
1. Polynomials (C&C 18.1-18.5 – “skim” only)
a) n + 1 equations & n + 1 unknowns
b) Lagrange Polynomials
c) Newton's Divided Difference (NDD) Polynomials
2. Splines (C&C 18.6 – assigned reading)
The first 3 approaches find the same polynomial. We will only
cover superficially, and concentrate on Splines.
Engineering Computation
Curve Fitting: Interpolation 4
Curve Fitting by Interpolation:
General Scheme
Given:
Set of points (xi,yi), not necessarily evenly spaced or
in ascending order.
Assume: x = independent variable;
y = dependent variable.
Find:
y = f(x) at some value of x not in the set (xi).
Method: Determine the function f(x) which passes through all
(or most) points.
Engineering Computation
Curve Fitting: Interpolation 5
n+1 Equations and n+1 Unknowns (C&C 18.3)
Given n+1 data points (xi,yi), find an nth order polynomial:
y = a0 + a1x + a2x2 + … + anxn
that will pass through all the points.
Too much work!
• Equations are notoriously ill-conditioned for large n.
• Equations are not diagonally dominant.
• Method is rarely used.
Engineering Computation
Curve Fitting: Interpolation 6
Lagrange Interpolating Polynomials (C&C 18.2)
Given n+1 data points (xi,yi), find the nth order polynomial:
y = pn(x) = a0 + a1x + a2x2 + … + anxn
that passes through all of the points. The Lagrangian polynomials
approach employs a set of nth order polynomials, Li(x), such that:
n
p n (x) 
y
i * L i (x)

i 0
where Li(x) satisfies the condition:

1 at x = x i
Li (x)  0 at x = x where j  i
j
Engineering Computation
Curve Fitting: Interpolation 7
Newton's Divided Difference (NDD) Polynomial (C&C 18.1)
Gives the same polynomial as the Lagrange method but is
computationally easier. General form for n+1 data points:
pn(x) = b0 + b1(x–x0) + b2(x–x0)(x–x1)
+ … + bn(x–x0)(x–x1)(x–x2)…(x–xn-1)
with b0, b1, ... , bn all unknown.
Note that ith term is zero at xj for j < i. Each term insures that
the polynomial correctly interpolates at one new point. The
algorithm is recursive and readily suited for spreadsheet or
other programmed calculation.
i
xi
f(xi)
First
Second
Third
0
x0
f(x0)
f[x1,x0]
f[x2,x1,x0]
f[x3,x2,x1,x0]
1
x1
f(x1)
f[x2,x1]
f[x3,x2,x1]
2
x2
f(x2)
f[x3,x2]
3
x3
f(x2)
Engineering Computation
Curve Fitting: Interpolation 9
Newton's Divided Differences (NDD)
versus
Lagrange Polynomials:
1. Both methods give the same results.
2. Comparison based on a count of the FLOPS:
Evaluate coefficients: Interpolate for one x:
3.
Lagrange:
(n+1) (n+1)  n2
(n+1) (n+1)  n2
NDD:
(n+1) (n+1)/2  n2/2
n
Easy to add a node with NDD.
Need to start over with Lagrange.
4. Both methods share a major problem: as the number of points
increases, so does the order of polynomial. This may cause
excessive "wiggles" or "waves" between points.
Engineering Computation
Curve Fitting: Interpolation 10
Splines (C&C 18.6)
Issue:
Need to overcome the "wiggle" or "wave" problem
Idea:
Use a piecewise polynomial approximation
Simplest idea:
Straight line on each segment.
4
The problem is that
3
g'(x) and g"(x) are discontinuous
2
1
0
0
1
2
3
4
5
Engineering Computation
Curve Fitting: Interpolation 11
Splines (cont.)
Most frequently used: Cubic Splines
==> Separate Cubic polynomial on each interval.
This is the analytical/numerical analog of a flexible edge
(spline) which is used by draftsmen. With this tool the first
and second derivatives (curvature) are continuous, and the
function appears "smooth".
Engineering Computation
Curve Fitting: Interpolation 12
Cubic Splines
Objective:
Define a 3rd-order polynomial for each interval:
fi(x) = aix3 + bix2 + cix + di
For n+1 data points, (x0,y0), (x1,y1), … , (xn,yn), there are
n intervals with
4 unknowns per interval (ai, bi, ci, and di).
==> Total 4n unknowns.
4
3
We need 4n equations to
compute all 4 unknowns for
each interval
2
1
0
0
1
2
3
4
5
n=3, 3 segments, n+1 points
Engineering Computation
Curve Fitting: Interpolation 13
Cubic Splines
How do we obtain the required 4n equations?
• functions must pass through fi(x) at knots (points):
yi-1 = ai(xi-1)3 + bi(xi-1)2 + ci(xi-1) + di
yi = ai(xi)3 + bi(xi)2 + ci(xi) + di
[2 equations/interval = 2n]
• 1st and 2nd derivatives must be equal at interior knots (xi,yi):
3ai(xi)2 + 2bixi + ci = 3ai+1(xi)2 + 2bi+1xi + ci+1
6aixi+ 2bi = 6ai+1xi+ 2bi+1
[ 2 equations/interior knot = 2n-2 ]
TOTAL: 2n + (2n - 2) = 4n - 2
Engineering Computation
Curve Fitting: Interpolation 14
Cubic Splines
We need an extra 2 conditions for Cubic Splines
Natural Splines
Setting the 2nd derivatives at exterior knots equal to zero
allows the function to "relax":
0 = 6a1x0 + 2b1
&
0 = 6anxn+1 + 2bn
[ 1 equations/exterior knot = 2 ]
TOTAL: 2n + (2n - 2) +2 = 4n
Engineering Computation
Curve Fitting: Interpolation 15
Cubic Splines
Alternatives for two extra conditions
(instead of setting 2nd derivative = 0)
1) Specify 1st derivatives at exterior knots.
f '(x0) = 3a1x02 + b1x0
f '(xn+1) = 3anxn+12 + bn xn+1
2) Add an extra point to the first and last intervals through
which spline must pass: not-a-knot splines.
Engineering Computation
Curve Fitting: Interpolation 16
Cubic Splines Computation
If set up cleverly, the 4n x 4n system of equations can be reduced to
solving an (n-1) x (n-1) tridiagonal system of equations.
Define a new set of unknowns: Let si = f "(xi) be the second
derivative of the cubic spline at interior point i, i = 1, ..., n-1. First
set up the n-1 equations to solve for curvatures, f "(x) at each of the
interior knots (see C&C Box 18.3):
(xi – xi-1) s i-1 + 2 (x i+1– x i-1) s i + (x i+1– x i) s i+1
6[f (x i1 )  f (x i )] 6[f (x i )  f (x i1)]


(x i1  xi )
(x i  x i1 )
Engineering Computation
Curve Fitting: Interpolation 17
Natural Cubic Splines
Convenient tridiagonal equations for natural splines:
These basic equations for the second derivatives can also be written
in terms of the distances hi = (xi+1–xi) between the points or knots
and of yi = f(xi):
For i = 0:
s0 = 0 [natural spline condition]
For i = 1:
2 (h0 +h1) s1 + h1s2
For i = 2 to n–2: hi-1s i-1
+ 2 (h i-1 + h i) s i
= RHS1
+ h is i+1 = RHSi
For i = n–1:
h n-2 sn-2 + 2 (h n-2 + h n-1) s n-1
For i = n:
sn = 0 [natural spline condition]
= RHSn-1
y
 yi yi  yi  1 
i

1
RHSi  6 


hi
hi  1 

Engineering Computation
Curve Fitting: Interpolation 18
Natural Cubic Splines
Noting that this is a triangular banded system of equations of
order n-1, we solve with a method which takes advantage of
this, i.e., the Thomas algorithm given in C&C Section 11.1.
Solve the triangular banded system, e.g., for n = 7 looks like:
s1
x
x
-
s2
x
x
x
-
s3
x
x
x
-
s4
x
x
x
-
s5
x
x
x
s6
x
x
Engineering Computation
Curve Fitting: Interpolation 19
Cubic Spline Interpolation
If we want to find specific interpolants, we do not need to
determine the cubic functions in all of the intervals, rather just the
interval in which the x lies. For the ith interval spanning [xi-1,xi]:
fi(x) = ai(x–xi-1)3 + bi(x–xi-1)2 + ci(x–xi-1) + di
in which:
si  si1
si 1
ai 
bi 
6h i1
2
yi  yi1 (2si1  si )hi1
ci 

h i1
6
di = yi-1
Obtain a different cubic polynomial for each interval [xi-1, xi].
However, first need to solve for the values of all the si .
Engineering Computation
Curve Fitting: Interpolation 20
Natural Cubic Splines: Example
Set up equations to solve for the unknown curvatures at each
interior knot for the following data:
i
xi
y = f(xi)
si
0
1
4.75
0
1
2
4
?
2
3
5.25
?
3
5
19.75
?
4
6
42
0
(xi – xi-1) si-1 + 2 (xi+1 – xi-1) si + (xi+1 – xi) si+1

6[f (x i1 )  f (x i )] 6[f (x i )  f (x i1)]

(x i1  xi )
(x i  x i1 )
Engineering Computation
Curve Fitting: Interpolation 21
Splines: Example (cont'd)
i
xi
y = f(xi)
0
1
4.75
0
1
2
4
?
2
3
5.25
?
3
5
19.75
?
4
6
42
0
For i = 1: (2 – 1) s0 + 2 (3 – 1) s1 + (3 – 2) s2
6[5.25  4] 6[4  4.75]


(3  2)
(2  1)
or
4 s1 +
For i = 2:
s1 +
For i = 3:
s2
= 12
6 s2 + 2 s3
= 36
2 s2 + 6 s3
= 90
si
Engineering Computation
Curve Fitting: Interpolation 22
Splines: Example (cont'd)
 4 1 0   s1  12 
1 6 2  s   36 

 2  
 0 2 6  s3  90 
Solve by any appropriate method (e.g., Thomas algorithm):
 s1   2.85 
  

s

0.59
 2 

s  14.80 
 3 

Engineering Computation
Curve Fitting: Interpolation 23
Splines: Example (cont'd)
Now estimate f(2.4) by using the spline:
i
xi
y = f(xi)
si
0
1
4.75
0
1
2
4
-
2.4
?
2
3
5.25
0.59
3
5
19.75
14.80
4
6
42
2.85
-
0
We need only solve for the ith=2nd interval, i.e.,
(x1 = 2) < x < (x2 = 3):
f2(x) = a2(x–x1)3 + b2(x–x1)2 + c2(x–x1) + d2
Engineering Computation
Curve Fitting: Interpolation 24
Splines: Example (cont'd)
Solving: f2(x) = a2(x–x1)3 + b2(x–x1)2 + c2(x–x1) + d2
s 2  s1 )
0.59  2.85
a2 

 -0.377
6(x 2  x1 )
6  32 
s1 2.85
b2  
 1.43
2
2
f (x 2 )  f (x1) [s1  s 2 ](x 2  x1)
c2 

(x 2  x1 )
6
5.25  4 [2(2.85)  0.59](3  2)


 0.202
(3  2)
6
d2 = f(x1) = 4.00
f2(2.4) = – 0.377(2.4–2)3 + 1.43(2.4–2)2 + 0.202(2.4–2) + 4.00
= 4.29
Engineering Computation
Curve Fitting: Interpolation 25
Splines: Example in Matlab
1
x = 0:10;
y = sin(x);
xx = 0:0.25:10;
yy = spline(x,y,xx);
plot(x,y,'o',xx,yy)
0.5
0
-0.5
-1
0
2
4
6
8
10
1
x = 0:10;
y = sin(x);
xx = 0:0.25:10;
yy = spline(x,[0 y 0],xx);
plot(x,y,'o',xx,yy)
0.5
0
-0.5
-1
0
2
4
6
8
10
Engineering Computation
Curve Fitting: Interpolation 26
npts = 10;
xy = [randn(1,npts); randn(1,npts)];
plot(xy(1,:),xy(2,:),'ro','LineWidth',2);
for n = 1:npts,
text(xy(1,n),xy(2,n),[' ' num2str(n)])
end
set(gca,'XTick',[],'YTick',[])
cv = cscvn(xy);
fnplt(cv,'r',2)
fnplt(cscvn(xy),'r',2)
1.5
5
7
1
3
2
6
4
0.5
0
9
10
1
-0.5
-1
8
-1.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2