第5章〓样条性质
5．1引言

本章讨论样条曲线的几个常见性质，以及如何从样条方程计算这些性质［Mathews，2004］。第一个性质是称为最小点和最大点的样条曲线的关键点(有时也称临界点或极值点)。此外，对于3次或以上的样条曲线，拐点(POI)很重要，这是曲率从凹变为凸或从凸变为凹的位置。第二个性质是样条曲线的切线和法线。曲线的切线是曲线方程的导数，而法线是垂直于切线的直线。切线和法线可以用直线方程表示，或者若它们是针对曲线上的特定点计算的，则可以表示为特定的向量。对参数曲线和隐式曲线都讨论了切线和法线的计算。第三个性质是任意两个给定点之间的样条曲线的长度，这可以通过空间曲线和参数曲线来计算。曲线段的长度通过对许多小线段求和来近似。第四个性质是曲线下的面积，本章对于y=f(x)和x=f(y)形式的曲线都进行了讨论。该面积的计算方法是考虑曲线下的一个非常薄的矩形区域，然后对所有这些矩形区域求和。对此的扩展是计算由两条曲线界定的面积，这是上曲线下的面积减去下曲线下的面积。第五个性质是曲线下区域的质心，它是使用矩计算的。接下来讨论各种插值和曲线拟合的方法。当我们对某些给定数据点之间的某个中间值感兴趣时，就会进行插值。通过考虑连接相邻数据点的直线，插值可以是线性的，也可以通过考虑连接数据点的更高次曲线来非线性实现。本章还为此目的提到了内置的MATLAB函数。曲线拟合试图将多项式曲线拟合到给定的数据点并估计多项式的系数。这样可以通过多项式函数表示任意数据。本章最后回顾了MATLAB中包含的几个二维绘图函数。这些函数大致分为两种类型： 一种使用符号变量，另一种使用值的集合。这些函数也可以是显式或隐式或参数化的。其中，一些函数的参数可用于指定颜色和透明度值。

5．2关键点

关键点这里表示样条的最大点和最小点以及3次或以上样条的拐点（POI，参见图5．1)。对于最小点和最大点，曲线的斜率为零，因为线是水平的。曲线y=f(x)的这两种点都可以通过计算f′(x)=0的根来找到。



图5．1关键点


令r为方程f′(x)=0的根。要确定根对应于最小点还是最大点，请考虑r左右两侧的小位移δ。对于最小点，左侧(LHS)斜率应为负，右侧(RHS)斜率应为正。这意味着： 
f′(r-δ)<0和f′(r+δ)>0(5．1)
或者，由于斜率从负变为正，则斜率变化率为正，即
f″(r)>0(5．2)
对于最大点，LHS斜率应为正，RHS斜率应为负。这意味着： 
f′(r-δ)>0和f′(r+δ)<0(5．3)
或者，由于斜率从正变为负，则斜率变化率为负，即
f″(r)<0(5．4)
最小点或最大点的坐标由［r，f(r)］给出。





拐点是曲率从正(凸)变为负(凹)的位置，反之亦然。要满足的必要条件是在拐点处的曲率应为零，并且两侧的曲率符号应相反。令r为方程f″(x)=0的根。考虑r左右两侧的小位移δ。拐点存在的必要条件是
f″(x)=0
sign{f″(r-δ)}≠sign{f″(r+δ)}(5．5)
例5．1找到三次曲线的关键点： y=2+13x-31x2+18x3。

解： 

从给定的方程： 
f(x)=2+13x-31x2+18x3f′(x)=13x-62x+54x2f″(x)=-62+108x
设f(x)=13-62x+54x2=0的根为r1=0．87、r2=0．28。令δ=0．1。

现在，f(r1-δ)=f(0．87-0．1)=13-62(0．77)+54(0．77)2=-2．72<0。

此外，f(r1+δ)=f(0．87+0．1)=13-62(0．97)+54(0．97)2=+3．67>0。

因此，在x=0．87处有一个最小点。

由于f(0．87)=1．7，最小点的坐标为(0．87，1．7)。

验证： f(0．87)=-62+108(0．87)=-62+93．96=31．96>0。

另外，f(r2-δ)=f(0．28-0．1)=13-62(0．18)+54(0．18)2=+3．59>0。

并且，f(r2+δ)=f(0．28+0．1)=13-62(0．38)+54(0．38)2=-2．76<0。

因此，在x=0．28处有一个最大点。

由于f(0．28)=3．6，最大点坐标为(0．28，3．6)。

验证： f(0．28)=-62+108(0．28)=-62+30．24=-31．76<0。

