Introduction to Scientific
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter3 Piecewise Polynomial
Interpolation
? Piecewise Linear Interpolation
? Piecewise Cubic Hermite Interpolation
? Cubic Splines
0 1 2 3 4 5 6 7-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Piecewise Linear Interpolation
)()( iiii zzbazL ???),( ii yx
ii
ii
iii xx
yybya
?
???
?
?
1
1,
?
?
?
?
?
?
?
??
??
??
?
?? nnn
xzxzL
xzxzL
xzxzL
zL
11
322
211
if)(
if)(
if)(
)(
??
Piecewise Linear
function [a,b] = pwL(x,y)
n = length(x);
a = y(1:n-1);
b = diff(y),/ diff(x);
z=linspace(0,1,9);
[a,b]=pwL(z,sin(2*pi*z));
for i=1:n-1
b(i)=(y(i+1)- y(i))/ (x(i+1)- x(i));
end
b=(y(2:n)-y(1:n-1))./(x(2:n)-x(1:n-
1))
b=diff(y)./diff(x)
Test code
Evalution
],[ ???z?)( ?zL
Problem,
],,[ 1 ?? ? ixxz ii
if z==x(n)
i =n-1;
else
i=sum(x<=z);
end
Binary search
mid = floor((Left+Right)/2);
If z<x(mid)
Right=mid;
Else
Left =mid;
end
10000?N 1 3, 3 l o g 2 ?? Nn
function i = Locate(x,z,g)
% g (1<=g<=n-1) is an optional input parameter
% search for i begins,the value i=g is tried.
if nargin==3
if (x(g)<=z) & (z<=x(g+1))
i = g; return
end;end
n = length(x);
if z = = x(n)
i = n-1;
else
Left = 1; Right = n;
while Right > Left+1
Binary_search
end
g,guss
function LVals = pwLEval(a,b,x,zVals)
% Evaluates the piecewise linear polynomial defined by the column
%(n-1)-vectors
m = length(zVals);
LVals = zeros(m,1);
g = 1;
for j=1:m
i = Locate(x,zVals(j),g);
LVals(j) = a(i) + b(i)*(zVals(j)-x(i));
g = i;
end
m-vector
6
04.2)^9.(
1
01.)3.(
1
2 ?????? xx
A priori determination of breakpoints
],[ 1?? ii xxz
))((
2
)()()(
1
)2(
????? ii xzxz
fzLzf ?
8
|)()(|
2
2 hMzLzf ?? ],[ ???z
??? ??
?
??
?
?
?
???? 2222
188
|)()(|
n
MhMzLzf
??? 8/)(1 2Mn ???
Static
Adaptive Piecewise Linear Interpolation
Problem
Adaptive Piecewise Linear Interpolation
],[ xRxLacceptable
?????
?
??
?
? ?
2
)()(
2
xRfxLfxRxLf
m i nhxLxR ??
or
Piecewise Cubic Hermit
Interpolation
2
3,
2
3,,0
4321
???? ????? xxxx
0??
Hermit cubic interpolant
)()()()()( 22 RLLL xzxzdxzcxzbazq ????????
RRLL
RRLL
sxqsxq
yxqyxq
??
??
)(',)('
)(,)(
LR
LRLR
L
LLR
LL
LL
xx
yy
x
yy
y
x
yss
d
x
sy
c
sbya
?
?
?
?
?
?
?
??
?
?
?
?
??
'
'2
'
],[)( 4)4( RL xxzMzf ??
LR xxhh
Mzqzf ????,
3 8 4
)()( 44
410,10
e r r o re r r o rhh ??
)()()()()( 122 ????????? iiiiiiiii xzxzdxzcxzbazq
?
?
?
?
?
?
?
??
??
??
?
?? nnn
xzxzq
xzxzq
xzxzq
zC
11
322
211
if)(
if)(
if)(
)(
??
Piecewise cubic polynomial
nisxCyxC iiii,1,)(',)( ???
For i=1:n-1
[a(i),b(i),c(i),d(i)=Hcubic(x(i),y(i),s(i),x(i+1),y(i+1),s(i+1))
End
function [a,b,c,d] = pwC(x,y,s)
% Piecewise cubic Hermite interpolation.
n = length(x);
a = y(1:n-1);
b = s(1:n-1);
Dx = diff(x);
Dy = diff(y);
yp = Dy,/ Dx;
c = (yp - s(1:n-1)),/ Dx;
d = (s(2:n) + s(1:n-1) - 2*yp),/ (Dx.* Dx);
x,y,s,vector
Evaluation
function Cvals = pwCEval(a,b,c,d,x,zVals)
m = length(zVals);
Cvals = zeros(m,1);
g=1;
for j=1:m
i = Locate(x,zVals(j),g);
Cvals(j) = d(i)*(zVals(j)-x(i+1)) + c(i);
Cvals(j) = Cvals(j)*(zVals(j)-x(i)) + b(i);
Cvals(j) = Cvals(j)*(zVals(j)-x(i)) + a(i);
g = i;
end
Locate
Cubic version of
HornerN
Cubic spline interpolant
Given with
find a piecewise cubic interpolant with the
property that S,S' and S'' are continuous.
),(),.,,,,( 11 nn yxyx
,...1 ?? ???? nxx
Why
need
spline?
Continuity at the interior knots
? ?
)()(
'2
)(
'
)()(
1
2
2
1
2
?
?
??
?
??
?
?
?
?
????
ii
i
iii
i
i
ii
iiii
xzxz
x
yss
xz
x
sy
xzsyzq
ii
ii
i xx
yyy
?
??
?
?
1
1'
? ? )](2)(4[
'2'2)(''
12
1
?
? ???
?
???
?
??
ii
i
iii
i
ii
i xzxzx
yss
x
syzq
? ? )](2)(4[
22)(
212
1
'
121
1
1
'
1''
1 ??
?
???
?
??
? ????
???
?
??
ii
i
iii
i
ii
i xzxzx
yss
x
syzq
)32(2)( '11'' iii
i
ii yssxxq ???? ??
)23(2)( 21' 1
1
1
''
1 ???
?
?? ???? iii
i
ii ssyxxq
? ? 2:1,3
)(2
'
1
'
1
2111
??????
???????
??
????
niyxyx
sxsxxsx
iiii
iiiiiii
)(3)(21 '21'123122112 yxyxsxsxxsxi ?????????????
)(3)(22 '32'234233223 yxyxsxsxxsxi ?????????????
)(3)(23 '43'345344334 yxyxsxsxxsxi ?????????????
)(3)(24 '54'456455445 yxyxsxsxxsxi ?????????????
)(3)(25 '65'567566556 yxyxsxsxxsxi ?????????????
r
sxyxyx
yxyx
yxyx
yxyx
sxyxyx
s
s
s
s
s
TTs ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
???
???
???
?????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
75
'
65
'
56
'
54
'
45
'
43
'
34
'
32
'
23
12
'
21
'
12
6
5
4
3
2
)(3
)(3
)(3
)(3
)(3
)6:2(
? ?
? ?
? ?
? ?
? ??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
?????
?????
?????
????
?
655
4545
3434
2323
121
2
2
2
2
2
xxx
xxxx
xxxx
xxxx
xxx
T
tridiagonal
? ?
? ?
)2()2(
2,23,2
3232
2323
1211
00
20
2
00
???
????
????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
?????
?
nn
nnnn
nnnn
tt
xxxx
xxxx
tt
T
?
?????
?
?
nss,1
choice
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
?
?
????
2
'
23
'
32
'
32
'
23
1
)(3
)(3
n
nnnn
r
yxyx
yxyx
r
r ?
rTs ?
n=length(x);
Dx=diff(x);
yp=diff(y)./Dx;
T=zeros(n-2,n-2);
r=zeros(n-2,1);
for i=2:n-3
T(i,i)=2(Dx(i)+Dx(i+1));
T(i,i-1)=Dx(i+1);
T(i,i+1)=Dx(i);
r(i)=3(Dx(i+1)*yp(i)+Dx(i)*yp(i+1));
end
Help diag
The complete spline
RnL ss ?? ??,1
? ?'21'12312212 3)(2 yxyxsxsxxx L ??????????? ?
? ?' 12' 21
211221
3
)(2
????
??????
????
???????
nnnn
Rnnnnnn
yxyx
xsxxsx ?
? ? Lxyxyxsxsxx ?2'21'1231221 3)(2 ???????????
? ? Rnnnnn
nnnnn
xyxyx
sxxsx
?2' 12' 21
11221
3
)(2
?????
?????
??????
?????
T(1,1)=2*(Dx(1)+Dx(2));
T(1,2)=Dx(1);
r(1)=3*(Dx(2)*yp(1)+Dx(1)*yp(2))-Dx(2)*muL;
T(n-2,n-2)=2*(Dx(n-2)+Dx(n-1));
T(n-2,n-3)=Dx(n-1);
r(n-2)=3*(Dx(n-1)*yp(n-2)+Dx(n-2)*yp(n-1))-Dx(n-2)*muR;
s=[muL;T\r(1:n-2);muR];
The Natural Spline
)( 1''1 xqL ??
? ? )](2)(4[
'2'2)(''
12
1
?
? ???
?
???
?
??
ii
i
iii
i
ii
i xzxzx
yss
x
syzq
1
121
1
11 '22'2
x
yss
x
sy
L ?
???
?
???
?
?
??
?
? ????
12
'
11 232
1 xsys L?
?
?
??
?
? ????
??? 11
'
1 232
1
n
R
nnn xsys
?
)('' 1 nnR xq ??? 0?? RL ??
T(1,1) = 2*Dx(1) + 1.5*Dx(2);
T(1,2) = Dx(1);
r(1) = 1.5*Dx(2)*yp(1) + 3*Dx(1)*yp(2) + Dx(1)*Dx(2)*muL/4;
T(n-2,n-2) = 1.5*Dx(n-2)+2*Dx(n-1);
T(n-2,n-3) = Dx(n-1);
r(n-2) = 3*Dx(n-1)*yp(n-2) + 1.5*Dx(n-2)*yp(n-1)- …
Dx(n-1)*Dx(n-2)*muR/4;
stilde = T\r;
s1 = (3*yp(1) - stilde(1) - muL*Dx(1)/2)/2;
sn = (3*yp(n-1) - stilde(n-2) + muR*Dx(n-1)/2)/2;
s = [s1;stilde;sn];
The Not-a-Knot Spline
)()( 2'''22'''1 xqxq ?
? ? ? ? 22
'
232
2
1
'
121 22
x
yss
x
yss
?
??
?
?
??
)2(2 '232
2
2
1'
121 yssx
x
yss ????
?
?
??
?
?
?
?
????
)2(2 ' 222
2
2
1'
11 ???
?
?
?? ????
?
?
??
?
?
?
?
???? nnn
n
n
nnn yssx
x
yss
i d e n t i ca l a r e )( a n d )( 21 xqxq
q = Dx(1)*Dx(1)/Dx(2);
T(1,1) = 2*Dx(1) +Dx(2) + q;
T(1,2) = Dx(1) + q;
r(1) = Dx(2)*yp(1) + Dx(1)*yp(2)+2*yp(2)*(q+Dx(1));
q = Dx(n-1)*Dx(n-1)/Dx(n-2);
T(n-2,n-2) = 2*Dx(n-1) + Dx(n-2)+q;
T(n-2,n-3) = Dx(n-1)+q;
r(n-2) = Dx(n-1)*yp(n-2) + Dx(n-2)*yp(n-1) …
+2*yp(n-2)*(Dx(n-1)+q);
stilde = T\r;
s1 = -stilde(1)+2*yp(1);
s1 = s1 + ((Dx(1)/Dx(2))^2)*(stilde(1)+stilde(2)-2*yp(2));
sn = -stilde(n-2) +2*yp(n-1);
sn = sn+((Dx(n-1)/Dx(n-2))^2)*(stilde(n-3) …
+stilde(n-2)-2*yp(n-2));
s = [s1;stilde;sn];
function [a,b,c,d] = CubicSpline(x,y,derivative,muL,muR)
% Cubic spline interpolation with prescribed end conditions.
% Usage:
% [a,b,c,d] = CubicSpline(x,y,1,muL,muR)
% S'(x(1)) = muL,S'(x(n)) = muR
% [a,b,c,d] = CubicSpline(x,y,2,muL,muR)
% S''(x(1)) = muL,S''(x(n)) = muR
% [a,b,c,d] = CubicSpline(x,y)
% S'''(z) continuous at x(2) and x(n-1)
0 0.5 110
-12
10 -10
10 -8
10 -6
10 -4
Knot Spacing = 0.050
0 0.5 110
-12
10 -10
10 -8
10 -6
10 -4
Knot Spacing = 0.005
Not-a-knot spline error
0 1 2 3 4 5 6 7 8 9 10-1
-0.5
0
0.5
1
1.5
Bad end conditions
-5 -4 -3 -2 -1 0 1 2 3 4 5-1.5
-1
-0.5
0
0.5
1
1.5
n = 9 Spline Interpolant of atan(x)
Matlab Spline Tools
z=linspace(-5,5);
x=linspace(-5,5,9);
y=atan(x);
Svals=spline(x,y,z);
plot(z,Svals)
Not-a-knot spline
pp-representation
x=linspace(-5,5,9);
y=atan(x);
S=spline(x,y);
z=linspace(-5,5);
Svals=ppval(S,z)
plot(z,Svals)
pp-representation
[x,rho,L,k]=unmkpp(S)
The coefficients of the local
polynomials are assembled in an
L-by-k matrix rho
3
3,
2
2,3,4,)()()()( iiiiiii xzxzxzzS ??????? ????
1
,)(j)r h o ( i,
????? jk
iji xx?
3
3,
2
2,3,4,)()()()( iiiiiii xzxzxzzS ??????? ????
2
3,2,3,)(3)(2)(' iiiii xzxzzS ????? ???
drho = [3*rho(:,1) 2*rho(:,2) rho(:,3)];
dS = mkpp(x,drho);
V = PPVAL(PP,XX) returns the value at the points XX
of the piecewise
polynomial contained in PP,as constructed by
SPLINE or the spline utility
z = linspace(-5,5);
Svals = ppval(S,z);
dSvals = ppval(dS,z);
-5 -4 -3 -2 -1 0 1 2 3 4 50
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Derivative of n = 9 Spline Interpolant of atan(x)
n = 9;x = linspace(-5,5,n);y = atan(x);S = spline(x,y);
[x,rho,L,k] = unmkpp(S);
drho = [3*rho(:,1) 2*rho(:,2) rho(:,3)];
dS = mkpp(x,drho);
z = linspace(-5,5);
Svals = ppval(S,z);
dSvals = ppval(dS,z);
atanvals = atan(z);
Figure;plot(z,atanvals,z,Svals,x,y,'*');
title(sprintf('n = %2.0f Spline Interpolant of atan(x)',n))
datanvals = ones(size(z))./(1 + z.*z);
figure
plot(z,datanvals,z,dSvals)
title(sprintf('Derivative of n = %2.0f Spline Interpolant of atan(x)',n))
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter3 Piecewise Polynomial
Interpolation
? Piecewise Linear Interpolation
? Piecewise Cubic Hermite Interpolation
? Cubic Splines
0 1 2 3 4 5 6 7-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Piecewise Linear Interpolation
)()( iiii zzbazL ???),( ii yx
ii
ii
iii xx
yybya
?
???
?
?
1
1,
?
?
?
?
?
?
?
??
??
??
?
?? nnn
xzxzL
xzxzL
xzxzL
zL
11
322
211
if)(
if)(
if)(
)(
??
Piecewise Linear
function [a,b] = pwL(x,y)
n = length(x);
a = y(1:n-1);
b = diff(y),/ diff(x);
z=linspace(0,1,9);
[a,b]=pwL(z,sin(2*pi*z));
for i=1:n-1
b(i)=(y(i+1)- y(i))/ (x(i+1)- x(i));
end
b=(y(2:n)-y(1:n-1))./(x(2:n)-x(1:n-
1))
b=diff(y)./diff(x)
Test code
Evalution
],[ ???z?)( ?zL
Problem,
],,[ 1 ?? ? ixxz ii
if z==x(n)
i =n-1;
else
i=sum(x<=z);
end
Binary search
mid = floor((Left+Right)/2);
If z<x(mid)
Right=mid;
Else
Left =mid;
end
10000?N 1 3, 3 l o g 2 ?? Nn
function i = Locate(x,z,g)
% g (1<=g<=n-1) is an optional input parameter
% search for i begins,the value i=g is tried.
if nargin==3
if (x(g)<=z) & (z<=x(g+1))
i = g; return
end;end
n = length(x);
if z = = x(n)
i = n-1;
else
Left = 1; Right = n;
while Right > Left+1
Binary_search
end
g,guss
function LVals = pwLEval(a,b,x,zVals)
% Evaluates the piecewise linear polynomial defined by the column
%(n-1)-vectors
m = length(zVals);
LVals = zeros(m,1);
g = 1;
for j=1:m
i = Locate(x,zVals(j),g);
LVals(j) = a(i) + b(i)*(zVals(j)-x(i));
g = i;
end
m-vector
6
04.2)^9.(
1
01.)3.(
1
2 ?????? xx
A priori determination of breakpoints
],[ 1?? ii xxz
))((
2
)()()(
1
)2(
????? ii xzxz
fzLzf ?
8
|)()(|
2
2 hMzLzf ?? ],[ ???z
??? ??
?
??
?
?
?
???? 2222
188
|)()(|
n
MhMzLzf
??? 8/)(1 2Mn ???
Static
Adaptive Piecewise Linear Interpolation
Problem
Adaptive Piecewise Linear Interpolation
],[ xRxLacceptable
?????
?
??
?
? ?
2
)()(
2
xRfxLfxRxLf
m i nhxLxR ??
or
Piecewise Cubic Hermit
Interpolation
2
3,
2
3,,0
4321
???? ????? xxxx
0??
Hermit cubic interpolant
)()()()()( 22 RLLL xzxzdxzcxzbazq ????????
RRLL
RRLL
sxqsxq
yxqyxq
??
??
)(',)('
)(,)(
LR
LRLR
L
LLR
LL
LL
xx
yy
x
yy
y
x
yss
d
x
sy
c
sbya
?
?
?
?
?
?
?
??
?
?
?
?
??
'
'2
'
],[)( 4)4( RL xxzMzf ??
LR xxhh
Mzqzf ????,
3 8 4
)()( 44
410,10
e r r o re r r o rhh ??
)()()()()( 122 ????????? iiiiiiiii xzxzdxzcxzbazq
?
?
?
?
?
?
?
??
??
??
?
?? nnn
xzxzq
xzxzq
xzxzq
zC
11
322
211
if)(
if)(
if)(
)(
??
Piecewise cubic polynomial
nisxCyxC iiii,1,)(',)( ???
For i=1:n-1
[a(i),b(i),c(i),d(i)=Hcubic(x(i),y(i),s(i),x(i+1),y(i+1),s(i+1))
End
function [a,b,c,d] = pwC(x,y,s)
% Piecewise cubic Hermite interpolation.
n = length(x);
a = y(1:n-1);
b = s(1:n-1);
Dx = diff(x);
Dy = diff(y);
yp = Dy,/ Dx;
c = (yp - s(1:n-1)),/ Dx;
d = (s(2:n) + s(1:n-1) - 2*yp),/ (Dx.* Dx);
x,y,s,vector
Evaluation
function Cvals = pwCEval(a,b,c,d,x,zVals)
m = length(zVals);
Cvals = zeros(m,1);
g=1;
for j=1:m
i = Locate(x,zVals(j),g);
Cvals(j) = d(i)*(zVals(j)-x(i+1)) + c(i);
Cvals(j) = Cvals(j)*(zVals(j)-x(i)) + b(i);
Cvals(j) = Cvals(j)*(zVals(j)-x(i)) + a(i);
g = i;
end
Locate
Cubic version of
HornerN
Cubic spline interpolant
Given with
find a piecewise cubic interpolant with the
property that S,S' and S'' are continuous.
),(),.,,,,( 11 nn yxyx
,...1 ?? ???? nxx
Why
need
spline?
Continuity at the interior knots
? ?
)()(
'2
)(
'
)()(
1
2
2
1
2
?
?
??
?
??
?
?
?
?
????
ii
i
iii
i
i
ii
iiii
xzxz
x
yss
xz
x
sy
xzsyzq
ii
ii
i xx
yyy
?
??
?
?
1
1'
? ? )](2)(4[
'2'2)(''
12
1
?
? ???
?
???
?
??
ii
i
iii
i
ii
i xzxzx
yss
x
syzq
? ? )](2)(4[
22)(
212
1
'
121
1
1
'
1''
1 ??
?
???
?
??
? ????
???
?
??
ii
i
iii
i
ii
i xzxzx
yss
x
syzq
)32(2)( '11'' iii
i
ii yssxxq ???? ??
)23(2)( 21' 1
1
1
''
1 ???
?
?? ???? iii
i
ii ssyxxq
? ? 2:1,3
)(2
'
1
'
1
2111
??????
???????
??
????
niyxyx
sxsxxsx
iiii
iiiiiii
)(3)(21 '21'123122112 yxyxsxsxxsxi ?????????????
)(3)(22 '32'234233223 yxyxsxsxxsxi ?????????????
)(3)(23 '43'345344334 yxyxsxsxxsxi ?????????????
)(3)(24 '54'456455445 yxyxsxsxxsxi ?????????????
)(3)(25 '65'567566556 yxyxsxsxxsxi ?????????????
r
sxyxyx
yxyx
yxyx
yxyx
sxyxyx
s
s
s
s
s
TTs ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
???
???
???
?????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
75
'
65
'
56
'
54
'
45
'
43
'
34
'
32
'
23
12
'
21
'
12
6
5
4
3
2
)(3
)(3
)(3
)(3
)(3
)6:2(
? ?
? ?
? ?
? ?
? ??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
?????
?????
?????
????
?
655
4545
3434
2323
121
2
2
2
2
2
xxx
xxxx
xxxx
xxxx
xxx
T
tridiagonal
? ?
? ?
)2()2(
2,23,2
3232
2323
1211
00
20
2
00
???
????
????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
?????
?
nn
nnnn
nnnn
tt
xxxx
xxxx
tt
T
?
?????
?
?
nss,1
choice
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
?
?
????
2
'
23
'
32
'
32
'
23
1
)(3
)(3
n
nnnn
r
yxyx
yxyx
r
r ?
rTs ?
n=length(x);
Dx=diff(x);
yp=diff(y)./Dx;
T=zeros(n-2,n-2);
r=zeros(n-2,1);
for i=2:n-3
T(i,i)=2(Dx(i)+Dx(i+1));
T(i,i-1)=Dx(i+1);
T(i,i+1)=Dx(i);
r(i)=3(Dx(i+1)*yp(i)+Dx(i)*yp(i+1));
end
Help diag
The complete spline
RnL ss ?? ??,1
? ?'21'12312212 3)(2 yxyxsxsxxx L ??????????? ?
? ?' 12' 21
211221
3
)(2
????
??????
????
???????
nnnn
Rnnnnnn
yxyx
xsxxsx ?
? ? Lxyxyxsxsxx ?2'21'1231221 3)(2 ???????????
? ? Rnnnnn
nnnnn
xyxyx
sxxsx
?2' 12' 21
11221
3
)(2
?????
?????
??????
?????
T(1,1)=2*(Dx(1)+Dx(2));
T(1,2)=Dx(1);
r(1)=3*(Dx(2)*yp(1)+Dx(1)*yp(2))-Dx(2)*muL;
T(n-2,n-2)=2*(Dx(n-2)+Dx(n-1));
T(n-2,n-3)=Dx(n-1);
r(n-2)=3*(Dx(n-1)*yp(n-2)+Dx(n-2)*yp(n-1))-Dx(n-2)*muR;
s=[muL;T\r(1:n-2);muR];
The Natural Spline
)( 1''1 xqL ??
? ? )](2)(4[
'2'2)(''
12
1
?
? ???
?
???
?
??
ii
i
iii
i
ii
i xzxzx
yss
x
syzq
1
121
1
11 '22'2
x
yss
x
sy
L ?
???
?
???
?
?
??
?
? ????
12
'
11 232
1 xsys L?
?
?
??
?
? ????
??? 11
'
1 232
1
n
R
nnn xsys
?
)('' 1 nnR xq ??? 0?? RL ??
T(1,1) = 2*Dx(1) + 1.5*Dx(2);
T(1,2) = Dx(1);
r(1) = 1.5*Dx(2)*yp(1) + 3*Dx(1)*yp(2) + Dx(1)*Dx(2)*muL/4;
T(n-2,n-2) = 1.5*Dx(n-2)+2*Dx(n-1);
T(n-2,n-3) = Dx(n-1);
r(n-2) = 3*Dx(n-1)*yp(n-2) + 1.5*Dx(n-2)*yp(n-1)- …
Dx(n-1)*Dx(n-2)*muR/4;
stilde = T\r;
s1 = (3*yp(1) - stilde(1) - muL*Dx(1)/2)/2;
sn = (3*yp(n-1) - stilde(n-2) + muR*Dx(n-1)/2)/2;
s = [s1;stilde;sn];
The Not-a-Knot Spline
)()( 2'''22'''1 xqxq ?
? ? ? ? 22
'
232
2
1
'
121 22
x
yss
x
yss
?
??
?
?
??
)2(2 '232
2
2
1'
121 yssx
x
yss ????
?
?
??
?
?
?
?
????
)2(2 ' 222
2
2
1'
11 ???
?
?
?? ????
?
?
??
?
?
?
?
???? nnn
n
n
nnn yssx
x
yss
i d e n t i ca l a r e )( a n d )( 21 xqxq
q = Dx(1)*Dx(1)/Dx(2);
T(1,1) = 2*Dx(1) +Dx(2) + q;
T(1,2) = Dx(1) + q;
r(1) = Dx(2)*yp(1) + Dx(1)*yp(2)+2*yp(2)*(q+Dx(1));
q = Dx(n-1)*Dx(n-1)/Dx(n-2);
T(n-2,n-2) = 2*Dx(n-1) + Dx(n-2)+q;
T(n-2,n-3) = Dx(n-1)+q;
r(n-2) = Dx(n-1)*yp(n-2) + Dx(n-2)*yp(n-1) …
+2*yp(n-2)*(Dx(n-1)+q);
stilde = T\r;
s1 = -stilde(1)+2*yp(1);
s1 = s1 + ((Dx(1)/Dx(2))^2)*(stilde(1)+stilde(2)-2*yp(2));
sn = -stilde(n-2) +2*yp(n-1);
sn = sn+((Dx(n-1)/Dx(n-2))^2)*(stilde(n-3) …
+stilde(n-2)-2*yp(n-2));
s = [s1;stilde;sn];
function [a,b,c,d] = CubicSpline(x,y,derivative,muL,muR)
% Cubic spline interpolation with prescribed end conditions.
% Usage:
% [a,b,c,d] = CubicSpline(x,y,1,muL,muR)
% S'(x(1)) = muL,S'(x(n)) = muR
% [a,b,c,d] = CubicSpline(x,y,2,muL,muR)
% S''(x(1)) = muL,S''(x(n)) = muR
% [a,b,c,d] = CubicSpline(x,y)
% S'''(z) continuous at x(2) and x(n-1)
0 0.5 110
-12
10 -10
10 -8
10 -6
10 -4
Knot Spacing = 0.050
0 0.5 110
-12
10 -10
10 -8
10 -6
10 -4
Knot Spacing = 0.005
Not-a-knot spline error
0 1 2 3 4 5 6 7 8 9 10-1
-0.5
0
0.5
1
1.5
Bad end conditions
-5 -4 -3 -2 -1 0 1 2 3 4 5-1.5
-1
-0.5
0
0.5
1
1.5
n = 9 Spline Interpolant of atan(x)
Matlab Spline Tools
z=linspace(-5,5);
x=linspace(-5,5,9);
y=atan(x);
Svals=spline(x,y,z);
plot(z,Svals)
Not-a-knot spline
pp-representation
x=linspace(-5,5,9);
y=atan(x);
S=spline(x,y);
z=linspace(-5,5);
Svals=ppval(S,z)
plot(z,Svals)
pp-representation
[x,rho,L,k]=unmkpp(S)
The coefficients of the local
polynomials are assembled in an
L-by-k matrix rho
3
3,
2
2,3,4,)()()()( iiiiiii xzxzxzzS ??????? ????
1
,)(j)r h o ( i,
????? jk
iji xx?
3
3,
2
2,3,4,)()()()( iiiiiii xzxzxzzS ??????? ????
2
3,2,3,)(3)(2)(' iiiii xzxzzS ????? ???
drho = [3*rho(:,1) 2*rho(:,2) rho(:,3)];
dS = mkpp(x,drho);
V = PPVAL(PP,XX) returns the value at the points XX
of the piecewise
polynomial contained in PP,as constructed by
SPLINE or the spline utility
z = linspace(-5,5);
Svals = ppval(S,z);
dSvals = ppval(dS,z);
-5 -4 -3 -2 -1 0 1 2 3 4 50
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Derivative of n = 9 Spline Interpolant of atan(x)
n = 9;x = linspace(-5,5,n);y = atan(x);S = spline(x,y);
[x,rho,L,k] = unmkpp(S);
drho = [3*rho(:,1) 2*rho(:,2) rho(:,3)];
dS = mkpp(x,drho);
z = linspace(-5,5);
Svals = ppval(S,z);
dSvals = ppval(dS,z);
atanvals = atan(z);
Figure;plot(z,atanvals,z,Svals,x,y,'*');
title(sprintf('n = %2.0f Spline Interpolant of atan(x)',n))
datanvals = ones(size(z))./(1 + z.*z);
figure
plot(z,datanvals,z,dSvals)
title(sprintf('Derivative of n = %2.0f Spline Interpolant of atan(x)',n))
