飞机机翼图设计

伯恩思坦多项式与Bezier曲线
一、引言
1971年法国雷诺汽车公司的工程师Bezier提出了一种新的参数曲线表示法。

这种方法能方便地控制输入参数(控制点)以改变曲线的形状。

被称为Bezier曲线，数学原理使用了伯恩思坦多项式。

设f(x)是定义在[0，1]上的连续函数，称表达式
∑=
--
≈
n
k
k
n
k
k
n
t
t
C
n
k
f
x
f
)
1(
)
/
(
)
(
右端为函数的伯恩思坦逼近多项式。

下面是函数)
sin(
)
(x
x
fπ
=的伯恩思坦多项式逼近实验程序
n=input('input n=');
x=[0:n]/n; f=sin(x*pi);
for i=1:n+1
y=f;t=x(i);
for k=n:-1:1
for j=1:k
y(j)=t*y(j)+(1-t)*y(j+1);
end
end
p(i)=y(1);
end
max(abs(f-p))
plot(x,f,'b',x,p,'o',x,p,'r')
下面两图分别是取不同点数的伯恩思坦多项式逼近。

n=10逼近n=20逼近
二、Bezier曲线
Bezier曲线的形状是通过一组多边折线(控制多边形)的各顶点P0，P1，…，P m所定义出来的。

在多边折线的各顶点中，只有第一点P0和最后一点P m在曲线上，其余的点则用以定义曲线的阶次。

给定控制多边形顶点P0，P1，…，P m的坐标
(x0，y0)，(x1，y1)，……，(x m，y m)
曲线参数方程为
∑=
--
=
m
k
k
k
m
k
k
m
x
t
t
C
t
x
) 1(
)
(，∑
=
--
=
m
k
k
k
m
k
k
m
y
t
t
C
t
y
) 1(
) (
其中，k
m C 为组合数，其计算公式为
!
)(!!k m k m C k
m -=
在工程实际中，人们常用矢量函数的形式来表示平面曲线，若记
⎥⎦⎤
⎢⎣⎡=)()()(t y t x t P ， ⎥⎦
⎤⎢⎣⎡=k k k
y x P (k = 0，1，…，m )
则有数学表达式
∑=--=
m
k k k
m k k m
P t t C
t P 0
)
1()(
作为特殊情形，下面分别给出 m = 1，2，3时的Bezier 曲线数学表达式。

1、一次Bezier 曲线是通过平面上两点P x y P x y 000111(,),(,)的直线段，其数学表达式为
10)1()(P t P t t P +-=，（10≤≤t ）
2、 二次Bezier 曲线是由平面上三个点P x y P x y P x y 000111222(,),(,),(,)组成控制多边形，并由此确定的抛物线。

数学表达式为 P t t P t t P t P t ()()()=-+-+≤≤121012
012
2
3、三次Bezier 曲线是由平面上四个点P x y P x y P x y P x y 000111222333(,),(,),(,),(,)组成控制多边形，并由此所确定的三次曲线。

数学表达式为 P t t P t t P t t P t P t ()()()()=-+-+-+≤≤13131013
02
12
23
3
绘制曲线的算法常用加权平均法。

这一算法的根据是组合数计算中常用的递推公式
1
11---+=k m k
m k
m C C C
以二次Bezier 曲线为例说明原理，由于
])1[(])1)[(1()1(2)1()(211022
102
tP P t t tP P t t P t P t t P t t P +-++--=+-+-= 由于 (1－t ) 和 t 都是介于 0 和 1 之间的数，所以认为是加权平均的权值。

给定三个
点P 0(x 0，y 0)、P 1(x 1，y 1)、P 2(x 2，y 2)作为控制多边形顶点后，二次Bezier 曲线上对应于参数 t 的点P (t )的坐标可用如下方法计算 （1）取P P 01,点坐标的加权平均得P 01的坐标 1001)1(tx x t x +-=
1001)1(ty y t y +-=
（2）取P P 12,点的坐标的加权平均计算P 12的坐标 2112)1(tx x t x +-= 2112)1(ty y t y +-=
（3）取P P 0112,点的坐标的加权平均计算P (t )的坐标 1201)1()(tx x t t x +-= 1201)1()(ty y t t y +-=
曲线的矩阵表示算法
由于控制多边形各顶点的坐标可表示为矩阵形式的数据结构，借用矩阵的表示也可以实现算法。

仍以计算二次Bezier 曲线上对应于t t t 01100011===,.,, 的11个点的坐标数据为例说明算法原理。

由于二次Bezier 曲线的数学表达式可以写作矩阵形式
⎥⎥⎥⎦
⎤
⎢⎢⎢
⎣⎡--=21022
])1(2)1[()(P P P t t
t t t P 我们希望最后所得的曲线上11个点的坐标数据以11×2的矩阵形式给出，以第一列表
示X 坐标，以第二列表示Y 的坐标。

由上式
⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡22
1100210y x y x y x P P P 是一个3×2矩阵(它是算法的初始输入数据)，而])1(2)1[(2
2
t t t t --应该是11×3矩阵。

三、应用
例1、飞机机翼图设计。

利用贝塞尔曲线设计机翼剖面图，提取曲线数据，绘制机翼柱面图
MATLAB 程序：
p=[0 0;0.01 0.5;0.5 1;2 0]; %输入控制多边形顶点
图4-5
n=10;
t=(0:n)'/n;t1=1-t;
z=[t1.^3 3*t.*t1.^2 3*t1.*t.^2 t.^3]*p; %矩阵法计算曲线坐标数据px=p(:,1);py=p(:,2); %提取控制多边形顶点坐标x=z(:,1);y=z(:,2); %提取曲线坐标
figure(1),plot(px,py,px,py,'o',x,y) %绘机翼剖面图
x=[x(n+1:-1:2);x]; %利用对称性扩充曲线坐标y=[y(n+1:-1:2);-y];
E=ones(1,11);X=x*E;Y=y*E; %制做柱面坐标数据
Z=ones(2*n+1,1)*[0:10];
figure(2),mesh(Z,X,Y),hold on %绘柱面
colormap([0 0 1])
axis off
fill3(zeros(21,1),x,y,'c') %填充机翼剖面
view(-50,46)。