对于拐点，设f(x)=-62+108x=0，根r=0．574。

现在sign{f(r-δ)}=sign{f(0．574-0．1)}=sign{-62+108(0．474)}=sign(-10．8)=-1。

并且sign{f(r+δ)}=sign{f(0．574+0．1)}=sign{-62+108(0．674)}=sign(10．8)=+1。

由于符号相反，在r=0．574处有一个拐点。

拐点的坐标： ［0．574，f(0．574)］即(0．574，2．662)(见图5．2)。



图5．2例5．1的绘图


MATLAB Code 5．1

clear all; clc;



%ax^3 + bx^2 + cx + d

p = ［18, -31, 13, 2］;

dp = polyder(p);

d2p = polyder(dp);

r = roots(dp);

if polyval(d2p, r(1)) > 0, m1 = 0; else m1 = 1; end;

b1 = polyval(p, r(1));

if m1 = = 0, fprintf('minimum point : (%.2f, %.2f)＼n', r(1), b1); end;

if m1 = = 1, fprintf('maximum point : (%.2f, %.2f)＼n', r(1), b1); end;

if polyval(d2p, r(2)) > 0, m2 = 0; else m2 = 1; end;

b2 = polyval(p, r(2));

if m2 = = 0, fprintf('minimum point : (%.2f, %.2f)＼n', r(2), b2); end;

if m2 = = 1, fprintf('maximum point : (%.2f, %.2f)＼n', r(2), b2); end;

s = roots(d2p);

s1 = polyval(d2p, s - 0.1);

s2 = polyval(d2p, s + 0.1);

b3 = polyval(p, s);

if (sign(s1) = = sign(s2)), fprintf('no poi');

else fprintf('poi : (%.2f, %.2f)＼n', s, b3); end;



%plotting

syms x;

y = p(1)*x^3 + p(2)*x^2 + p(3)*x + p(4);

X = ［r(1), r(2), s］; Y = ［b1, b2, b3］;

xx = linspace(min(［r(1), r(2), s］), max(［r(1), r(2), s］));

yy = subs(y, x, xx);

plot(xx, yy, 'b-', 'LineWidth', 1.5); hold on;

scatter(X, Y, 20, 'r', 'filled');

plot(X, Y, 'ko');

grid; xlabel('x'); ylabel('y');



%labeling

if m1 = = 0, text(r(1), b1+0.1, 'Minimum');

else text(r(1), b1+0.1, 'Maximum');

end;

if m2 = = 0, text(r(2), b2+0.1, 'Minimum');

else text(r(2), b2+0.1, 'Maximum');

end;

text(s, b3+0.1, 'POI');

hold off;

注解

polyder： 微分多项式。

roots： 找到多项式方程的根。

polyval： 在指定值处计算多项式。

sign： 返回参数+1、0或-1的符号。



5．3切线和法线

参数曲线可以表示为有序对，即
C(t)={x(t),y(t)}(5．6)
切向量由一阶导数获得： 
T(t)=C′(t)={x′(t),y′(t)}(5．7)
单位切向量通过除以其大小获得： 
T(t)/|T(t)|=C′(t)={x′(t),y′(t)}/|x′(t),y′(t)|(5．8)
法向量通过将切向量旋转90°得到： 
N(t)=R(90°)T(t)(5．9)
展开： 
N(t)=0-10100001x′(t)y′(t)1=-y′(t)x′(t)1(5．10)
写成有序对： 
N(t)={-y′(t),x′(t)}(5．11)
单位法向量通过除以其大小获得： 
N(t)/|N(t)|={-y′(t),x′(t)}/|{-y′(t),x′(t)}|(5．12)
通过曲线上一点(x1，y1)的切线方程： 
y-y1x-x1=y′(t)x′(t)(5．13)
通过曲线上一点(x1，y1)的法线方程： 
y-y1x-x1=-x′(t)y′(t)(5．14)
例5．2对于圆C(t)={cost，sint}，求点P(1/2，1/2)处的单位切向量、单位法向量、切线和法线。

解： 

给定曲线： C(t)={cost，sint}，

切向量： T(t)=C′(t)={x′(t)，y′(t)}={-sint，cost}=单位切向量； 

法向量： R(90°)C′(t)={-y′(t)，x′(t)}={-cost，-sint}=单位法向量。

现在在点P(1/2，1/2)处，求解cost=1/2，必须有t=π/4，

P处的单位切向量： {-sin(π/4)，cos(π/4)}={-1/2，1/2}； 

P处的单位法向量： {-cos(π/4)，-sin(π/4)}={-1/2，-1/2}。

