三次样条插值算法详解.ppt

1 三次样条插值 鉴于高次插值不收敛又不稳定的特点 低次插值既具有收敛性又具有稳定性 因此低次值更具有实用价值 但是低次插值的光滑性较差 比如分段线性插值多项式在插值区间中仅具有连续性 在插值节点处有棱角 一阶导数不存在 分段三次Hermite插值多项式在插值区间中仅具有一阶导数即一阶光滑性但不具备二阶光滑性 不能满足某些实际应用比如汽车 轮船 飞机等的外形中流线形设计 样条是在二十世纪初期经常用于图样设计的一种富有弹性的细长条 多个样条互相弯曲连接后沿其边缘画出的曲线就是三次样条曲线 后来数学上对其进行了抽象 定义了m次样条函数 并成为数值逼近的重要研究分枝 进一步扩大了样条函数的应用范围 2 样条函数的定义 定义4 1设区间 a b 上给定一个节点划分a x0 x1 xn 1 xn b如果存在正整数k使得 a b 上的分段函数s x 满足如下两条 1 在 a b 上有直到k 1阶连续导数 2 在每个小区间 xi xi 1 上是次数不大于k的多项式 则称分段函数s x 是以 2 6 为节点集的k次样条函数 3 三次样条插值函数的定义 并且关于这个节点集的三次样条函数s x 满足插值条件 则称这个三次样条函数s x 为三次样条插值函数 4 三次样条插值函数的边界条件 插值条件 连续性条件 一阶导数连续条件 二阶导数连续条件 5 1 因为s x 在每个小区间上是一个次小于三次的多项式 故有四个未知系数 2 因为s x 有n分段 从而共有4n个未知系数 3 但插值条件与样条条件仅给出4n 2个条件 无法定出4n个未知系数 还差2个条件 这2个条件我们用边界条件给出 6 通常我们对插值多项式在两端点的状态加以要求也就是所谓的边界条件 第一边界条件 由区间端点处的一阶导数给出即 7 第二边界条件 由区间端点处的二阶导数给出即 特殊情况为自然边界条件 由区间端点处的二阶导数恒为0给出即 8 这样三次样条插值问题就分成三类 其实不止这三类 第三类又称周期边界条件 由区间端点处的函数值或导数值满足周期条件给出 9 样条函数的例子 容易验证 是满足如下数据的第一类边界样条插值问题解 10 样条函数的例子 11 通常有三转角法 三弯矩法 B样条基函数法 三次样条插值函数的求法 这三种方法的基本思想是类似的 都是通过待定某些参数来确定插值函数 但肯定不是待定4n个参数 而是利用已知条件将待定参数减小到最少 比如 待定一阶导数 待定二阶导数 采用基函数方法来确定插值函数 12 三转角法 待定一阶数 为了确定三次样条插值函数的表达式S x 我们采用待定系数法来求解 我们待定什么系数呢 考虑到带一阶导数的分段三次Hermite插值多项式 13 我们采用待定一阶导数的方法即设 因为分段三次Hermite插值多项式已经至少是一阶连续可导了 为了让它成为三次样条函数只需确定节点处的一阶导数使这些节点处的二阶导数连续即可 14 15 由于在内部节点处二阶导数连续条件 整理化简后得 16 称为三转角法基本方程组 以上推导还没有考虑边界条件 针对不同类型的三次样条问题 就可以导出不同的方程组 17 第一类三次样条插值问题方程组 基本方程组化为n 1阶方程组 由于已知 化为矩阵形式 18 这是一个严格对角占优的三对角方程组 用追赶法可以求解 19 第二类三次样条插值问题的方程组 由于已知 故得 20 稍加整理得 联合基本方程组得一个n 1阶三对角方程组 化成矩阵形式为 仍然是严格对角占优 21 第三类样条插值问题的方程组 立即可得下式 由于 22 其中 联合基本方程得一个广义三对角或周期三对角方程组 这个方程组的系数矩阵仍然是严格对角占优阵 23 求解这些不同类型的样条插值问题的方程组 我们可得所要待定的一阶导数 称为三次样条插值问题三转角公式 再代入S x 的每一段表达式 就求得三次样条函数的表达式 24 例1 对于给定的节点及函数值 解 这是自然边界条件下的样条问题 25 我们可以将上述计算列于表中 26 由些得如下方程组 利用三转角公式 27 同样可以求得第三段表达式 28 29 30 三弯矩法 待定二阶导数 选择二阶导数作为待定参数 由于三次样条S x 是三次多项式 故它的二阶导数是一次多项式 从而 思考 1 的原因 31 32 从而推导出了三次样条S x 在第k个小区间 xk xk 1 上的表达式为 它的系数都是用二阶导数与函数值表示 33 对所有中间节点xk k 1 2 n 1 左边小区间与右边小区间上的三次多项式的一阶导数应当连续 确定二阶导数 34 三弯矩法基本方程 注意到这个基本方程只包括了n 1个方程 但却有n个二阶导数需要待定 这是一个欠定方程组 还需要根据边界条件再确定两个方程 35 曲率调整样条 这种样条的边界条件是已知两端点的二阶导数值 这样从三弯矩基本方程可以导数确定其它n 2个待定参数的方程组 36 自然样条 这种样条的边界条件是 已知两端点的二阶导数值为0 这样从三弯矩基本方程可以导数确定其它n 2个待定参数的方程组 37 固支样条 这种样条的边界条件是 已知两端点的一阶导数值 根据前面推导过程中得到的样条函数S x 的一阶导数的表达式 2 11 得方程 38 固支样条 这样从三弯矩基本方程可以导数确定n个待定参数的方程组 39 非扭结样条 这种样条的边界条件是 要求样条S x 在开始的两个小区间 x0 x1 x1 x2 上的三阶导数相同 在最后两个小区间 xn 2 xn 1 xn 1 xn 上的三阶导数相同 对表达式 2 9 再求一次导数得方程 40 非扭结样条 再由三弯矩基本方程 可得 41 周期样条 这种样条的边界条件是 要求样条S x 及其导数是以区间长度xn x0为周期的函数即 这些条件可以确定如下两个方程 42 再由三弯矩基本方程 可得 周期样条 43 在Matlab中数据点称之为断点 如果三次样条插值没有边界条件 最常用的方法 就是采用非扭结 not a knot 条件 这个条件强迫第1个和第2个三次多项式的三阶导数相等 对最后一个和倒数第2个三次多项式也做同样地处理 Matlab中三次样条插值也有现成的函数 y interp1 x0 y0 x spline y spline x0 y0 x pp csape x0 y0 conds pp csape x0 y0 conds valconds y ppval pp x 其中x0 y0是已知数据点 x是插值点 y是插值点的函数值 对于三次样条插值 我们提倡使用函数csape csape的返回值是pp形式 要求插值点的近似函数值 必须调用函数ppval MATLAB中三次样条函数法实现 44 pp csape x0 y0 conds valconds conds指定插值的边界条件 其值可为 complete 边界为一阶导数 一阶导数的值在valconds参数中给出 not a knot 非扭结条件 periodic 周期条件 second 边界为二阶导数 二阶导数的值在valconds参数中给出 若忽略valconds参数 二阶导数的缺省值为 0 0 MATLAB中三次样条函数法实现 45 例2 第一边界条件的例题 clear x 1 2 4 5 y 1 3 4 2 pp csape x y complete 17 8 19 8 pp coefs MATLAB代码 46 第一边界条件的例题 pp form pp breaks 1245 coefs 3x4double pieces 3order 4dim 1pp coefs 0 125002 12501 0000 0 1250 0 37501 75003 00000 3750 1 1250 1 25004 0000 分段多项式结构与系数矩阵 47 如下代码求解上述样条问题 例3 自然边界条件例题 48 x 4 3 2 101234 y 00 151 122 362 361 460 490 060 pp csape x y second xx 4 0 01 4 yy ppval pp xx holdon plot x y ok plot xx yy k holdoff MATLAB程序 49 pp form pp breaks 4 3 2 101234 coefs 8x4double pieces 8order 4dim 1 分段多项式结构 50 pp coefs 0 180856038291610 00000000000000 0 030856038291610 0 084280191458030 542568114874820 511712076583210 15000000000000 0 393735272459500 289727540500741 344007731958761 120000000000000 14922128129602 0 891478276877760 742256995581742 360000000000000 13685014727540 0 44381443298969 0 593035714285712 360000000000000 13337812960236 0 03326399116348 1 070114138438881 46000000000000 0 060362665684830 36687039764359 0 736507731958760 49000000000000 0 061927466863030 18578240058910 0 183854933726070 06000000000000 分段多项式的系数矩阵 51 图像 52 例4 对如下Runge现象中的函数 求用n分点作节点的三次样条插值多项式s3 x 的图象 取n 5 10 15 20等 将区间等分成n份 取分点作为插值节点 利用matlab函数csape 作出它的三次样条插值函数 并作出这些插值函数与被插函数的图像 但由于插值节点数目较多 故不列出插值函数的表达式 只画出它们的图象与被插函数的图象的复合 从这些图象可以看出 随着节点数目的增多 样条插值函数图象越来越符合被插函数 因此避免Runge现象 也就是说三次样条插值函数具有收敛性 53 54 55 5