第四章 MATLAB多项式与插值.ppt

第四章多项式与插值 4 1MATLAB与多项式 一 多项式的建立1 MATLAB中多项式用行向量表示 其元素为多项式的系数 且从左至右按降幂排列 2 已知一个多项式的全部根X 求多项式系数的函数是poly X 该函数返回以X为全部根的一个多项式P 首项系数为1 当X是一个长度为m的向量时 P是一个长度为m 1的向量 3 给定n 1个点可以唯一确定一个n阶多项式 利用polyfit可以确定多项式的系数 调用格式为 p polyfit x y n 其中x y是同维向量 代表数据点的横 纵坐标 n是多项式的阶数 二 多项式计算1 多项式求根求多项式p x 的根的函数是roots P 这里 P是p x 的系数向量 该函数返回方程p x 0的全部根 含重根 复根 2 多项式求值求多项式p x 在某点或某些点的函数值的函数是polyval P x 若x为一数值 则求多项式在该点的值 若x为向量或矩阵 则对向量或矩阵中的每个元素求其多项式的值 例1已知一个多项式 1 计算f x 0的全部根 2 由方程f x 0的根构造一个多项式g x 并与f x 进行对比 3 计算f 5 f 7 8 f 9 6 f 12 3 的值 P 3 0 4 5 7 2 5 X roots P 求方程f x 0的根G poly X 求多项式g x X0 5 7 8 9 6 12 3 f polyval P X0 求多项式f x 在给定点的值 3 多项式的四则运算 1 多项式的加减法 注 多项式求值还有一个函数是polyvalm 其调用格式与polyval相同 但含义不同 polyvalm函数要求x为方阵 它以方阵为自变量求多项式的值 functionp3 poly add p1 p2 n1 length p1 n2 length p2 ifn1 n2p3 p1 p2 endifn1 n2 p3 p1 zeros 1 n1 n2 p2 endifn1 n2 p3 zeros 1 n2 n1 p1 p2 end 加法 c poly add a b 减法 c poly add a b 2 多项式的乘法函数conv P1 P2 用于求多项式P1和P2的乘积 3 多项式的除法函数 Q r deconv P1 P2 用于对多项式P1和P2作除法运算 其中Q返回多项式P1除以P2的商式 r返回P1除以P2的余式 这里 Q和r仍是多项式系数向量 deconv是conv的逆函数 即有P1 conv P2 Q r 例2设有两个多项式 计算 1 求f x g x f x g x 2 求f x g x f x g x f 3 5 2 7 5 6 g 3 5 3 poly add f g 求f x g x poly add f g 求f x g x conv f g 求f x g x Q r deconv f g 求f x g x 商式送Q 余式送r 4 多项式的微分与积分 1 对多项式求导数的函数是 p polyder P 求多项式P的导函数p polyder P Q 求P Q的导函数 p q polyder P Q 求P Q的导函数 导函数的分子存入p 分母存入q 2 对多项式的积分函数 d poly itg c d是多项式c积分后的系数 但不包括积分常数 functionpy poly itg p n length p py p n 1 1 1 0 例3求有理分式的导数 P 3 5 0 8 1 5 Q 10 5 0 0 6 0 0 7 1 0 100 p q polyder P Q 若求多项式P的积分 c poly itg P 4 2MATLAB插值 通常取为多项式函数 代数插值 多项式插值 已知f x 在点xi上的函数值yi f xi i 0 1 2 n 则称P x 为f x 的n次代数插值多项式 称x0 x1 xn为插值结点 P x 为插值函数 条件P xk yk k 0 1 n 为插值条件 f x 为被插值函数 如果P x a0 a1x anxn满足 P xk yk k 0 1 n 设f x C a b 取点a x0 x1 xn b 代数插值问题 定理 若插值结点x0 x1 xn是 n 1 个互异点 则满足插值条件P xk yk k 0 1 n 的n次插值多项式P x a0 a1x anxn存在而且是唯一的 方程组系数矩阵取行列式 故方程组有唯一解 从而插值多项式P x 存在而且是唯一的 构造3次多项式P x 逼近Erf x 设P x a0 a1x a2x2 a3x3 令P xk Erf xk 得 求解 得a0 0 a1 1 293 a2 0 5099 a3 0 0538所以 P x 1 293x 0 5099x2 0 0538x3 MATLAB计算程序x 0 6 1 8 y erf x x x y y A ones 4 1 xx 2x 3 p A y a0 p 1 a1 p 2 a2 p 3 a3 p 4 t 0 2 2 u a0 a1 t a2 t 2 a3 t 3 plot x y o t u 注 一次多项式插值 过两点直线 二次多项式插值 过三点抛物线 不用待定系数法 1 计算量大 2 不易讨论误差 几何意义 两条曲线有交点 公共点 一 线性插值 线性插值是两个数据点的直线拟合 或 误差估计 在MATLAB中 命令interp1可做线性插值 调用格式为 yi interp1 x y xi 其中x表示数据点横坐标的列数组 y表示数据纵坐标的列数组 可以有多列 另外 interp1命令有三种可选参数 yi interp1 x y xi linear 线性插值 缺省 yi interp1 x y xi spline 三次样条yi interp1 x y xi cubic 三次插值 例3已知数据表如下 分别求y 0 9 0 7 0 6 0 5处x的值 x 0 0 0 25 0 5 0 75 1 0 y 0 9126 0 8109 0 6931 0 5596 0 4055 yi 0 9 0 7 0 6 0 5 xi interp1 y x yi linear yi xi ans 0 90000 03100 70000 48540 60000 67430 50000 8467 二 用幂级数做多项式插值 给定n 1个数据点 过n 1个点的n阶多项式可写为幂级数形式 注 过n 1个数据点的n阶插值多项式是唯一的 对n 1个数据点 设 则得到n 1个线性方程 可以表示为矩阵形式 求解该方程组可确定系数 或用polyfit x y n 确定 例3求下列数据点拟合多项式的系数 并求当x 2 101和4 234处y的值 并画出数据点和曲线 x 1 1 2 3 3 9 5 1 y 3 887 4 276 4 651 2 117 n length x 1 a n 1 ones size x a n x forj n 1 1 1a j a j 1 x endcoeff a y xi 2 101 4 234 yi zeros size xi fork 1 n 1yi yi coeff k xi n 1 k endxp 1 1 0 05 5 1 yp zeros size xp fork 1 n 1yp yp coeff k xp n 1 k endplot xp yp x y ro 三 Lagrange插值多项式 1 插值基函数 形状函数的图形如下 n 8 2 Lagrange插值多项式 functionfi Lagran x f xi fi zeros size xi np1 length f fori 1 np1z ones size xi forj 1 np1ifi jz z xi x j x i x j endendfi fi z f i endreturn 3 MATLAB程序实现 调用格式 yi Lagran x y xi clearx 1 1 2 3 3 9 5 1 y 3 887 4 276 4 651 2 117 xi 2 101 4 234 yi Lagran x y xi 例4 写出拟合下面三个数据点的Lagrange插值公式 并计算x 2 101 4 234时y的值 4 截断误差 例5用的5个等距点对函数进行插值估计 插值结果及误差分布如下图 可见 误差峰值出现在端点附近的区间里 这是由于的局部峰值在端点附近 Runge反例 5 x 5 取xk 5 k计算 f xk k 0 1 10 构造L10 x 取 tk 5 0 05k k 0 1 200 计算 L10 tk x 5 5 y 1 1 x 2 t 5 0 05 5 y1 1 1 t 2 n length t fori 1 nz t i s 0 fork 1 11Lk 1 u x k forj 1 11ifj k Lk Lk z x j u x j endends s Lk y k endy2 i s endplot x y ko t y1 t y2 r 减小误差的方案 1 减小插值区域 即b a 2 增加数据点个数 3 使用可变间距的数据点 Chebyshev点 总结 1 尽可能在小区间上使用多项式插值 2 只能在一定范围内依靠增加插值点个数提高插值精度 如果插值点个数过多往往会适得其反 5 Lagrange插值公式的微分与积分 插值公式 微分 实际上 是一个拟合如下数据点的n阶多项式 一般地 是拟合数据点的多项式 所以 多项式可通过拟合待定数据点的n阶多项式表示为幂级数形式 对所有i 的幂级数形式可用函数shape pw计算 functionp shape pw x np length x forj 1 npy zeros 1 np y j 1 p j polyfit x y np 1 end 如例4也可求解如下 x 1 1 2 3 3 9 5 1 y 3 887 4 276 4 651 2 117 xi 2 101 4 234 np length x p shape pw x s 0 fori 1 nps s p i y i endsyi polyval s xi 结果为 yi 4 14574 3007 为计算Lagrange插值多项式的一阶导数 可用polyder函数将p的每一行转换为一阶导数的系数数组 x 1 1 2 3 3 9 5 1 y 3 887 4 276 4 651 2 117 xi 2 101 4 234 np length x p shape pw x s 0 fori 1 nps s polyder p i y i endyi polyval s xi 结果为 yi 0 6292 1 4004 取x0 x1 x2 求二次函数P x a0 a1 x x0 a2 x x0 x x1 满足条件P x0 f x0 P x1 f x1 P x2 f x2 插值条件引出关于a0 a1 a2方程 四 牛顿插值问题 解下三角方程组过程中引入符号 a0 f x0 a1 f x1 x2 a2 f x0 x1 x2 P x a0 a1 x x0 a2 x x0 x x1 定义 若已知函数f x 在点x0 x1 xn处的值f x0 f x1 f xn 如果i j 则 n阶均差 例由函数表求各阶均差 解 按公式计算一阶差商 二阶差商 三阶差商如下 MATLAB程序计算 56 16 2 24 56401403 5640 13 71 5640 1322 5640 1320 x 2 1013 y 56 16 2 24 f yn length x fork 2 nforj n 1 kf j f j f j 1 x j x j 1 k endfend 另外一个程序 x 2 1013 y 56 16 2 24 n length x A zeros n n A 1 y forj 2 nfori j nA i j A i j 1 A i 1 j 1 x i x i j 1 endendA 560000 1640000 214 1300 20 72043120 在编写牛顿插值多项式程序之前 复习几个命令 例6 1 2求三个一次多项式积 它们的零点分别依次为0 4 0 8 1 2 解 X1 0 4 0 8 1 2 l1 poly X1 L1 poly2sym l1 运行后输出结果为l1 1 0000 2 40001 7600 0 3840L1 x 3 12 5 x 2 44 25 x 48 125 P1 poly 0 4 P2 poly 0 8 P3 poly 1 2 C conv conv P1 P2 P3 L1 poly2sym C 运行后输出的结果与方法1相同 牛顿插值多项式的函数文件 function A C L newploy X Y n length X A zeros n n A 1 Y forj 2 nfori j nA i j A i j 1 A i 1 j 1 X i X i j 1 endendC A n n fork n 1 1 1C conv C poly X k d length C C d C d A k k endL k poly2sym C 例 给出节点X 2 15 1 000 011 022 033 25 Y 17 037 241 052 0317 0623 05 作五阶牛顿插值多项式和差商 解 X 2 15 1 000 011 022 033 25 Y 17 037 241 052 0317 0623 05 A C L newploy X Y 五 Chebyshev多项式 等距Lagrange插值的误差在中间部分最小 越靠近端点处误差越大 降低误差的一种方法是改变插值点的分布 即增大定义域中间部分的插值步长 同时减少定义域两端的插值步长 但是最优的插值点分布取决于插值多项式的目的 若目的是近似一个函数 则最优插值点Chebyshev多项式的零点 因为 1 此时L x 分布是最平坦的分布 2 误差值不会像等距插值那样随不同插值区间而变化 Lagrange插值余项定理 取插值结点a x0 x1 xn b则满足Ln xk f xk 的n次插值多项式Ln x 误差 其中 选择 x0 x1 xn 使 结论 Chebyshev多项式Tn 1 x 的全部零点 Chebyshev多项式 当时 Chebyshev多项式的图形如下 幂级数形式Chebyshev多项式的系数可由函数Cheby pw计算 functionpn Cheby pw n pbb 1 ifn 0 pn pbb break endpb 10 ifn 1 pn pb break endfori 2 n pn 2 pb 0 0 0 pbb pbb pb pb pn end 调用格式 p Cheby pw n n为多项式阶数 p是系数的行数组 Chebyshev多项式的根系数可由roots命令计算roots Cheby pw n 也可如下求解 令 得 有k个根 都在 1 1 也可用下面公式计算 若插值区间为 a b 可建立映射 下图比较了区间上9个Chebyshev点和等距点条件下L x 的分布 例 函数 取等距插值结点 5 4 3 2 1 0 1 2 3 4 5 x 5 5 w11 x x 5 x 4 x 3 x 2 x 1 x x 1 x 2 x 3 x 4 x 5 w11 x 的图形 4 9491 4 5482 3 7787 2 7032 1 40870 00001 40872 70323 77874 54824 9491 w11 x w11 x x x0 x x1 x x2 x x10 插值函数L10 x 取Chebyshev结点 插值函数L10 x 取等距结点插值 六 Lobatto点 Lobatto点是Chebyshev多项式的零点加上x 1和x 1对区间 a b k 1个Lobatto点可由k阶Chebyshev多项式得到 七 Legendre多项式 Legendre多项式 n阶Legendre多项式的系数可由函数Legen pw计算 functionpn Legen pw n pbb 1 ifn 0 pn pbb break endpb 10 ifn 1 pn pb break endfori 2 n pn 2 i 1 pb 0 i 1 0 0 pbb i pbb pb pb pn end 调用 p Legen pw n n阶Legendre多项式的根可由roots计算 已知节点x0和x1处的函数值及导数值 求三次插值函数 八 三次Hermite插值问题 例 已知插值条件 求3次插值函数 解 设 得a0 0 a1 0 列出方程组 求解 得a2 3 a3 2 所以 有H x 3x2 2x3 3 2x x2 利用基函数表示Hermite插值 定理 两点Hermite插值的误差估计式 反复应用Roll定理 得F 4 t 有一个零点设为 由 得 所以 显然 F t 有三个零点x0 x x1 由Roll定理知 存在F t 的两个零点t0 t1满足x0 t0 t1 x1 而x0和x1也是F x 的零点 故F x 有四个相异零点 当插值节点含有奇点时的处理方法 以三次Hermite多项式为例 三次Hermite多项式 可以拟合两点函数值和导数值 已知两点和及它们的函数值和导数值 可建立方程组 若用有限差分近似表示导数边界条件 可通过Lagrange插值确定该方程组的系数 当y x x y 含有奇点时 常需引入一个参数s 并将x和y表示为x s 和y s 假设点A处s 0 B处s 1 则x y表示为关于s的三次多项式 x s 和y s 的边界条件可写为 其中a b是任意参数 在一定程度上影响曲线形状 例6确定一条经过点A B的曲线 使满足 解 引入s 按上述方法可得到两方程组 分别求解 a 3 b 3 c 0001 0010 1111 3210 1 a 4 0 d 0001 0010 1111 3210 1 0 2 b s 0 0 01 1 x polyval c s y polyval d s plot x y 结果 c 3331 d 1001 若用差分近似 设置四个数据点 取 求解 z 0 01 a 3 b 3 s 1 0 x 1 1 y 1 1 s 2 z x 2 1 z a y 2 1 s 3 1 z x 3 4 y 3 2 z b s 4 1 x 4 4 y 4 2 c polyfit s x length s 1 d polyfit s y length s 1 ss 0 0 1 1 xp polyval c ss yp polyval d ss plot xp yp ylabel y xlabel x 结果是 c 3 90213 12312 96911 d 1 0307 0 03090 00021 九 分段插值问题 结点增多 多项式次数增高 逼近精度越好 未必 多结点高次插值往往在局部误差更大 Runge 龙格 现象 实用 采用分段低次插值有分段线性 分段二次插值等缺点 分段插值函数只能保证连续性 不能保证光滑性 分段插值法可以克服在大范围内使用高次插值带来的误差放大的问题 所谓的分段插值法就是首先在插值区间 a b 上插入节点然后在每个小的子区间 xi 1 xi 上构造低次插值多项式 再将每个子区间上的多项式连接 作为插值区间的插值函数 1 分段线性插值 区间 a b 上的分段线性插值函数为 分段线性插值在MATLAB中的函数调用 yi interp1 x0 y0 xi 1 横纵坐标都是向量 2 x是n维向量 y为矩阵时 行数必须是x的维数 例 x0 6 6 y0 sin x0 xi 6 25 6 yi interp1 x0 y0 xi x 6 0 001 6 y sin x plot x0 y0 o xi yi x y legend 节点 xi yi 分段线性插值函数 被插值函数y sinx title y sinx及其分段线性插值函数和节点的图形 分段Hermite插值 即在每个小区间上用分段的三次Hermite插值 在matlab里有专门的调用函数 y1 interp1 x y x1 pchip 例 x0 6 6 y0 sin x0 xi 6 25 6 yi interp1 x0 y0 xi pchip x 6 0 001 6 y sin x plot x0 y0 o xi yi x y legend 节点 xi yi 分段埃尔米特插值函数 被插值函数y sinx title y sinx及其分段埃尔米特插值函数和节点的图形 x 5 5 y 1 1 x 2 plot x y x y o x 5 5 y 1 1 x 2 xi 5 05 5 yi spline x y xi plot xi yi b x y ro 被插值函数 5 x 5 十 样条插值 分段线性插值可以得到整体连续函数 但在连接点处一般不光滑 而Hermite插值虽然在连接点处一阶光滑 但整体插值由于结点多次数高而有可能发生龙格现象 若用分段三次Hermite插值 节点处一阶导数连续光滑 但在曲线的凹凸性变化比较大的地方 误差较大 既想分段插值 又想在结点处保持光滑 甚至二阶光滑 三次样条 希望 定义 设有n 1个点的点列 若函数满足 3 在整个区间内具有二阶连续导数 2 在每个小区间 上是三次多项式 此时叫插值函数 则称为点列的三次样条插值函数或三次样条多项式 简称三次样条 定义中的一阶导数连续意味着曲线没有急弯 二阶导数连续意味着曲线每一点的曲率半径有意义 则样条函数可写为 考虑区间 令 列方程组 求解得 所以 则在区间上 同理 在区间上 由在节点处连续 式 相等 消去 得 式除了对第一个点和最后一个点不成立外 其它各点均成立 所以共有n 2个方程 而待定的个数为n 所以还需两个额外条件 可由下述边界条件提供 常见到的边界条件有三种 a 在端点处指定 若已知 则方程组为 该方程组为含n 2个未知量的n 2个方程 其系数矩阵为三对角阵 可用追赶法求解 大多数情形中 可令 等价于在几何上趋于端点的曲线变为直线 b 从内部外推 利用对进行外推 将其代入方程组的第一个方程中 整理得 类似 利用对进行外推 将其代入方程组的最后一个方程 整理得 此时 同样得到一三对角方程组 求解即可 c 循环边界条件适用于封闭曲线 即第一个点和最后一个点相同 它们的导数值也相同 此时相当于有n 1个点 相应有n 1个方程 可以求解 在MATLAB中 三次样条插值运算实现如下 yi interp1 x y xi spline 或yi spline x y xi 其中x y都是向量形式的点 xi是进行插值的点的横坐标向量 yi为插值函数值 例7求三次样条插值并作图比较 以s为参数 分别做x y的样条函数xx 1 0 866 0 500 50 866111 04021 15001 3 1 541 82802 17362 58833 0860 yy 0 0 25 0 433 0 5 0 433 0 25000 150 25980 3 0 30 30 30 30 3 s 1 length xx sp 1 length xx 100 length xx xp spline s xx sp yp spline s yy sp plot xp yp holdonplot xx yy ro xlabel x ylabel y axis 13 5 11 若不取s为参数 直接以x y做样条函数 曲线不光滑 xi 1 1 100 3 5 yi spline xx yy xi plot xi yi xx yy ro x 0 0 0155 0 1485 0 3493 0 6480 1 0547 2 0 y 0 0 1242 0 3654 0 4975 0 5472 0 4781 0 n length x t 0 n 1 tt 0 25 n 1 xx spline t x tt yy spline t y tt plot xx yy x y o 二维插值是基于与一维插值同样的基本思想 然而 正如名字所隐含的 二维插值是对两变量的函数z f x y 进行插值 为了说明这个附加的维数 考虑一个问题 设人们对海水的温度分布估计感兴趣 给定的温度值取自海水表面均匀分布的格栅 采集了下列的数据 width 1 5 indexforwidthofplate i e thex dimension depth 1 3 indexfordepthofplate i e they dimension 十一 二维插值 temps 8281808284 7963616581 8484828586 temperaturedatatemps 828180828479636165818484828586矩阵temps表示整个海水的温度分布 temps的列与下标depth或y 维相联系 行与下标width或x 维相联系 为了估计在中间点的温度 我们必须对它们进行辨识 MESHGRID命令 调用格式一 X Y meshgrid x y X Y meshgrid 3 2 3 3 2 3 Z 7 3 X 4 exp X 2