P处的切线： y-y1x-x1=y′(t)x′(t)，

替换值： y-1/2x-1/2=cos(π/4)-sin(π/4)=-1，

化简： x+y-2=0。

P处的法线： y-y1x-x1=-x′(t)y′(t)，

替换值： y-1/2x-1/2=sin(π/4)cos(π/4)=1，



图5．3例5．2的绘图

化简： x-y=0。

例5.2的绘图如图5．3所示。

MATLAB Code 5．2

clear all; clc; format compact;

P = ［1/sqrt(2), 1/sqrt(2)］;

syms t;

x = cos(t); y = sin(t);

C = ［x, y］;

fprintf('Tangent vector ＼n');

T = ［diff(x), diff(y)］

fprintf('Normal vector ＼n');

N = ［-diff(y), diff(x)］

R = solve(sin(t)== 1/sqrt(2), cos(t)== 
1/sqrt(2));

fprintf('Tangent vector at P ＼n');

TP = subs(T, 't', R)

fprintf('Normal vector at P ＼n');

NP = subs(N, 't', R)

syms X Y;

fprintf('Tangent Line at P＼n');

Y1 = (TP(2)/TP(1))*(X - P(1))+P(2)

fprintf('Normal Line at P＼n');

Y2 = (-TP(1)/TP(2))*(X - P(1))+P(2)



%plotting

ezplot(x, y);hold on;

axis(［-2 2 -2 2］);

xx = linspace(-2,2);

yy1 = subs(Y1, X, xx);

yy2 = subs(Y2, X, xx);

plot(xx, yy1, 'r-', xx, yy2, 'r-');

scatter(P(1), P(2), 20, 'r', 'filled');

grid; hold off;

注解

diff： 计算导数和偏导数。

solve： 生成方程的解。



若曲线方程以隐式形式f(x，y)=0表示，则曲线方程在P(x0，y0)处的切线方程可通过首先计算在给定点的偏导数来进行，它为我们提供了该点的法向量:
N(x,y)=fx,fy(5．15)
因此，点P的法线由下式给出： 
Np=fpx,fpy(5．16)
P处的法线方程计算为
Tp： =fpx(x-x0)+fpy(y-y0)(5．17)
例5．3在点P(1，2)处找到曲线x3+2xy+y2=9的法向量和切线。

解： 

这里，f=x3+2xy+y2-9，

偏导数： fx(x，y)=f/x=3x2+2y和fy(x，y)=f/y=2x+2y； 

法向量： N(x，y)=［fx(x，y)，fy(x，y)］=［3x2+2y，2x+2y］； 

P点的法向量： N(1，2)=［7，6］； 

点P处的切线方程： 7(x-1)+6(y-2)=0，即7x+6y-19=0。

例5.3的绘图如图5.4所示。



图5．4例5．3的绘图


MATLAB Code 5．3

clear all; clc;

syms x y;

f = x^3 + 2*x*y + y^2 - 9;

df = ［diff(f, x), diff(f, y)］;

p = ［1, 2］;

fprintf('Normal vector : ＼n');

n = subs(df, ［x, y］, ［p(1), p(2)］)

fprintf('Tangent line : ＼n');

t = dot(n, ［x-1, y-2］)



% plotting

ezplot(f); hold on; grid;

ezplot(t);

quiver(p(1), p(2), n(1), n(2));

scatter(p(1), p(2), 20, 'r', 'filled');

hold off;

注解

dot： 计算向量点积。

quiver： 将向量描述为带有方向和大小的箭头。



5．4曲线长度

空间域中点x=a和x=b之间的曲线y=f(x)的长度由下式给出： 
L=∫ba1+dydx2dx(5．18)
这个表达式的推导假设一条曲线可以用许多小线段来近似(见图5．5)。对于如此小的线段，若δx是水平距离，δy是垂直距离，则线段的长度可以近似为：
δr=(δx)2+(δy)2(5．19)


图5．5曲线长度的推导


整条曲线的长度是x=a和x=b之间这些小段距离的总和。积分： 
L=∫ba(dx)2+(dy)2(5．20)
化简： 
L=∫badxdx2+dydx2dx=∫ba1+dydx2dx(5．21)

例5．4求x=1和x=5之间曲线y=3x+2的长度。

解： 

这里，dy/dx=3，根据式(5．21)，可以得到
L=∫ba1+dydx2dx=∫ba1+(3)2dx=∫ba10dx=410=12．65
验证： y(1)=5、y(5)=17、P1=(1，5)、P2=(5，17)、距离(P1，P2)=16+144=12．65。

