1 四种插值方法
上图中同时显示了4种插值方法的图像,分别是:
1 单段多项式(绿色),六次;
2 分段多项式(蓝色),每段都是三次的;
3 B样条曲线 (红色),三次的;
4 贝塞尔曲线(黑色),三次的;
2 特点
2.1 控制点
从图中可以看出,在这四种插值曲线中单段多项式和分段多项式都通过了控制点(红色点),而B样条曲线和贝塞尔曲线没有通过。
在控制点处曲线都是连续且光滑的。
在用鼠标改变一个控制点的时候,单段多项式和分段多项式插值曲线的所有片段都变化了,而B样条曲线和贝塞尔曲线只在控制点所在的片段有变化,其它的片段不变。
两个控制点很接近时,单段多项式出现了剧烈的震荡,而其它几种方法没有,这也是实际中基本不使用单段多项式插值的原因之一。
2.2 Mathematica程序
画图的程序如下:
DynamicModule[{pts={{0,0},{1,1},{2,-1},{3,0},{5,2},{6,-1},{7,3}}},
LocatorPane[Dynamic[pts],Dynamic[
p1=Plot[InterpolatingPolynomial[pts,x],{x,0,7},PlotStyle->Darker@Green,PlotRange->7];
f1=Interpolation[Transpose@{Range[1,7,1],pts[[;;,1]]},InterpolationOrder->3,Method->"Spline"];
f2=Interpolation[Transpose@{Range[1,7,1],pts[[;;,2]]},InterpolationOrder->3,Method->"Spline"];
Show[Graphics[{BezierCurve[pts,SplineDegree->3],Red,BSplineCurve[pts,SplineDegree->3],Blue,Line@Transpose@{f1@Range[1,7,0.01],f2@Range[1,7,0.01]},Red,PointSize[0.01],Point[pts]},PlotRange->{{-1,8},{-3,5}},ImageSize->600],p1]
],Appearance->None]]
2 B样条曲线
B样条曲线最早是为汽车、船舶等工业产品做外观造型设计而发明的,但是由于它有一些非常好的性质(光滑、能灵活地改变曲线形状、表达精确),所以被用到越来越多的问题上,比如求微分方程的数值解。我们都知道,积分比微分难,对于一些复杂的问题,只会用普通的数值积分(比如欧拉积分、龙格库塔积分法)是不够的。一般来说,想求解稍微复杂一点的微分方程问题时(例如IVP、BVP、偏微分方程),你发现简单的打靶法(shooting method)不好用,因为方程可能对初值非常敏感(例如敏感到需要精确到小数点后20位),所以经常采用配置法(collocation method)。要使用配置法就得选择样条曲线方程并需要计算曲线的导数。下面谈谈如何计算B样条曲线的导数。
2.1 基础知识
控制点决定了B样条曲线的形状,给定控制点和次数(degree),可以构造一条唯一的B样条曲线,例如由以下$6$个控制点生成的$3$次B样条曲线如下图所示。如果我们改变控制点,曲线的形状也随之改变。具体来说,这条曲线是一个二维的曲线(曲线上的点有 $x y$ 两个坐标),它也是个参数曲线,给定一个参数值可以得到一个二维点,把所有参数下的函数值画出来就是图中看到的黑色曲线。
controlPts = {{0, 0}, {1, 1}, {2, -1}, {3, 0}, {4, -2}, {5, 1}};
假如你想知道在某个参数坐标 $u$ 处的二维点坐标应该如何计算呢?最高效的计算方法不是直接用B样条曲线的定义式,而是采用Cox-de Boor公式,它是一个递归形式的公式,从最底层的控制点开始经过 $p$ 步构造出最终的点, $p$ 是曲线次数。$3$次B样条曲线(third-degree (cubic) B-spline)是最常用的B样条曲线,这可能是因为大多数工程和物理问题都是用二阶光滑微分方程描述的,选择$3$次样条足矣。$3$次B样条曲线的方程总是$3$次多项式。
2.2 B样条的导数
B样条曲线不仅连续而且光滑,所以存在导数,如何求导呢?B样条曲线的导数仍然是一个B样条曲线,只不过次数降低一次,控制点减少一个。所以仍然可以像以前计算B样条曲线一样计算B样条曲线的导数。而且,既然B样条曲线的导数还是B样条曲线,利用这个递归关系,更高阶的导数依然是B样条曲线。下面的例子展示了B样条曲线的一阶导数和二阶导数,你可以把它们看成速度和加速度。图中的红色箭头表示一阶导数,黑色箭头表示二阶导数,为了方便展示,这里只画出了导数的方向,实际不同点处的导数幅值(箭头长度)是不一样的。当然还可以很容易地继续画出更高阶的导数,这里就不展示了。
2.3 解微分方程
我们来用B样条解几个微分方程看看效果怎么样。
3.1 线性常微分方程BVP
先试试简单点的线性常微分方程,我们来解边界值问题(BVP),即给定两个端点的函数值。
$$
y''+y'-6y=x \ \ \ \ \ \ \ \ \ \ x\in[0,1]
\newline y(0)=0,\ \ y(1)=1
$$
这个方程存在解析解,如下$^{[1]}$
$$
y(x)= \frac{\left(43-e^2\right) e^{-3
x}-\left(43-\frac{1}{e^3}\right) e^{2 x}}{36
\left(\frac{1}{e^3}-e^2\right)}-\frac{x}{6}-\frac{1}{36}
$$
B样条的数值解如下图(黄色),作为对比,也画出了解析解(黑色),相对误差最大值为$3.72797\times10^{-7}$,求解中使用的参数为collocation point数量为101,控制点数量为100。可以看到数值求解的精度还是很高的。
3.2 非线性常微分方程BVP
再试试稍微复杂点的非线性常微分方程,同样考虑一个边界值问题(BVP),即被称为Bratu's problem的
$$
y''+ \lambda e^{y}=0 \ \ \ \ \ \ \ \ \ \ x\in[0,1]
\newline y(0)=0,\ \ y(1)=0
$$
这个方程的解析解$^{[2]}$如下
$$
y(x)=-2\text{ln}\left[\frac{\text{cosh}((x-\frac{1}{2})\frac{\theta}{2}
)}{\text{cosh}(\frac{\theta}{4})}\right]
$$
参数$\lambda=3.513830719$、$\theta=4.79871456$时的B样条的数值解如下图(黄色),作为对比,也画出了解析解(黑色),相对误差最大值为$2.77846\times10^{-3}$,求解中使用的参数为collocation point数量为101,控制点数量为100。
3 Clothoids曲线
如果你想用Clothoids曲线给机器人规划路径,它最大的缺点是没有封闭的表达形式(closed-form expression)。什么意义呢?就是说,如果你想得到曲线上的一个点,你得从头开始算积分,一直算到这个点,显然算积分太慢了,所以直接使用Clothoids曲线很难满足实时性的要求,一般只能用在离线的场合下,也就是预先规划好。
在Mathematica软件中画出Clothoids曲线的方法如下:
Show[ParametricPlot[{FresnelC[x],FresnelS[x]},{x,0,4}],Graphics[{Red,PointSize[0.01],Point/@Table[{FresnelC[x],FresnelS[x]},{x,0,4,0.1}]}]]
[1] Use of Cubic B-Spline in Approximating Solutions of Boundary Value Problems, Maria Munguia, 2015.
[2] On exact and numerical solutions of the one-dimensional planar Bratu problem, Ron Buckmire, 2003.