Matlab Math
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
多项式和样条
在平面上给定 n个点 (xk,yk),可以唯一确定一个最多 n-
1次的多项式通过这些点，这个多项式叫插值多项式
插值多项式
P(xk ) = yk,k = 1,2,…,n
Lagrange插值形式
? ? ??
?
?
?
?
?
?
?
?
?
?k
k
kj jk
j y
xx
xx
xP )(
插值多项式例子
x =0:3;
y = [-5 –6 –1 16];
disp([x; y])
)16(
)6(
)2)(1(
)1(
)2(
)3)(1(
)6(
)2(
)3)(2(
)5(
)6(
)3)(2)(1(
)(
??
??
?
??
??
??
??
?
???
?
xxxxxx
xxxxxx
xP
52)( 3 ??? xxxP
Monomial基
nnnn cxcxcxcxP ????? ??? 12211,..)(
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
nnn
n
n
n
n
nn
nn
y
y
y
c
c
c
xxx
xxx
xxx
??
?
?????
?
?
2
1
2
1
21
2
2
2
1
2
1
2
1
1
1
1
1
1
Vandermonde矩阵
jnkjk xv ??,
x =0:3; y = [-5 –6 –1 16];
V = vander(x)
c = V\y’
52)( 3 ??? xxxP
Vandermonde矩阵是非
奇异的，但条件数是非
常坏的
Polyinterp(Lagrange插值形式 )
function v = polyinterp(x,y,u)
n = length(x);
for k = 1:n
w = ones(size(u));
for j = [1:k-1 k+1:n]
w = (u – x(j))./(x(k) – x(j)).*w;
end
v = v + w*y(k);
end
-0.5 0 0.5 1 1.5 2 2.5 3 3.5-10
-5
0
5
10
15
20
25
u = -.25:.01:3.25; v = polyinterp(x,y,u);
plot(x,y,’o’,u,v,’-’)
Polyinterp(符号运算 )
symx =sym(‘x’)
P = polyinterp(x,y,symx)
pretty(P)
P = simplify(P)
P =
x^3-2*x-5
0 1 2 3 4 5 6 76
8
10
12
14
16
18
20
22
Polyinterp(另外的例子 )
x = 1:6; y = [16 18 21 17 15 12];
u =,75:.05:6.25; v = polyinterp(x,y,u);
plot(x,y,’o’,u,v,’-’);
0 1 2 3 4 5 6 710
12
14
16
18
20
22
分片线性插值
x = 1:6; y = [16 18 21 17 15 12];
plot(x,y,’o’,u,v,’-’);
kk
kk
kk xx
yyxxyxL
?
????
?
?
1
1)()(
function [v,sigma] = piecelin(x,y,u)
d = diff(y)./diff(x); % First divided difference
% Find subinterval indices,x(k) <= u < u(k+1)
n = length(x);
k = ones(size(u));
for j = 2:n-1
k(u >= x(j)) = j;
end
% Evaluate interpolant
s = u - x(k);
v = y(k) + s.*d(k);
分片三次插值
kkkk dh
hssd
h
hssy
h
shshy
h
shsxP
2
2
12
2
3
223
13
32 )()(2323
)( ????????? ??
设 s = x-xk,h = hk
11
11
)(',)('
,)(,)(
??
??
??
??
kkkk
kkkk
dxPdxP
yxPyxP
Hermite插值（ osculatory插值）
pchiptx.m
splinetx.m
pchiptx
pchip,piecewise cubic Hermit interpolating polynomial,
Matlab中 pchip算法基于 Fritsch和 Carlson
1、如果左右导数是相反符号或有 0，则 dk=0
2,如果同号且区间长度相等，则是调和平均
???
?
???
? ??
? kkkd ??
11
2
11
1
3、如果同号且区间长度不相等，则是加权调和平均
1211
2
1
121
2,2
,
2
1
??
?
????
??
?
?
??
?
?
??
?
kkkk
kkk
hhwhhw
ww
d
ww
??
k
kk
k h
yy ?? ? 1?
0 1 2 3 4 5 6 710
12
14
16
18
20
22
三次样条
三次样条也是分片三次插值函数。
物理上的样条在满足插值限制的前提下，最小化势能。
数学上的样条必须满足二次导数连续，且满足插值限制。
参考文献：
A Practical Guide to Spline,Carl de Boor。
他也是 Matlab的 spline函数和 spline工具箱的作者。
三次样条
2
1 )46()26()126()(''
h
dhsdhsshxP kkk ?????? ??
0,,426)('' 1 ?????? ? sxxh ddxP k
k
kkk
k
?
kk
k
kkk
k hsxxh
ddxP ???????
?
?
?,,
246)(''
1
1
1
?
1
1
11,,246)(''
?
?
?? ???????
kk
k
kkk
k hsxxh
ddxP ?
kkkkkkkkkkk hhdhdhhdh ?? 111111 33)(2 ?????? ?????
kkkkk ddd ?? 334 111 ???? ???
二阶导
数连续
等距
,Not-a-knot”
在边界上，把两段合成一段，在等距的情况下：
2121 2
1
2
52 ?? ??? dd
121 2
5
2
12
??? ??? nnnn dd ??
这样，可以导致一个线性代数方程组,rAd ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
n
nnnnn
r
hh
hh
hh
r
r
d
d
d
d
1221
3223
2112
1
2
1
3,
??
??
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
????
221
2211
2233
1122
122
)(2
)(2
)(2
nnn
nnnn
hhh
hhhh
hhhh
hhhh
hhh
A
???
三次样条
rAd \?
然后我们可以计算出各个节点的一阶导数
0 2 4 6 810
12
14
16
18
20
22
分片线性插值
0 2 4 6 85
10
15
20
25
full degree多项式插值
0 2 4 6 810
12
14
16
18
20
22
保形状 Hermit插值
0 2 4 6 810
12
14
16
18
20
22
样条插值
分 析
上面的图描述了光滑性和局部单调性（保形状）的一种
折衷。
分片线性插值：保持单调性的，但光滑性比较差。
Full degree多项式插值：无限可微，但不保持形状，特
别是在端点的地方。
Pchip和 spline插值在这两个极端之间。样条比 pchip光滑，
样条的两阶导数连续，而 pchip一阶导数连续。不连续的
两阶导数隐含着不连续的曲率。人的眼睛可以检测出图
形上曲率的不连续。另一方面,pchip是保形状的，而样
条不一定保形状。
pchiptx 和 splinetx
pchiptx和 splinetx都是基于分片三次 Hermite插值，在每个区间
kkkk dh
hssd
h
hssy
h
shshy
h
shsxP
2
2
12
2
3
323
13
32 )()(2323
)( ????????? ??
kkkk dscssdyxP 32)( ????
2
11 2,23
h
ddb
h
ddc kkk
k
kkk
k
?? ?????? ??
用 monomial基：
pchiptx.m,
splinetx.m
请比较一下与 matlab中的 pchip
与 spline的区别
1 2 3 4 5 6-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
Interpolatation
interpgui
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
多项式和样条
在平面上给定 n个点 (xk,yk),可以唯一确定一个最多 n-
1次的多项式通过这些点，这个多项式叫插值多项式
插值多项式
P(xk ) = yk,k = 1,2,…,n
Lagrange插值形式
? ? ??
?
?
?
?
?
?
?
?
?
?k
k
kj jk
j y
xx
xx
xP )(
插值多项式例子
x =0:3;
y = [-5 –6 –1 16];
disp([x; y])
)16(
)6(
)2)(1(
)1(
)2(
)3)(1(
)6(
)2(
)3)(2(
)5(
)6(
)3)(2)(1(
)(
??
??
?
??
??
??
??
?
???
?
xxxxxx
xxxxxx
xP
52)( 3 ??? xxxP
Monomial基
nnnn cxcxcxcxP ????? ??? 12211,..)(
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
nnn
n
n
n
n
nn
nn
y
y
y
c
c
c
xxx
xxx
xxx
??
?
?????
?
?
2
1
2
1
21
2
2
2
1
2
1
2
1
1
1
1
1
1
Vandermonde矩阵
jnkjk xv ??,
x =0:3; y = [-5 –6 –1 16];
V = vander(x)
c = V\y’
52)( 3 ??? xxxP
Vandermonde矩阵是非
奇异的，但条件数是非
常坏的
Polyinterp(Lagrange插值形式 )
function v = polyinterp(x,y,u)
n = length(x);
for k = 1:n
w = ones(size(u));
for j = [1:k-1 k+1:n]
w = (u – x(j))./(x(k) – x(j)).*w;
end
v = v + w*y(k);
end
-0.5 0 0.5 1 1.5 2 2.5 3 3.5-10
-5
0
5
10
15
20
25
u = -.25:.01:3.25; v = polyinterp(x,y,u);
plot(x,y,’o’,u,v,’-’)
Polyinterp(符号运算 )
symx =sym(‘x’)
P = polyinterp(x,y,symx)
pretty(P)
P = simplify(P)
P =
x^3-2*x-5
0 1 2 3 4 5 6 76
8
10
12
14
16
18
20
22
Polyinterp(另外的例子 )
x = 1:6; y = [16 18 21 17 15 12];
u =,75:.05:6.25; v = polyinterp(x,y,u);
plot(x,y,’o’,u,v,’-’);
0 1 2 3 4 5 6 710
12
14
16
18
20
22
分片线性插值
x = 1:6; y = [16 18 21 17 15 12];
plot(x,y,’o’,u,v,’-’);
kk
kk
kk xx
yyxxyxL
?
????
?
?
1
1)()(
function [v,sigma] = piecelin(x,y,u)
d = diff(y)./diff(x); % First divided difference
% Find subinterval indices,x(k) <= u < u(k+1)
n = length(x);
k = ones(size(u));
for j = 2:n-1
k(u >= x(j)) = j;
end
% Evaluate interpolant
s = u - x(k);
v = y(k) + s.*d(k);
分片三次插值
kkkk dh
hssd
h
hssy
h
shshy
h
shsxP
2
2
12
2
3
223
13
32 )()(2323
)( ????????? ??
设 s = x-xk,h = hk
11
11
)(',)('
,)(,)(
??
??
??
??
kkkk
kkkk
dxPdxP
yxPyxP
Hermite插值（ osculatory插值）
pchiptx.m
splinetx.m
pchiptx
pchip,piecewise cubic Hermit interpolating polynomial,
Matlab中 pchip算法基于 Fritsch和 Carlson
1、如果左右导数是相反符号或有 0，则 dk=0
2,如果同号且区间长度相等，则是调和平均
???
?
???
? ??
? kkkd ??
11
2
11
1
3、如果同号且区间长度不相等，则是加权调和平均
1211
2
1
121
2,2
,
2
1
??
?
????
??
?
?
??
?
?
??
?
kkkk
kkk
hhwhhw
ww
d
ww
??
k
kk
k h
yy ?? ? 1?
0 1 2 3 4 5 6 710
12
14
16
18
20
22
三次样条
三次样条也是分片三次插值函数。
物理上的样条在满足插值限制的前提下，最小化势能。
数学上的样条必须满足二次导数连续，且满足插值限制。
参考文献：
A Practical Guide to Spline,Carl de Boor。
他也是 Matlab的 spline函数和 spline工具箱的作者。
三次样条
2
1 )46()26()126()(''
h
dhsdhsshxP kkk ?????? ??
0,,426)('' 1 ?????? ? sxxh ddxP k
k
kkk
k
?
kk
k
kkk
k hsxxh
ddxP ???????
?
?
?,,
246)(''
1
1
1
?
1
1
11,,246)(''
?
?
?? ???????
kk
k
kkk
k hsxxh
ddxP ?
kkkkkkkkkkk hhdhdhhdh ?? 111111 33)(2 ?????? ?????
kkkkk ddd ?? 334 111 ???? ???
二阶导
数连续
等距
,Not-a-knot”
在边界上，把两段合成一段，在等距的情况下：
2121 2
1
2
52 ?? ??? dd
121 2
5
2
12
??? ??? nnnn dd ??
这样，可以导致一个线性代数方程组,rAd ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
n
nnnnn
r
hh
hh
hh
r
r
d
d
d
d
1221
3223
2112
1
2
1
3,
??
??
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
????
221
2211
2233
1122
122
)(2
)(2
)(2
nnn
nnnn
hhh
hhhh
hhhh
hhhh
hhh
A
???
三次样条
rAd \?
然后我们可以计算出各个节点的一阶导数
0 2 4 6 810
12
14
16
18
20
22
分片线性插值
0 2 4 6 85
10
15
20
25
full degree多项式插值
0 2 4 6 810
12
14
16
18
20
22
保形状 Hermit插值
0 2 4 6 810
12
14
16
18
20
22
样条插值
分 析
上面的图描述了光滑性和局部单调性（保形状）的一种
折衷。
分片线性插值：保持单调性的，但光滑性比较差。
Full degree多项式插值：无限可微，但不保持形状，特
别是在端点的地方。
Pchip和 spline插值在这两个极端之间。样条比 pchip光滑，
样条的两阶导数连续，而 pchip一阶导数连续。不连续的
两阶导数隐含着不连续的曲率。人的眼睛可以检测出图
形上曲率的不连续。另一方面,pchip是保形状的，而样
条不一定保形状。
pchiptx 和 splinetx
pchiptx和 splinetx都是基于分片三次 Hermite插值，在每个区间
kkkk dh
hssd
h
hssy
h
shshy
h
shsxP
2
2
12
2
3
323
13
32 )()(2323
)( ????????? ??
kkkk dscssdyxP 32)( ????
2
11 2,23
h
ddb
h
ddc kkk
k
kkk
k
?? ?????? ??
用 monomial基：
pchiptx.m,
splinetx.m
请比较一下与 matlab中的 pchip
与 spline的区别
1 2 3 4 5 6-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
Interpolatation
interpgui