MATLAB Code 5．4

clear all; clc;

syms x;

y = 3*x+2;

d = diff(y);

e = sqrt(1 + d^2);

f = int(e, 1, 5);

fprintf('Length of curve : %f＼n', eval(f));

注解

int： 整数符号表达。


若曲线方程以参数形式表示，即x(t)和y(t)，则由式(5．19)得出： 
dr=(dx)2+(dy)2
对t的微分： 
drdt=dxdt2+dydt2
积分： 
L=∫badxdt2+dydt2dt(5.22)
例5．5对0≤t≤2π，确定参数曲线x=3sint，y=3cost的长度。

解： 

这里，dx/dt=3cost，dy/dt=-3sint。

根据式(5．22)，
L=∫badxdt2+dydt2dt=∫2π0(3cost)2+(-3sint)2dt=∫2π03dt=6π
验证： 由于这代表一个半径为3的圆，所以圆的周长=2π×3=6π。

MATLAB Code 5．5

clear all; clc;

syms t;

x = 3*sin(t);

y = 3*cos(t);

dx = diff(x);

dy = diff(y);

r = sqrt((dx)^2 + (dy)^2);

f = int(r, 0, 2*pi);

fprintf('Length : %f＼n', eval(f));

5．5曲线下面积




图5．6曲线下面积

x=a和x=b之间的曲线y=f(x)下方的面积是通过考虑一个宽度为dx且高度为y=f(x)的非常小的矩形区域来计算的(见图5．6)。这个矩形区域的面积是f(x)dx。

为了找到整个区域，需要将小矩形区域在指定的范围内积分。
A=∫baf(x)dx(5．23)
例5．6求曲线y=10-x2在值x=-1和x=2之间的面积。

解： 

从式(5．23)可以得到(见图5．7)： 
A=∫baf(x)dx=∫2-1(10-x2)dx=27


图5．7例5．6的绘图


MATLAB Code 5．6

clear; clc;

syms x;

y = 10 - x^2;

f = int(y, -1, 2);

eval(f);

fprintf('Area : %f＼n', eval(f));



%plotting

xx = linspace(-4, 4);

yy = subs(y, x, xx);

plot(xx, yy);

xlabel('x'); ylabel('y');

grid; hold on;

axis(［-5 5, 0 12］);

text(-3, 9, 'y = 10 - x^2');



%filling

x = linspace(-1, 2);

y1 = 10 - x.^2;

y2 = zeros(1,100);

X = ［x,fliplr(x)］;

Y = ［y1,fliplr(y2)］;

fill(X,Y,'g');

alpha(0.25);

hold off;

注解

zeros： 生成一个用零填充的矩阵。

fliplr： 左右方向翻转数组。

fill： 用颜色填充多边形。

alpha： 设置透明度值。


上面计算的面积实际上是曲线f(x)和x轴之间的面积。但是，若曲线本身位于x轴下方，则由上述公式计算得出的面积为负值。在这种情况下，我们需要取面积的绝对值。因此，最好在找到曲线下的区域之前绘制曲线。

例5．7求值x=-1和x=2之间的曲线y=x3下的面积。

解： 

从式(5．23)可以得到(见图5．8)： 
A=∫baf(x)dx=∫0-1x3dx+∫20x3dx=x440-1+x4420=14+164=174=4．25


图5．8例5．7的绘图


或者，在积分之前先取函数的绝对值： 
A=∫2-1|x3|dx=x4sign(x)42-1=164-1(-1)4=174
与不正确的、没考虑符号而直接取积分的结果进行比较： 
A=∫2-1x3dx=x442-1=164-14=154

MATLAB Code 5．7

clear; clc;

syms x;

y = x^3;

f1 = abs(int(y, -1, 0));

f2 = int(y, 0, 2);

f = f1 + f2;

fprintf('Area : %f＼n', eval(f));



%alternative

f3 = int(abs(y), -1, 2);

fprintf('Alternatively : %f＼n', eval(f3));



%plotting

xx = linspace(-3, 3);

yy = subs(y, x, xx);

plot(xx, yy);

xlabel('x'); ylabel('y');

hold on;

X = ［-3 3］; Y = ［0, 0］; plot(X, Y, 'k--');

X = ［0 0］; Y = ［-10, 10］; plot(X, Y, 'k--');

axis(［-3 3, -10 10］); grid;

text(0.5, 4, 'y = x^3');

x = linspace(-1, 2);

y1 = x.^3;

y2 = zeros(1,100);

X = ［x,fliplr(x)］;

Y = ［y1,fliplr(y2)］;

fill(X,Y,'g'); alpha(0.25);

hold off;

若需要曲线和y轴之间的面积，则通过扩展先前的想法，我们首先以x=f(y)的形式表示曲线，然后通过绘制一系列宽度为dy和高度x=f(y)的非常小的矩形区域来计算，我们在指定的限制范围c~d上积分。在这种情况下,
A=∫dcf(y)dy(5．24)
例5．8求值y=3和y=5之间的曲线y=x+1下的面积。

解： 

重写x=f(y)=y2-1。

从式(5．24)可以得到(见图5．9)： 
A=∫dcf(y)dy=∫53(y2-1)dy=30．67


图5．9例5．8的绘图



MATLAB Code 5．8

clear all; clc;

syms y;

x = y^2 - 1;

f = int(x, 3, 5);

fprintf('Area : %f＼n', eval(f));



%plotting

yy = linspace(0, 6);

xx = subs(x, y, yy);

plot(yy, xx);

xlabel('y'); ylabel('x');

hold on;

X = ［3, 3］; Y = ［0, subs(x, y, 3)］; plot(X, Y, 'r');

X = ［5, 5］; Y = ［0, subs(x, y, 5)］; plot(X, Y, 'r');

plot(［0, 6］, ［0, 0］, 'k-'); grid;

x = linspace(3, 5);

y1 = x.^2 - 1;

y2 = zeros(1,100);

X = ［x,fliplr(x)］;

Y = ［y1,fliplr(y2)］;

fill(X,Y,'g'); alpha(0.25);

text(3, 20, 'x = y^2 - 1');

hold off;

由两条曲线f(x)和g(x)界定的面积是上曲线与x轴之间的面积减去下曲线与x轴之间的面积： 
A=∫ba|f(y)-g(x)|dx(5．25)
例5．9在值x=-1和x=1之间找到曲线y=x3和y=x所包围的区域面积。

解： 

根据式(5．25)，可以得到(见图5．10)：
A=∫ba|f(y)-g(x)|dx=∫1-1|x3-x|dx=0．5


图5．10例5．9的绘图



MATLAB Code 5．9

clear; clc;

syms x;

y1 = x;

y2 = x^3;

f = int(abs(y1 - y2), -1, 1);

fprintf('Area : %f＼n', eval(f));



%plotting

xx = linspace(-2,2);

yy1 = subs(y1, x, xx);

yy2 = subs(y2, x, xx);

plot(xx, yy1, xx, yy2);

xlabel('x'); ylabel('y');

axis(［-2 2 -2 2］); hold on; grid;

plot(［-2, 2］, ［0, 0］, 'k--');

plot(［0, 0］, ［-2, 2］, 'k--');

text(-1.5, -0.5, 'y = x^3');

text(-0.5, -0.75, 'y = x');

x = linspace(-1, 1);

y1 = x;

y2 = x.^3;

X = ［x,fliplr(x)］;

Y = ［y1,fliplr(y2)］;

fill(X,Y,'y'); alpha(0.25);

hold off;

5．6质心

地球对实体物体的每个粒子都施加引力。若所有这些力都被一个等效的力代替，则这个力将通过一个称为重心的点起作用。假设物体的密度是均匀的，重心将与质心重合。若密度均匀的物体是一块薄板，则质心将与该区域的质心重合。因此，板的质心是一个点，通过它可以平衡板的整个重量。很明显，对于矩形或圆形，质心正好在中心。在找到任意形状的质心之前，让我们首先考虑一些随机维度的多边形，如图5．11所示。



图5．11推导多边形的质心


为了找到多边形的质心，将它分成两个矩形，每个矩形的质心都在它的中心。左侧矩形的宽度为4，高度为4，因此其面积为A1=4 × 4=16，其中心位于(x1，y1)=(-3，1)。右侧矩形的宽度为5，高度为6，因此其面积为A2=5 × 6=30，其中心位于(x2，y2)=(1．5，2)。为了找到多边形的质心，我们相对于X轴和Y轴计算各个组成矩形的矩，并使它们的总和等于整个区域的矩。令x和y表示质心的坐标，A是它的整个区域面积。通过使关于X轴和Y轴的矩相等，可以得到： 
A1x1+A2x2=Ax-A1y1+A2y2=Ay-
代入数值，整个区域的质心由下式给出： 
x-=A1x1+A2x2A=16(-3)+30(1．5)16+30=-0．065
y-=A1y1+A2y2A=16(1)+30(2)16+30=1．652


图5．12导出曲线下区域的质心


综上所述，我们可以制定一个一般规则，即一个区域的质心等于其组成部分的矩之和除以整个区域面积。

图5．12描绘了函数y=f(x)的图形，需要找到由线x=a和x=b界定的曲线下方区域的质心。考虑一条宽度为dx的非常薄的条带，其与Y轴距离为x，高度为y=f(x)，因此条带的面积为ydx。因此，条带沿x方向的力矩为(ydx)x，条带沿y方向的力矩为(ydx)(y/2)，因为条带的中心正好位于沿条带高度的中间，因此，使区域内所有此类条带的矩相等并除以区域面积本身，我们得到区域的质心，如下所示： 
x-=∫baxf(x)dx∫baf(x)dx
y-=∫ba{f(x)}2dx2∫baf(x)dx(5．26)
例5．10求曲线y=x3、x=0和x=2所围区域的质心。

解： 

根据式(5．23)，面积： ∫baf(x)dx=∫20x3dx=4。

根据式(5．26)，质心坐标： 
x-=∫20xx3dx4=1．6
y-=∫20x6dx2×4=2．29
MATLAB Code 5．10

clear all; clc;

syms x y;

y = x^3;

a = int(y, 0, 2);

eval(a);

fprintf('Area : %f＼n', eval(a));

m1 = int(x*y, 0, 2);

m2 = int(0.5*y^2, 0, 2);

xc = m1/a;

yc = m2/a;

fprintf('Centroid : (%.2f, %.2f) ＼n', eval(xc), eval(yc));

若一个区域由区间［a，b］上的两条曲线f(x)和g(x)界定，则有界区域的质心如下所示： 
x-=∫bax［f(x)-g(x)］dx∫ba［f(x)-g(x)］dx
y-=∫ba{［f(x)］2-［g(x)］2}dx2∫ba［f(x)-g(x)］dx(5．27)
5．7插值和曲线拟合

为了结束关于样条的话题，我们最后看一下MATLAB提供的几个函数，这些函数用于执行与样条相关的两个重要任务，即插值和曲线拟合。这些任务仅专门使用编程工具完成，因此不会出现可以手动计算相关数值的问题。

当存在一些粗略的数据点并且我们有兴趣在点之间找到一些中间值时可以进行插值，但这些值不是直接提供的［Mathews，2004］。因此，可以对给定的数据进行插值以找到这些值。MATLAB提供的第一个函数是interp1，它代表一维插值，默认情况下在给定数据之间使用线性插值。它还有许多其他选项，如图5．13所示，其中使用了两种类型的数据模式，分别是阶跃函数和正弦函数。第一行表示使用“线性”选项进行插值，第二行表示使用“最近邻居”选项进行插值，第三行表示使用“前一个邻居”选项进行插值，第四行表示使用“下一个邻居”选项进行插值。



图5．13使用带有选项的interp1进行插值


生成显示图的相应代码如下： 

MATLAB Code 5．11

clear; clc;

x = -3:3;

y = ［-1 -1 -1 0 1 1 1］;

t = -3:.01:3;

subplot (221)

plot(x,y,'o',t,interp1(x,y,t, 'linear'), 'r-', 'LineWidth', 2); title('linear');

subplot (223)

plot(x,y,'o',t,interp1(x,y,t, 'nearest'), 'r-', 'LineWidth', 2); title('nearest');

x = 0:2*pi;

y = sin(x);

t = 0:.01:2*pi;

subplot (222)

plot(x,y,'o',t,interp1(x,y,t, 'linear'), 'r-', 'LineWidth', 2); title('linear');

subplot (224)

plot(x,y,'o',t,interp1(x,y,t, 'nearest'), 'r-', 'LineWidth', 2); title('nearest');

figure

x = -3:3;

y = ［-1 -1 -1 0 1 1 1］;

t = -3:.01:3;

subplot (221)

plot(x,y,'o',t,interp1(x,y,t, 'previous'), 'r-', 'LineWidth', 2); title('previous');

subplot (223)

plot(x,y,'o',t,interp1(x,y,t, 'next'), 'r-', 'LineWidth', 2);title('next');

x = 0:2*pi;

y = sin(x);

t = 0:.01:2*pi;

subplot (222)

plot(x,y,'o',t,interp1(x,y,t, 'previous'), 'r-', 'LineWidth', 2); title('previous');

subplot (224)

plot(x,y,'o',t,interp1(x,y,t, 'next'), 'r-', 'LineWidth', 2);title('next');

注解

interp1： 执行一维插值。


第二个函数是pchip，代表“分段三次厄米特插值多项式”。它提供了数值的保形分段三次厄米特插值的分段多项式形式。对于每个子区间，它在端点之间进行插值，并且还保持端点处的斜率是连续的。如图5．14所示，图(a)是阶跃函数，图(b)是正弦函数。



图5．14使用pchip进行插值


生成绘图的相应代码如下： 

MATLAB Code 5．12

clear; clc;

x = -3:3; y = ［-1 -1 -1 0 1 1 1］; t = -3:.01:3;

subplot (121)

plot(x,y,'o',t, pchip(x,y,t), 'r-', 'LineWidth', 2);

axis tight; axis square;

x = 0:2*pi; y = sin(x); t = 0:.01:2*pi;

subplot (122)

plot(x,y,'o',t, pchip(x,y,t), 'r-', 'LineWidth', 2);

axis tight; axis square;

第三个函数是spline，代表“分段三次样条插值多项式”。它在数据值上提供了三次样条的分段多项式形式。阶跃函数和正弦函数如图5．15所示。



图5．15使用spline插值


生成绘图的相应代码如下： 

MATLAB Code 5．13

clear; clc;

x = -3:3; y = ［-1 -1 -1 0 1 1 1］; t = -3:.01:3;

subplot (121)

plot(x,y,'o',t, spline(x,y,t), 'r-', 'LineWidth', 2);

axis tight; axis square;

x = 0:2*pi; y = sin(x); t = 0:.01:2*pi;

subplot (122)

plot(x,y,'o',t,spline(x,y,t), 'r-', 'LineWidth', 2);

axis tight; axis square;

图5．16提供了三者之间的比较。interp1函数通过直线连接数据点，spline函数提供最平滑的曲线，而pchip函数能在端点减少振荡。



图5．16插值方法之间的比较


为了将曲线拟合到一组给定的数据点，MATLAB提供了一个名为polyfit的函数，该函数创建一个指定次数的拟合多项式，其与数据点的平方误差之和最小，并返回多项式的系数。下面的代码说明了这个过程。其中，生成了两组数据点，并使用了不同阶的曲线来拟合数据。通过改变d的值绘制结果。多项式的系数由函数返回。图5．17上一行显示阶跃函数，下一行显示正弦函数。各列对应用于曲线拟合的多项式的次数，即1、3和9。



图5．17使用polyfit的曲线拟合


生成绘图的相应代码如下： 

clear; clc;

x = -10:10;

y = ［-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 1 1 1 1 1 1 1 1 1 1］;

t = -10:.01:10;

subplot (231)

d = 1; pf = polyfit(x, y, d);

pv = polyval(pf, t);

plot(x,y,'o',t, pv, 'LineWidth', 2);

axis tight; axis square;

title('degree = 1');

subplot (232)

d = 3; pf = polyfit(x, y, d);

pv = polyval(pf, t);

plot(x,y,'o',t, pv, 'LineWidth', 2);

axis tight; axis square;

title('degree = 3');

subplot (233)

d = 9; pf = polyfit(x, y, d);

pv = polyval(pf, t);

plot(x,y,'o',t, pv, 'LineWidth', 2);

axis tight; axis square;

title('degree = 9');

clear x y t pf pv;

x = -pi:pi;

y = sin(x);

t = -pi:.1:pi;

subplot (234)

d = 1; pf = polyfit(x, y, d);

pv = polyval(pf, t);

plot(x,y,'o',t, pv, 'LineWidth', 2);

axis tight; axis square;

subplot (235)

d = 3; pf = polyfit(x, y, d);

pv = polyval(pf, t);

plot(x,y,'o',t, pv, 'LineWidth', 2);

axis tight; axis square;

subplot (236)

d = 9; pf = polyfit(x, y, d);

pv = polyval(pf, t);

plot(x,y,'o',t, pv, 'LineWidth', 2);

axis tight; axis square;

axis(［-pi pi -1 1］);

注解

polyfit： 生成多项式以拟合给定数据。



5．8关于二维绘图函数的说明

在我们将主要关注点从二维域转移到三维域之前，本节总结了所使用的MATLAB 二维绘图函数和一些附加函数［Marchand，2002］。鼓励读者从MATLAB文档中探索有关这些函数的更多详细信息。

(1) ezplot： 此函数可用于使用符号变量进行绘图。

① 一个变量(见图5．18)： 

ezplot('1 - 2.25*t + 1.25*t^2');

② 两个变量(见图5．19)： 

ezplot('x^4 + y^3 = 2*x*y');




图5．18使用ezplot绘制一个变量




图5．19使用ezplot绘制两个变量



③ 参数变量(见图5．20)： 

ezplot('cos(t)', 'sin(t)');

(2) plot： 此函数可用于使用值向量进行绘图。

① 参数方程(见图5．21)： 

t = 0:pi/50:10*pi; plot(t.*sin(t),t.*cos(t));




图5．20使用ezplot绘制参数变量




图5．21使用plot绘制一个变量





图5．22使用plot绘制显式方程

② 显式方程(见图5．22)： 

x = -pi:.1:pi; y = tan(sin(x)) - sin(tan(x)); plot(x,y);

(3) ezcontour： 等高线图，一个版本只有边缘，另一个是填充版本(见图5．23)。

ezcontour ('x*exp(-x^2 - y^2)');

ezcontourf ('x*exp(-x^2 - y^2)');



图5．23使用ezcontour绘制的等高线图


(4) fimplicit： 绘制隐式函数(从MATLAB 2016版引入)。已安装MATLAB包的版本可以通过在命令行输入ver来检查(见图5．24)。

fimplicit(@(x,y) sin(x.^2 + y.^2) + cos(x.^2 + y.^2) - 1, ［-10 10］);




图5．24使用fimplicit绘制的隐式函数图

(5) patch： 根据顶点和颜色创建填充多边形(见图5．25)。

x = ［4 6 11 9］;

y = ［2 7 9 4］;

c = ［0 4 6 8］;

colormap(jet);

patch(x,y,c);

colorbar;

hold;

v1 = ［2 4; 2 10; 8 4］;

patch('Vertices', v1, 'FaceColor', 'red', 'FaceAlpha', 0.3, 'EdgeColor', 'red');

axis(［0 12 0 12］);

(6) fplot： 绘制连续或分段函数(见图5．26)。

fplot(@(x) sin(x),［-2*pi pi］,'b');

hold on;

fplot(@(x) cos(x),［-pi 2*pi］,'r');

hold off;












图5．25使用patch绘图




图5．26使用fplot绘图



注解

colormap： 使用预定义的颜色查找表指定颜色方案。

colorbar: 通过在颜色图中附加颜色来创建颜色条。



5．9本章小结

以下几点总结了本章讨论的主题： 

 样条曲线的关键点是最大点、最小点和拐点。

 对于最小点，斜率变化率为正。

 对于最大点，斜率变化率为负。

 对于拐点，曲率应为零且两侧应有相反的符号。

 曲线的切线由曲线方程的一阶导数获得。

 曲线的法线通过将切线旋转90°获得。

 曲线段的长度是通过将非常小的段的长度相加得到的。

 曲线下的面积是通过将非常小的矩形的面积相加而获得的。

 面积在x轴上方为正，在x轴下方为负。

 由两条曲线界定的面积是两条曲线下面积的绝对差。

 两个数据点之间的插值可以使用许多选项线性地完成。

 三次样条插值可以提供更平滑的插值曲线。

 可以使用特定次数的多项式对一组数据点进行曲线拟合。

 MATLAB为隐式、显式和参数表达式提供了许多绘图函数。

5．10复习题

1． 样条曲线的关键点是什么意思？

2． 如何用曲线的梯度区分最小点和最大点？

3． 曲线的拐点是如何确定的？为什么三次曲线总是有一个拐点？

4． 用有序对表示参数曲线，切线和法线如何计算？

5． 如何确定两个给定点之间的曲线长度的表达式？

6． 一条曲线与两条垂直线所包围的x轴之间的面积是如何确定的？

7． 计算x轴上方和下方的区域面积要遵循什么约定？

8． 两条给定曲线之间的面积是如何确定的？

9． 观察到的线性插值和三次样条插值有什么区别？

10． 为了获得一组数据点的最佳拟合曲线，需要满足哪些标准？

5．11练习题

1． 求y=x1．5在x=0和x=5之间的长度。

2． 对0≤x≤π/4，确定曲线y=ln(secx)的长度。

3． 从x=0到x=1找到y=x2和x=y2之间的区域。

4． 求曲线x=y+1和x=0．5y2-3之间的完全封闭区域。

5． 找出曲线的最小点、最大点和拐点： y=-2．67x+0．67x3。

6． 求曲线(t，t2)在点(1，1)处的切向量、法向量和切线方程。

7． 求参数曲线x=t3、y=t2在(-8，4)处的单位切向量和单位法向量。

8． 对于隐式曲线x3+y2-x=4，求点(-1，2)处的切线和法向量。

9． 对于摆线(t+sint，1-cost)，求t=π/2处的切向量和法向量。

10． 求曲线y=x和y=x3所围区域的质心。