第二章 函数的插值
第二章 插值法
证明：
yi , i = 0, ... , n 的 n 阶插值多
反证：若不唯一，则除了Ln(x) 外还有另一 n 阶多项 式 Pn(x) 满足 Pn(xi) = yi 。 考察 Qn ( x ) = Pn ( x ) - Ln ( x ) , 则 Qn 的阶数 n 而 Qn 有 n + 1 个不同的根 x0 … xn 注：若不将多项式次数限制为 n ，则插值多项式不唯一。 例如 P ( x ) = Ln ( x ) p( x ) ( x - x i ) 也是一个插值
sin 50 0 L2 ( 5 ) 0.76543 18
R2 ( x ) = - cos x ( x - )( x - )( x - ) ; 3! 6 4 3 1 cos 3 x 2 2
0.00044 R2 5 0.00077 18
=
x - x1 y + x 0 - x1 0
x - x0 y = x1 - x 0 1
l ( x) y
i =0 i
1
i
l0(x)
l1(x)
§1 Lagrange Polynomial
n1
希望找到li(x)，i = 0, …, n 使得 li(xj)=ij ；然后令
Pn ( x ) =
l (x) y
g(x) f(x)
x0
x1
x2
x
x3
x4
§1 拉格朗日多项式
Pn ( x i ) = y i ,
/* Lagrange Polynomial */
n 求 n 次多项式 Pn ( x ) = a0 a1 x an x 使得
i = 0 , ... , n
第二章插值法
lk ( xk 1 ) 0
n=2的情况，假定插值节点为
xk 1 , xk , xk 1 , 要求一个二次插值多项式L2 ( x)，使它满足 L2 ( x j ) y j ( j k 1, k , k 1)
y L2 ( x)在几何上就是通过三点(xk-1 , yk 1 ),(xk , yk ),(xk+1, yk 1 )的抛物线
插值法
§2.1 §2.2 §2.3 §2.4 §2.5 §2.6 §2.7 引言 拉格朗日插值 均差与牛顿插值公式 差分与等距节点插值 埃尔米特插值 分段低次插值 三次样条插值
一、插值问题
或者函数本身只是 一组实验数据，很 难对函数的性质进 行分析
对函数f (x),其函数形式可能很复杂且不利于在计算机上 ,
设函数
y f ( x ) 在区间 [a, b] 上有定义，且已知在
a x0 x1 x2 xn b
f ( xi ) yi , i 0,1,, n
如果存在一个简单函数 P ( x )，使得
P( xi ) f ( xi ) yi , i 0,1,, n
xx x x
如函数y sin x, 若给定 0, ]上5个等分点 [
其插值函数的图象如图
对于被插函数 ( x)和插值函数 ( x) f P
在节点xi处的函数值必然相等
但在节点外 ( x)的值可能就会偏离 ( x) P f 因此P( x)近似代替 ( x)必然存在着误差 f
整体误差的大小反映了插值函数的好坏
成立，则称 P ( x ) 为 f ( x ) 的插值函数
称点 xi , i 0,1,2,, n为插值节点
称区间 a , b]为插值区间 [
第2章_插值法
13.214 285 71
175 13.228756555322952...
考虑通过 + 1个节点0 < 1 < ⋯ < 的次插值
多项式 ()，满足条件
= ,
= 0,1, … ,
希望找到 li(x)，i = 0, …, n, 使得
= ; = ,
n次插值多项式, 插值节点为{ xi }in 0 [ a , b],则x [ a , b],有
f ( n 1) ( )
Rn (x )
n 1 ( x)
Lagrange型余项
(n 1)!
n
其中 n 1 ( x ) ( x xi ) , ( a , b) , 且依赖于 x.
满足条件P(xi) = f(xi) (i = 0, … n)。 P(x) 称
为f(x) 的插值函数。
P(x) f(x)
x0
x1
x2
x
x3
x4
定理1：设插值节点 ≠ ( ≠ )，则满足条件
= , = 0,1, … , 的插值多项式
= 0 + 1 + ⋯ +
− , , + 线性无关。
二次插值多项式
= − − + + + + ()
满足 = ( = − , , + )
例1:
已知 f ( x )满足 f (144) 12 , f (169) 13, f ( 225) 15
i 0
一次及二次差值余项
1 ′′
1 = − 0 − 1 ,
计算方法—插值法 (课堂PPT)
7
1 1
2 5
4 25
8 125
aa32
4
35
则，
解方程组得a0=10,a1=5,a2=-10,a3=2 即P3(x)=10+5x-10x2+2x3
当n=20,在109次/秒的计算机上计算需几万年！
.
2020/4/2
12
2.2 拉格朗日插值
2-2 线性插值与抛物插值
Chapter2 插值法
第二章 插 值 法
（ Interpolation） 2.1 引言
2.2 拉格朗日插值
2.3 均差与牛顿插值公式
Chapter2 插值法
2.4 埃尔米特插值
2.5 分段低次插值
2.6 三次样条插值
.
2020/4/2
1
2.1 引言
Chapter2 插值法
表示两个变量x,y内在关系一般由函数式 y=f(x)表达。但在实际问题中的函数是多种多 样的,有下面两种情况：
几何意义:L2(x)为过三点(x0,y0), (x1,y1), (x2,y2)的抛物线。
方法：基函数法，构造基函数l0(x), l1(x), l2(x) (三个二次式)
使L2(x)= y0l0(x)+y1l1(x)+y2l2(x)满足插值条件。 6 4 4 4 4 4 4 7 4 4 4 4 4 48
.
2020/4/2
15
2.2 拉格朗日插值
Chapter2 插值法
问题的提法： 已知y=f(x)的函数表，x0, x1, x2为互异节
x x0 x1 x2 y y0 y1 y2
点，求一个次数不超过2的多项式 L2(x)=a0+a1x+a2x2 ：L2(x0)=y0, L2(x1)=y1, L2(x2)=y2
第二章插值与拟合
1 不为零。
xn
n xn xn
实 用 测 量 数 据 处 理 方 法
中 南 大 学
三、线性插值
假定已知区间[xk, xk+1] 的端点处的函数值 yk=f(xk), yk+1=f(xk+1),要求线性插值多项式 L1(x),使它满足 L1(xk)=yk
L1(xk+1)=yk+1
则L1(x)的表达式可按下式给出：
实 用 测 量 数 据 处 理 方 法
中 南 大 学
l k 1 ( x k 1 ) 1, l k 1 ( x j ) 0( j k , k 1) l k ( x k ) 1, l k ( x j ) 0( j k 1, k 1) (28) l k 1 ( x k 1 ) 1, l k 1 ( x j ) 0( j k 1, k ) 满足（28 ）式的插值基函数很容 易求出的，例如求 l k 1 ( x)，因为它有两个零点 k 和x k 1，故可表达为： x l k 1 ( x) A( x x k )(x x k 1 ) 其中A为待定系数可由 k 1 ( x k 1 ) 1定出： l 1 A ( x k 1 x k )(x k 1 x k 1 ) ( x x k )(x x k 1 ) 于是l k 1 ( x)＝ ，同理可得 ( x k 1 x k )(x k 1 x k 1 ) ( x x k－1 )(x x k 1 ) ( x x k 1 )(x x k ) l k ( x)＝ ，l k 1 ( x)＝ ( x k x k－1 )(x k x k 1 ) ( x k＋1 x k 1 )(x k 1 x k )
解：2、抛物插值
计算方法(2)-插值法
2
2
yk1 2
f (xk
h
2
),
y
k
1 2
f (xk
h) 2
21
3.牛顿向后插值公式
Nn (xn
th)
yn
tyn
t(t 1) 2!
2
yn
t(t
1)
(t n!
n
1)
n
yn
(t 0)
插值余项
Rn
(xn
th)
t(t
1) (t (n 1)!
Nn (x0
th)
y0
ty0
t(t 1) 2!
2
y0Leabharlann 插值余项t(t
1)
(t n!
n
1)
n
y0
Rn (x0
th)
t(t
1) (t (n 1)!
n)
h n1
f
(n1) ( ),
(x0 , xn )
20
二.向后差分与牛顿向后插值公式
杂.
根据f(x)函数表或复杂的解析表达式构
造某个简单函数P(x)作为f(x)的近似.
2
2.问题的提法
1)已知条件 设函数y f (x)在区间[a,b]上
连 续, 且 在n 1个不 同点a x0 , x1, , xn b 上 分 别 取 值y0 , y1, , yn
第二章插值法
线性插值的几何意义:用 通过点 A(x0 , f (x0 )) 和 B(x1, f (x1 )) 的直线近似地代替曲线 y=f(x)由解析几何知道, 这条直线用点斜式表示为
y=f(x)
p(x)=ax+b
A(x.0,f(x.0)) B(x.1,f(x.1))
p(x)
y0
y1 x1
y0 x0
(x
x0 )
解: 用待定系数法, 将各节点值依次代入所求多项式, 得
解上述方程, 将求出的a0, a1, a2 代入 p(x) = a0 + a1x + a2x2 即得所求二次多项式
例2.5 求过点(0,1)、(1,2)、(2,3)的三点插值多项式
P(x) an x n an1 x n1 a1x a0
满足
P(x) an x n an1 x n1 a1 x a0 P(xi ) f (xi ) (i 0,1,2, , n)
则称P(x)为f(x)的n次插值多项式。这种插值法通常称
为代数插值法。其几何意义如下图所示
y y=P(x) y=f(x)
lk ( xi ) ki
1 0
(i k ) (i k )
l0 (x) 与 l1(x) 称为线性插值基函数。且有
lk (x)
1 j0
x xj , xk x j
jk
k 0,1
于是线性插值函数可以表示为与基函数的线性组合
p(x) l0 (x) y0 l1 (x) y1
例2.1 已知 100 10 , 121 11 , 求 y 115
为已知 f (x0 ), f (x1), , f (xn ) ,即 yi f (xi ) 若存在一个
f(x)的近似函数 (x),满足
插值法
第一节 Lagrange插值
一、问题提出
设 x0 , x1 xn 为给定的节点，yi f ( xi )，i 0,1,n
为相应的函数值，求一个次数不超过 n 的多项式 Pn (x)， 使其满足
Pn ( xi ) yi，
i 0,1,n .
这类问题称为插值问题。 f ( x) 称为被插值函数，Pn ( x) 称 为插值函数， x0 , x1 xn 称为插值节点
差商
二阶差商
三阶差商 四阶差商
x0 f ( x0 ) x1 f ( x1 )
x2 f ( x2 )
f [ x0 , x1 ]
f [ x1 , x2 ]
f [ x0 , x1 , x2 ]
f [ x0 , x1 , x2 , x3 ]
1 2 3 4
0 1 2 3 4
x3
f ( x3 ) f [ x2 , x3 ] f [ x1 , x2 , x3 ]
评价
优点: Lagrange基函数容易构造,结构紧凑,便于理 论研究. 缺点: 当增加或减少插值结点时,基函数需要重新 构造,不便于实际的计算使用
第二节 Newton插值
一、差商定义及性质
1.差商定义 f ( x ) f ( x ) i j f [ xi , x j ] , i j 为 f ( x) 在 xi , x j 称 两点处的一阶差商.xi x j
( n1) ( ) f ( n1) ( )
f ( x) Pn ( x) (n 1)! 0 ( x)
由此得
. f ( n1) ( ) Rn ( x) f ( x) Pn ( x) n1 ( x) (n 1)! 定理得证.
《数值分析》第二讲插值法PPT课件
1 xn xn2 xnn Vandermonde行列式
即方程组（2）有唯一解 (a0, a1, , an)
所以插值多项式
P (x ) a 0 a 1 x a 2 x 2 a n x n
存在且唯一
第二章：插值
§2.2 Lagrange插值
y
数值分析
1、线性插值
P 即(x)ykx yk k 1 1 x yk k(xxk)
l k ( x k 1 ) 0 ,l k ( x k ) 1 ,l k ( x k 1 ) 0 l k 1 ( x k 1 ) 0 ,l k 1 ( x k ) 0 ,l k 1 ( x k 1 ) 1
lk1(x)(x(k x 1 x xk k))x x ((k 1x k x 1k )1) lk(x)((xx k x xk k 1 1))((x xkxx k k1)1)
第二章：插值
数值分析
3、Lagrange插值多项式
令 L n ( x ) y 0 l 0 ( x ) y 1 l 1 ( x ) y n l n ( x )
其中，基函数
lk (x ) (x ( k x x x 0 ) 0 ) (( x x k x x k k 1 1 ) )x x k ( ( x x k k 1 ) 1 ) (( x x k x n x )n )
因此 P (x ) lk (x )y k lk 1 (x )y k 1
且
P (x k ) y k P (x k 1 ) y k 1
lk(x), lk1(x) 称为一次插值基函数
数值分析
第二章：插值
2、抛物线插值 令
y (xk , yk )
f (x)
lk1(x)(x(k x 1 x xk k))x x ((k 1x k x 1k )1) p( x) (xk1,yk1)
计算方法 第二章插值法_20191105
下面两种办法常用来确定经验函数y=g(x)
（1）插值法
（2）拟合法
根据问题的不同，有时要用插值技术来解决， 有时则应该采用拟合的方法才合理。
（1）插值法的基本思想
已知数据表
xi x1 x2 … xn f(xi) f(x1) f(x2) … f(xn)
求一个经验函数y=g(x)，使g(xi)=f(xi), i=1,…n.
x)
b0
a0 a1 x a2 x2 b1 x b2 x2 b3 x3
n
一般地：F( x) cii( x) i=0
例：F(x) a bx csin x span1, x,sin x,
当插值函数是代数多项式时，插值问题称为代
代 数插值。
数 插
设 Pn(x)=a0+a1x+…+anxn ,
y2
n 次插值多项式 :求次数≤n的多项式Ln(x), 使其满足
Ln(x0)=y0 , Ln(x1)=y1 , ...... , Ln(xn)=yn
..(7)
令 Ln(x)=l0(x)y0 + l1(x)y1 +… +ln(x)yn
求n 次多项式 lj(x) ,(j=0,1,…,n)使其满足条件
容易求得
三角插值：取
spani(
n
x) i=0
a1x a0
=spansin x,cos x,sin 2x,cos 2x, ,sin nx,cos nx
例：取 spansin x,cos x,
F ( x) a sin x bcos x
有理插值：F( x)= Pm ( x) Qn ( x)
例：F (
二次插值 (n=2) 求次数≤2 的多项式L2(x), 使其满足条件 L2(x0)=y0 , L2(x1)=y1 , L2(x2)=y2
数值分析--第2章 插值法
数值分析--第2章插值法第2章 插值法在科学研究与工程技术中，常常遇到这样的问题：由实验或测量得到一批离散样点，要求作出一条通过这些点的光滑曲线，以便满足设计要求或进行加工。
反映在数学上，即已知函数在一些点上的值，寻求它的分析表达式。
此外，一些函数虽有表达式，但因式子复杂，不易计算其值和进行理论分析，也需要构造一个简单函数来近似它。
解决这种问题的方法有两类：一类是给出函数)(x f 的一些样点，选定一个便于计算的函数)(x ϕ形式，如多项式、分式线性函数及三角多项式等，要求它通过已知样点，由此确定函数)(x ϕ作为)(x f 的近似，这就是插值法；另一类方法在选定近似函数的形式后，不要求近似函数过已知样点，只要求在某种意义下在这些样点上的总偏差最小。
这类方法称为曲线（数据）拟合法。
设已知函数f 在区间],[b a 上的1+n 个相异点ix 处的函数值(),0,,iif f x i n ==，要求构造一个简单函数()x ϕ作为函数()f x 的近似表达式()()f x x ϕ≈，使得()(),0,1,,iiix f x f i n ϕ=== (2-1) 这类问题称为插值问题。
称f 为被插值函数；()x ϕ为插值函数；nx x ,,0 为插值节点；(2-1)为插值条件。
若插值函数类{()}x ϕ是代数多项式，则相应的插值问题为代数插值。
若{()}x ϕ是三角多项式，则相应的插值问题称为三角插值。
若{()}x ϕ是有理分式，则相应的插值问题称为有理插值。
§1 Lagrange 插值1.1 Lagrange 插值多项式设函数f 在1+n 个相异点01,,,nx x x 上的值n i x f f ii ,,1,0),( ==是已知的，在次数不超过n 的多项式集合n P 中，求()nL x 使得(),0,1,,n i iL x f n n == (2-2) 定理2.1 存在惟一的多项式nn P L ∈满足插值条件(2-2)。
第二章-函数插值-1
大多数实际问题都可用函数来表示某种内在规律的数量关系 但函数通常无法给出，只有通过实验或观测得到的数据表 如何根据这些数据推测或估计其它点的函数值？ 例：已测得在某处海洋不同深度处的水温如下：
深度（M） 466 741 4.28 950 3.40 1422 2.54 1634 2.13 水温（oC） 7.04
常用插值方法
P(x)
x0
x1
x2
x
x3
x4
如果存在一个简单易算的函数 p(x) ，使得
p(xi) = f(xi)，i = 0, 1, ... , n
则称 p(x) 为 f(x) 的插值函数 [a, b] 为插值区间，xi 为插值节点，p(xi) = f(xi) 为插值条件 插值节点无需递增排列，但必须确保互不相同！ 求插值函数 p(x) 的方法就称为插值法
多项式插值
多项式插值 已知函数 y = f(x) 在 [a, b] 上 n + 1 个点 a x0 < x1 < ꞏ ꞏ ꞏ < xn b 处的函数值为 y0 = f(x0)，… ，yn = f(xn)
求次数 不超过 n 的多项式
p(x) = c0+c1x + ꞏ ꞏ ꞏ + cnxn，
设 p(x) = a0l0(x) + a1l1(x) + ꞏ ꞏ ꞏ + anln(x) 将 p(xi) = yi ，i = 0, 1, ... , n 代入，可得
ai = yi ，i = 0, 1, ... , n p(x) = y0l0(x) + y1l1(x) + ꞏ ꞏ ꞏ + ynln(x) Ln(x) 就称为 f(x) 的 Lagrange 插值多项式
1_多项式插值基本理论
yfnie@
10
二、插值多项式的构造方法
由于插值多项式的存在惟一性，无论是用何种 方法构造出的插值多项式，它们均恒等，进而 截断误差也都相同。 内容提要
Lagrange插值法 Newton插值法 等距节点插值公式 带导数的插值问题
August 6, 2012
August 6, 2012
yfnie@
13
1.2 Lagrange型插值公式
L n ( x ) f ( x i )li ( x )
i0
n
n
f ( xi )
n 1 ( x )
( x x i ) n 1 ( x i )
'
i0
不超过n次的多项式，满足所有的插值条件，因而是需 构造的插值多项式，称之为Lagrange插值多项式。
( x 1 )( x 0 ) ( 0 . 5 1 )( 0 . 5 0 )
2 3
x ( x 0 .5 )
2 ( x 1 )( x 0 . 5 )
4 3
x ( x 1)
二次插值多项式为
L 2 ( x ) f ( x 0 ) l 0 ( x ) f ( x 1 ) l1 ( x ) f ( x 2 ) l 2 ( x ) l 0 ( x ) 2 l1 ( x ) 3 l 2 ( x )
第二章 函数插值
问题提出
1 函数表达式过于复杂不便于计算, 而又需要计算许多点处的函数 值
2 仅有采样值, 而又需要知道非采样点处的函数值 ……
上述问题的一种解决思路：建立复杂函数或者未知函数的一个便于计 算的近似表达式.
内容提要
第二章：插值法
满足(2.1)式的 l i(x) 是否存在?若存在，具有什么形式呢?
先考虑 l0(x)。因 l0(x)是以 x1, x2 为零点的二次多项式,
所以它可写成 l0(x)＝ 0(x －x1)(x －x2), 其中0 是待定系 数。 又因为 l0( x0)=1，所以0(x0－x1)(x0－x2)＝1，则可有
n
| x - xi |
i=0
作为误差估计上限。
当 f(x) 为任一个次数 n 的多项式时， f (n1)( x) 0,
可知 Rn ( x) 0 ，即插值多项式对于次数 n 的多项式 是精确的。
例1 求经过A(0,1)，B(1,2)，C(2,3)三个插值点的插值多项式. 解：三个插值节点及对应的函数值为
-
3
);
1 2
cos x
3 2
0.00044
R2
5
18
0.00077
sin 50 = 0.7660444…
2次插值的实际误差 0.00061
高次插值通常优于 低次插值
但绝对不是次数越 高就越好，嘿 嘿……
例3 考虑下述的插值法问题：求二次多项式P(x)，满足 P(x0) = y0, P(x1) = y1，P(x2 ) = y2， 其中 x0 x2，y0、y1、y2 是已给的数据并给出使这一问题的解存在且唯一的条件.
x0 )(x -
x1 )，
[ x0 , x1 ]
当n = 2时 ， 抛 物 插 值 的 余 项 为
R2 ( x) =
1 6
f ( )( x -
x0 )(x -
x1 )(x -
x2 )，
[x0 , x2 ]
注： 通常不能确定 x , 而是估计 f (n1)( x) Mn1 , x(a,b)
第二章函数的插值方法
第⼆章函数的插值⽅法2.1 插值问题及其误差2.1.2 与插值有关的MATLAB 函数(⼀) POLY2SYM 函数调⽤格式⼀：poly2sym (C)调⽤格式⼆：f1=poly2sym(C,'V') 或 f2=poly2sym(C, sym ('V') ),(⼆) POLYVAL 函数调⽤格式：Y = polyval(P ,X)(三) POLY 函数调⽤格式：Y = poly (V)(四) CONV 函数调⽤格式：C =conv (A, B)例 2.1.2 求三个⼀次多项式)(x g 、)(x h 和)(x f 的积)()()(x h x g x f ??.它们的零点分别依次为0.4,0.8,1.2.解我们可以⽤两种MATLAB 程序求之.⽅法1 如输⼊MATLAB 程序>> X1=[0.4,0.8,1.2]; l1=poly(X1), L1=poly2sym (l1)运⾏后输出结果为l1 =1.0000 -2.4000 1.7600 -0.3840L1 =x^3-12/5*x^2+44/25*x-48/125⽅法2 如输⼊MATLAB 程序>> P1=poly(0.4);P2=poly(0.8);P3=poly(1.2);C =conv (conv (P1, P2), P3) , L1=poly2sym (C)运⾏后输出的结果与⽅法1相同.(五) DECONV 函数调⽤格式：[Q,R] =deconv (B,A)(六) roots(poly(1:n))命令调⽤格式：roots(poly(1:n))(七) det(a*eye(size (A)) - A)命令调⽤格式：b=det(a*eye(size (A)) - A)2.2 拉格朗⽇（Lagrange ）插值及其MATLAB 程序估计其误差.解输⼊程序>> X=[1,3];Y=[1,2]; l01= poly(X(2))/( X(1)- X(2)), l11=poly(X(1))/( X(2)- X(1)), l0=poly2sym (l01),l1=poly2sym(l11), P = l01* Y(1)+ l11* Y(2),L=poly2sym (P),x=1.5; Y = polyval(P,x)运⾏后输出基函数l 0和l 1及其插值多项式的系数向量P （略）、插值多项式L 和插值Y 为l0 = l1 = L = Y =-1/2*x+3/2 1/2*x-1/2 1/2*x+1/2 1.2500输⼊程序>> M=5;R1=M*abs((x-X(1))* (x-X(2)))/2运⾏后输出误差限为R1 =例2.2.2 求函数=)(x f e x -在]1,0[上线性插值多项式，并估计其误差.解输⼊程序>> X=[0,1]; Y =exp(-X) ,l01= poly(X(2))/( X(1)- X(2)),l11= poly(X(1))/( X(2)- X(1)), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),运⾏后输出基函数l 0和l 1及其插值多项式的系数向量P 和插值多项式L 为l0 = l1 = P =-x+1 x -0.6321 1.0000L =-1423408956596761/2251799813685248*x+1输⼊程序>> M=1;x=0:0.001:1; R1=M*max(abs((x-X(1)).*(x-X(2))))./2运⾏后输出误差限为R1 =0.1250.2.2.2 抛物线插值及其MATLAB 程序例 2.2.3 求将区间 [0, π/2] 分成n 等份)2,1(=n ，⽤x x f y cos )(==产⽣1+n 个节点，然后根据（6.9）和（6.13）式分别作线性插值函数)(1x P 和抛物线插值函数)(2x P .⽤它们分别计算cos (π/6) (取四位有效数字)，并估计其误差.解输⼊程序>> X=[0,pi/2]; Y =cos(X) ,l1=poly2sym (l11),P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=pi/6;Y = polyval(P,x)运⾏后输出基函数l 0和l 1及其插值多项式的系数向量P 、插值多项式和插值为l0 =-5734161139222659/9007199254740992*x+1l1 =5734161139222659/9007199254740992*xP =-0.6366 1.0000L =-5734161139222659/9007199254740992*x+1Y =0.6667输⼊程序>> M=1;x=pi/6; R1=M*abs((x-X(1))*(x-X(2)))/2运⾏后输出误差限为R1 =0.2742.（2）输⼊程序>> X=0:pi/4:pi/2; Y =cos(X) ,l01= conv (poly(X(2)),poly(X(3)))/(( X(1)- X(2))* ( X(1)- X(3))),l11= conv (poly(X(1)),poly(X(3)))/(( X(2)- X(1))* ( X(2)- X(3))),l21= conv (poly(X(1)),poly(X(2)))/(( X(3)- X(1))* ( X(3)- X(2))),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21),Y = polyval(P,x)运⾏后输出基函数l 01、l 11和l 21及其插值多项式的系数向量P 、插值多项式L 和插值Y 为l0 =228155022448185/281474976710656*x^2-2150310427208497/1125899906842624*x+1 l1 =-228155022448185/140737488355328*x^2+5734161139222659/2251799813685248*x228155022448185/281474976710656*x^2-5734161139222659/9007199254740992*xP =-0.3357 -0.1092 1.0000L=-6048313895780875/18014398509481984*x^2-7870612110600739/72057594037927936*x+1Y =0.8508输⼊程序>> M=1;x=pi/6; R2=M*abs((x-X(1))*(x-X(2)) *(x-X(3)))/6运⾏后输出误差限为R2 =0.0239.2.2.3 n 次拉格朗⽇（Lagrange ）插值及其MATLAB 程序例 2.2.4 给出节点数据00.17)00.2(=-f ，00.1)00.0(=f ，00.2)00.1(=f ，00.17)00.2(=f ，作三次拉格朗⽇插值多项式计算)6.0(f ，并估计其误差.解输⼊程序>> X=[-2,0,1,2]; Y =[17,1,2,17];p1=poly(X(1)); p2=poly(X(2));p3=poly(X(3)); p4=poly(X(4));l01= conv ( conv (p2, p3), p4)/(( X(1)- X(2))* ( X(1)- X(3))* ( X(1)- X(4))),l11= conv ( conv (p1, p3), p4)/(( X(2)- X(1))* ( X(2)- X(3))* ( X(2)- X(4))),l21= conv ( conv (p1, p2), p4)/(( X(3)- X(1))* ( X(3)- X(2))* ( X(3)- X(4))),l31= conv ( conv (p1, p2), p3)/(( X(4)- X(1))* ( X(4)- X(2))* ( X(4)- X(3))),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),运⾏后输出基函数l 0，l 1，l 2和l 3及其插值多项式的系数向量P （略）为l0 =-1/24*x^3+1/8*x^2-1/12*x ，l1 =1/4*x^3-1/4*x^2-x+1输⼊程序>> L=poly2sym (P),x=0.6; Y = polyval(P,x)运⾏后输出插值多项式和插值为L = Y =x^3+4*x^2-4*x+1 0.2560.输⼊程序>> syms M; x=0.6;R3=M*abs((x-X(1))*(x-X(2)) *(x-X(3)) *(x-X(4)))/24运⾏后输出误差限为R3 =91/2500*M即R 3 )(04.0)4(ξ≈f ， )00.2,00.2(-∈ξ.2.2.5 拉格朗⽇多项式和基函数的MATLAB 程序求拉格朗⽇插值多项式和基函数的MATLAB 主程序function [C, L,L1,l]=lagran1(X,Y)m=length(X); L=ones(m,m);for k=1: mV=1;for i=1:mif k~=iV=conv(V,poly(X(i)))/(X(k)-X(i));endendL1(k,:)=V; l(k,:)=poly2sym (V)endC=Y*L1;L=Y*l例2.2.5 给出节点数据03.17)15.2(=-f ，24.7)00.1(=-f ，05.1)01.0(=f ，03.2)02.1(=f ，06.17)03.2(=f ，05.23)25.3(=f ，作五次拉格朗⽇插值多项式和基函数，并写出估计其误差的公式.解在MATLAB ⼯作窗⼝输⼊程序>> X=[-2.15 -1.00 0.01 1.02 2.03 3.25];Y=[17.03 7.24 1.05 2.03 17.06 23.05];[C, L ,L1,l]= lagran1(X,Y)C =-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =1.0954-4.5745*x+3.3960*x^2+2.1076*x^3+0.0648*x^4-0.2169*x^5L1 =-0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.00040.0331 -0.1377 -0.0503 0.6305 -0.4852 0.0048-0.0693 0.2184 0.3961 -1.2116 -0.3166 1.00330.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097-0.0317 0.0358 0.2530 -0.0426 -0.2257 0.00230.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002l =[ -0.0056*x^5+0.0299*x^4-0.0323*x^3-0.0292*x^2+0.0382*x-0.0004] [ 0.0331*x^5-0.1377*x^4-0.0503*x^3+0.6305*x^2-0.4852*x+0.0048] [ -0.0693*x^5+0.2184*x^4+0.3961*x^3-1.2116*x^2-0.3166*x+1.0033] [ 0.0687*x^5-0.1469*x^4-0.5398*x^3+0.6528*x^2+0.9673*x-0.0097] [ -0.0317*x^5+0.0358*x^4+0.2530*x^3-0.0426*x^2-0.2257*x+0.0023] [ 0.0049*x^5+0.0004 *x^4-0.0266*x^3+0.0001*x^2+0.0220*x-0.0002]估计其误差的公式为)(5x R )25.3)(03.2)(02.1)(01.0()00.1)(15.2(!6)()6(----++=x x x x x x f ξ，)3.25,-2.15(∈ξ.2.2.6 拉格朗⽇插值及其误差估计的MATLAB 程序拉格朗⽇插值及其误差估计的MATLAB 主程序function [y,R]=lagranzi(X,Y,x,M)n=length(X); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:nfor j=1:nif j~=kp=p*(z-X(j))/(X(k)-X(j));endq1=abs(q1*(z-X(j)));c1=c1*j;s=p*Y(k)+s;endy(i)=s;endR=M*q1/c1;例 2.2.6 已知5.030sin = ，1707.045sin = ，0866.060sin = ，⽤拉格朗⽇插值及其误差估计的MATLAB 主程序求40sin 的近似值，并估计其误差.解在MATLAB ⼯作窗⼝输⼊程序>> x=2*pi/9; M=1; X=[pi/6 ,pi/4, pi/3];Y=[0.5,0.7071,0.8660]; [y,R]=lagranzi(X,Y,x,M)运⾏后输出插值y 及其误差限R 为y = R =0.6434 8.8610e-004.2.53 分段插值及其MATLAB 程序2.3.1 ⾼次插值的振荡和MATLAB 程序例2.3.1 作下列函数在指定区间上的n 次拉格朗⽇插值多项式)(x L n)10,8,6,4,2(=n 的图形，并讨论插值多项式)(x L n 的次数与误差)(x R n 的关系.（1）))432sin 3tan(cos()(2x x x f ++=，],[ππ-∈x ；（2）211)(xx f +=，]5,5[-∈x . 解将计算n 次拉格朗⽇插值多项式)(x L n 的值的MATLAB 程序保存名为lagr1.m 的M ⽂件.function y=lagr1(x0,y0,x)n=length(x0); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endend（1）现提供作n 次拉格朗⽇插值多项式)(x L n 的图形的MATLAB 程序，保存名为Li1.m 的M ⽂件m=150; x=-pi:2*pi/(m-1): pi;y=tan(cos((3^(1/2)+sin(2*x))./(3+4*x.^2)));plot(x,y,'k-'),gtext('y=tan(cos((sqrt(3)+sin(2x))/(3+4x^2)))'),pausen=3; x0=-pi:2*pi/(3-1):pi;y1=lagr1(x0,y0,x);hold on,plot(x,y1,'g<'),gtext('n=2'),pause,hold offn=5; x0=-pi:2*pi/(n-1):pi;y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));y2=lagr1(x0,y0,x);hold on,plot(x,y2,'b:'),gtext('n=4'),pause,hold offn=7; x0=-pi:2*pi/(n-1):pi;y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));y3=lagr1(x0,y0,x);hold on,plot(x,y3,'rp'),gtext('n=6'),pause,hold offn=9; x0=-pi:2*pi/(n-1):pi;y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));y4=lagr1(x0,y0,x);hold on,plot(x,y4,'m*'),gtext('n=8'),pause,hold offn=11; x0=-pi:2*pi/(n-1):pi;y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));y5=lagr1(x0,y0,x);hold on,plot(x,y5,'g:'),gtext('n=10')title('⾼次拉格朗⽇插值的振荡')在MA TLAB ⼯作窗⼝输⼊名为 Li1.m 的M ⽂件的⽂件名>> Li1.m回车运⾏后，便会逐次画出)(x f 在区间],[ππ-上的n 次拉格朗⽇插值多项式)(x L n )10,8,6,4,2(=n 的图形(略).（2）在MATLAB ⼯作窗⼝输⼊程序m=101; x=-5:10/(m-1):5; y=1./(1+x.^2);z=0*x;plot(x,z,'r',x,y,'k-'),gtext('y=1/(1+x^2)'),pausey1=lagr1(x0,y0,x);hold on,plot(x,y1,'g'),gtext('n=2'),pause,hold offn=5; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);y2=lagr1(x0,y0,x);hold on,plot(x,y2,'b:'),gtext('n=4'),pause,hold offn=7; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);y3=lagr1(x0,y0,x);hold on,plot(x,y3,'r'),gtext('n=6'),pause,hold offn=9; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);y4=lagr1(x0,y0,x);hold on,plot(x,y4,'r:'),gtext('n=8'),pause,hold offn=11; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);y5=lagr1(x0,y0,x);hold on,plot(x,y5,'m'),gtext('n=10')title('⾼次拉格朗⽇插值的振荡')回车运⾏后，便会逐次画出)1/(12x y +=在区间]5,5[-上的n 次拉格朗⽇插值多项式)(x L n )10,8,6,4,2(=n 的图形(略).2.3.4 作有关分段线性插值图形的MATLAB 程序作有关分段线性插值图形的MATLAB 主程序function s= xxczhjt1(x0,y0,xi,x,y)s= interp1(x0,y0,xi);Sn= interp1(x0,y0,x0);plot(x0,y0,'o',x0,Sn,'-',xi,s,'*',x,y,'-.')legend('节点(xi,yi)','分段线性插值函数Sn(x)','插值点(x,s)','被插值函数y')我们也可以直接在在MA TLAB ⼯作窗⼝编程序.例如，>> 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 及其分段线性插值函数和节点的图形')>> x0 =-6:6; y0 =cos(x0);xi = -6:.25:6;x=-6:0.001:6; y=cos(x); plot(x0,y0,'o',xi,yi,x,y,':'),legend('节点(xi,yi)','分段线性插值函数','被插值函数y=cosx ')title('y=cosx 及其分段线性插值函数和节点的图形')例2.3.5 设函数22511)(xx f +=，在区间]1,1[-上取等距节点),(i i y x 10,,2,1,0 =i , 构造分段线性插值函数)(x S n ，⽤MA TLAB 程序计算各⼩区间的中点i x 处)(x S n 的值，作出节点，插值点，)(x f 和)(x S n 的图形.解节点的横坐标和插值点等取值与例6.5.4相同.在MATLAB ⼯作窗⼝输⼊程序>>x0=-1:0.2:1; y0=1./(1+25.*x0.^2);xi=-0.9:0.2:0.9;b=max(x0);a=min(x0);x=a:0.001:b;y=1./(1+25.*x.^2);s=xxczhjt1(x0,y0,xi,x,y), title('y=1/(1+25 x^2)的分段线性插值的有关图形')运⾏后屏幕显⽰各⼩区间中点i x 处)(x S n 的值，出现节点，插值点，)(x f 和)(x S n 的图形(略)s =Columns 1 through 40.04864253393665 0.07941176470588 0.15000000000000 0.35000000000000Columns 5 through 80.75000000000000 0.75000000000000 0.35000000000000 0.15000000000000Columns 9 through 100.07941176470588 0.048642533936652.3.5 ⽤MATLAB 计算有关分段线性插值的误差例2.3.8 设函数))432sin 3tan(cos()(2xx x f ++=,在区间],[ππ-上取等距节点),(i i y x ，9,,2,1,0 =i ，构造分段线性插值函数)(x S n . (1)⽤MA TLAB 程序计算各⼩区间中点i x 处)(x S n 的值，作出节点，插值点，)(x f 和)(x S n 的图形；(2) ⽤MA TLAB 程序计算各⼩区间中点处)(x S n 的值及其相对误差；(3) ⽤MA TLAB 程序估计)(max "x f x ππ≤≤-和)(x S n 在区间],[ππ-上的误差限. 解在MATLAB ⼯作窗⼝输⼊程序>>h=2*pi/9; x0=-pi:h:pi;y0=tan(cos((sqrt(3)+sin(2*x0))./(3+4*x0.^2)));xi=-pi+h/2:h:pi-h/2;b=max(x0); a=min(x0);x=a:0.001:b; y=tan(cos((sqrt(3)+sin(2.*x))./(3+4*x.^2))); si=xxczhjt1(x0,y0,xi,x,y);xi,fi,si,Ri,i=[xi',fi',si',Ri']title('y=tan(cos((sqrt(3)+sin(2 x))/(3+4 x^2)))的分段线性插值的有关图形')运⾏后屏幕显⽰R i （略），并且作出节点，插值点，)(x f 和)(x S n 的图形(略).在MATLAB ⼯作窗⼝输⼊程序>> syms xy=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));yxx=diff(y,x,2),运⾏后屏幕显⽰函数)(x f 的⼆阶导数)(''x f （略）.在MATLAB ⼯作窗⼝输⼊程序>> h=2*pi/9; x=-pi:0.0001:-pi;yxx=2*tan(cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1/2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3^(1/2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2*cos(2.*x)./(3+4*x.^2)-8*(3^(1/2)+sin (2.*x))./(3+4.*x.^2).^2.*x).^2-(1+tan(cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2)).*(2.*c os(2.*x)./(3+4.*x.^2)-8*(3^(1/2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2-(1+tan(cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1/2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32*cos(2.*x)./(3+4.*x.^2).^2.*x+128*(3^(1/2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8*(3^(1/2)+sin(2.*x))./(3+4.*x.^2).^2);myxx=max(yxx), R=(h^2)* myxx/8运⾏后屏幕显⽰myxx =)(max "x f x π≤≤π-和)(x S n 在区间],[ππ-上的误差限R 如下 myxx = R =-0.02788637150664 -0.001698934907112.4 分段埃尔⽶特插值及其MATLAB 程序2.4.2 分段埃尔⽶特插值的MATLAB 程序调⽤格式⼀：YI=interp1(X,Y ,XI,'pchip')调⽤格式⼆：YI=interp1(X,XI,'pchip')例2.4.5 试⽤MA TLAB 程序计算例6.6.1中在各⼩区间中点处分段三次埃尔⽶特插值)(2/1+i n x H 及其相对误差.解在MATLAB ⼯作窗⼝输⼊程序>> h=0.2;x0=-1:h:1;y0=1./(1+25.*x0.^2); xi=-0.9:h:0.9;fi=1./(1+25.*xi.^2); yi=interp1(x0,y0,xi,'pchip');Ri=abs((fi-yi)./fi); xi,fi,yi,Ri,i=[xi',fi',yi',Ri'] 运⾏后屏幕显⽰各⼩区间中点x i 处的函数值f i ，插值s i ，相对误差值R i (略).2.4.3 作有关分段埃尔⽶特插值图形的MATLAB 程序作有关分段埃尔⽶特插值图形的MATLAB 主程序function H=hermitetx(x0,y0,xi,x,y)H= interp1(x0,y0,xi,'pchip');Hn= interp1(x0,y0,x,'pchip');plot(x0,y0,'o',x,Hn,'-',xi,H,'*',x,y,'-.')legend('节点(xi,yi)', '分段埃尔⽶特插值函数','插值点(x,H)','被插值函数y')我们也可以直接在在MATLAB ⼯作窗⼝编程序，例如，>> x0 =-6:6; y0 =sin(x0); xi = -6:.25:6;'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 及其分段埃尔⽶特插值函数和节点的图形')>> x0 =-6:6; y0 =cos(x0);xi = -6:.25:6;yi = interp1(x0,y0,xi,'pchip');x=-6:0.001:6; y=cos(x); plot(x0,y0,'o',xi,yi,x,y,':'),legend('节点(xi,yi)','分段埃尔⽶特插值函数','被插值函数y=cosx')title(' y=cosx 及其分段埃尔⽶特插值函数和节点的图形')例2.4.6 设函数211)(xx f +=定义在区间]5,5[-上，节点(X (i ),f (X (i )))的横坐标向量X 的元素是⾸项a =-5，末项b =5，公差h =1.5的等差数列，构造三次分段埃尔⽶特插值函数)(3,x H n .把区间]5,5[-分成20等份，构成20个⼩区间，⽤MA TLAB 程序计算各⼩区间中点i x 处)(3,x H n 的值，并作出节点，插值点，)(x f 和)(3,x H n 的图形.解在MATLAB ⼯作窗⼝输⼊程序>>x0=-5:1.5:5;y0=1./(1+x0.^2); x1=-4.75:0.5:4.75;x=-5:0.001:5; y=1./(1+x.^2); H= hermitetx(x0,y0,x1,x,y)title('函数y=1/(1+x^2)及其分段埃尔⽶特插值函数，插值，节点(xi,yi)的图形')运⾏后屏幕显⽰各⼩区间中点i x 处)(3,x H n 的值，出现节点，插值点，)(x f 和)(3,x H n 的图形(略).例2.4.7 设函数x x x f cos 5.0)(-=定义在区间],[ππ-上，取7=n ，按等距节点构造分段埃尔⽶特插值函数)(3,7x H ，⽤MATLAB 程序计算各⼩区间中点i x 处)(3,7x H 的值，作出节点，插值点，)(x f 和)(3,7x H 的图形.解记节点的横坐标7,,2,1,0,7/2, =π=+π-=i h ih x i 插值点)(21121+++=i i i x x x ，6,,2,1,0 =i .在MATLAB ⼯作窗⼝输⼊程序 >> h=2*pi/7; x0=-pi:h:pi;y0=0.5.*x0-cos(x0); xi=-pi+h/2:h:pi-h/2;b=max(x0); a=min(x0); x=a:0.001:b;y=0.5.*x-cos(x); H= hermitetx(x0,y0,xi,x,y)title('函数y=0.5x-cos(x)及其分段埃尔⽶特插值函数，插值，节点(xi,yi) 的图形')运⾏后屏幕显⽰各⼩区间中点i x 处)(3,7x H 的值，出现节点，插值点，)(x f 和)(3,7x H 的图形(略).2.4.4 ⽤MATLAB 计算有关分段埃尔⽶特插值的误差例2.4.8 设函数22511)(xx f +=定义在区间]1,1[-上，取10=n ，按等距节点构造分段埃尔⽶特插值函数)(3,x H n ，⽤MATLAB 程序在]1,1[-上计算)(max )4(11x f x ≤≤-和)(3,x H n 的误差公式和误差限.解在MATLAB ⼯作窗⼝输⼊程序yxxxx=150000000./(1+25.*x.^2).^5.*x.^4-4500000./(1+25.*x.^2).^4.*x.^2+15000./(1+25.*x.^2).^3;myxxxx=max(yxxxx), R=(h^4)* abs(myxxxx/384)运⾏后输出)(x f 的4阶导数在区间]1,1[-上绝对值的最⼤值myxxxx 和)(3,x H n 在区间]1,1[-上的误差公式myxxxx 为myxxxx = R =15000 625/16*h^4（4）在MATLAB ⼯作窗⼝输⼊程序>> h=0.2; R =625/16*h^4运⾏后输出误差限为R =0.06250000000000例2.4.9 设函数))432sin 3tan(cos()(2xx x f ++=定义在区间],[ππ-上，取9=n ，按等距节点构造分段埃尔⽶特插值函数)(3,x H n . (1)⽤MA TLAB 程序计算各⼩区间中点i x 处)(3,x H n 的值，作出节点，插值点，)(x f 和)(3,x H n 的图形；(2) 并⽤MATLAB 程序计算各⼩区间中点处)(3,x H n 的值及其相对误差；(3) ⽤MA TLAB 程序求)(max )4(x f x ππ≤≤-和)(3,x H n 在区间],[ππ-上的误差公式和各插值的误差限.解（1）记节点的横坐标9,,2,1,0,9/2, =π=+π-=i h ih x i ，插值点)(21121+++=i i i x x x ，8,,2,1,0 =i .在MATLAB ⼯作窗⼝输⼊程序>>h=2*pi/9; x0=-pi:h:pi;y0=tan(cos((sqrt(3)+sin(2*x0))./(3+4*x0.^2)));xi=-pi+h/2:h:pi-h/2;fi=tan(cos((3^(1/2)+sin(2*xi))./(3+4*xi.^2)));b=max(x0); a=min(x0); x=a:0.001:b;y=tan(cos((3^(1/2)+sin(2.*x))./(3+4*x.^2)));Hi= hermite tx (x0,y0,xi,x,y);Ri=abs((fi-Hi)./fi); xi,fi,Hi,Ri,i=[xi',fi',Hi',Ri']title('函数y=tan(cos((sqrt(3)+sin(2x))/(3+4x^2)))及其分段埃尔⽶特插值函数，插值，节点(xi,yi) 的图形')运⾏后屏幕显⽰各⼩区间中点x i 处的函数值f i ，插值H i ，相对误差值R i ，并且作出节点，插值点，)(x f 和)(3,x H n 的图形(略).（2）在MATLAB ⼯作窗⼝输⼊程序>> syms xy=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));yxxxx=diff(y,x,4),%simple(yxxxx)运⾏后屏幕显⽰函数)(x f 的4阶导数)()4(x f，然后将输出的)()4(x f 编程求)(max )4(x f x ππ≤≤-和)(3,x H n 及其在区间],[ππ-上的误差限的MA TLAB 程序如下>>syms h,x=-pi:0.0001:pi;yxxxx=-12.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).^2.*sin((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)).^3.*(2.*c os(2.*x)./(3.+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./(3.+4.*x.^2)-32.*cos(2.*x)./(3.+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^3.*x.^2.-8../2)+sin(2.*x))./(3.+4.*x.^2))).^2).^2.*sin((3.^(1./2)+sin(2.*x ))./(3.+4.*x.^2)).^4.*(2.*cos(2.*x)./(3.+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4.*tan(cos((3.^(1./2)+sin(2.* x))./(3.+4.*x.^2)))+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x. ^2))).^3.* (1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2) .*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^4.*(2.*cos(2.*x)./(3 +4.*x.^2)-8.* (3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4-8.*t an(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x)). /(3.+4.*x.^2)).^2.*(2.*cos(2.*x)./(3.+4.*x.^2)-8*(3.^(1./2)+sin(2.*x))./(3.+4*x.^2).^2.*x).^4+6.*(1+tan(cos((3.^(1./2)+sin(2.* x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)) .* (2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x. ^2).^2.*x).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8 .*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4-3.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^( 1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2).^2-4.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x) .*(-8.*cos(2.*x)./(3+4.*x.^2)+96.*sin(2.*x)./(3+4.*x.^2).^2.*x+768.*cos(2.*x)./(3+4.*x.^2).^3.*x.^2-48.*cos(2.*x)./(3+4.*x.^2) .^2-3072.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x)-(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x. ^2)).*(16.*sin(2.*x)./(3+4.*x.^2)+256.*cos(2.*x)./(3+4.*x.^2).^ 2.*x-3072.*sin(2.*x)./(3+4.*x.^2).^3.*x.^2+192.*sin(2.*x)./(3+4 .*x.^2).^2-24576.*cos(2.*x)./(3+4.*x.^2).^4.*x.^3+3072.*cos(2.* x)./(3+4.*x.^2).^3.*x+98304.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2) .^5.*x.^4-18432.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^2+38 4.* (3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3)-12.*(1+tan(cos((3.^(1 ./2)+sin(2.*x))./(3+4.*x.^2))).^2).^2*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin (2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+sin(2.*x))./(3+4 .*x.^2))-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*( 1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^3.*(2.*cos(2.*x)./(3+4.*x.^2)-8.* (3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin (2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*( 1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.* (3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))+36.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x .^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32. *cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+6.*ta n(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3x))./(3+4.*x.^2).^2.*x).^4+6.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+si n(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2).^2+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).*(-8.*cos(2.*x)./(3+4.*x.^2)+96.*sin(2.*x)./(3+4.*x.^2).^2.*x+768.*cos(2.*x)./(3+4.*x.^2).^3.*x.^2-48.*cos(2.*x)./(3+4.*x.^2).^2-3072.*(3.^(1./2)+s in(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x)myxx=max(yxxxx), R=(h.^4).* abs(myxx./384)将其保存为名为myxx.m 的M ⽂件，然后在MATLAB ⼯作窗⼝输⼊该⽂件名>> myxx运⾏后屏幕显⽰myxx =)(max )4(x f x π≤≤π-和)(3,x H n 在区间],[ππ-上的误差公式)(ma x 384)4(4x f h R x π≤≤π-=如下myxx = R =73.94706841647552 1734520780029061/9007199254740992*h^4最后在MATLAB ⼯作窗⼝输⼊>>h=2*pi/9; R =1734520780029061/9007199254740992*h^4运⾏后屏幕显⽰)(3,x H n 在区间],[ππ-上的误差限R =0.045744530299482.5 三次样条（Spline ）及其MATLAB 程序2.5.4 ⽤⼀阶导数计算的⼏种样条函数和MATLAB 程序计算压紧样条的MATLAB 主程序function[m,H,lambda,mu,D,A,dY,sk]=ClampedSp (X,Y,dyo,dyn)m=length(X);A=zeros(m,m);n=m-1;H=zeros(1,n);lambda=zeros(1,n);mu=zeros(1,n);lambda(1)=0;A(1,1)=2;A(1,2)=lambda(1);D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1-lambda(1);D(1)=2*dyo;for k=1:nhk=X(k+1)-X(k);H(k+1)=hk;endH=H(2:n+1);for k=1:n-1lambdak=H(k)/(H(k)+H(k+1)); lambda(k+1)=lambdak;muk=1-lambda(k+1);mu(k)=muk;dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)-Y(k+1))./H(k+1)));D(k+1)=dk;endD(m)=2*dyn;mu(n)=0; n;H;lambda;mu;D;for i=1:m-1A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i);enddY=A\D';m=length(X);S=zeros(m-1,1);for k=2:msk=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+Y(k)*((H(k-1)+2*X(k)-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H (k-1)^2)end例2.5.2 ⽤计算压紧样条的MATLAB 主程序ClampedSp.m 求例6.7.1的问题. 解保存ClampedSp.m 为M ⽂件，输⼊程序>> X=[0:2:6,10],Y=[0,16,36,54,82],dyo=8,dyn=7,[m,H,lambda,mu,D,A,dY,sk]=ClampedSp(X,Y,dyo,dyn)运⾏后输出压紧样条函数sk ，X 的维数m ，由1-i h ),,2,1(n i =、i λ、i µ和i d)1,,2,1(-=n i 的坐标组成的向量H 、lambda 和mu 、D ，（6.90）式的系数矩阵A ，线性⽅程组（6.90）的解向量dY 如下sk =2*(6-2*x)*x^2+2*x*(x-2)^2+9/4*(x-2)*x^2sk =2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+9/4*(x-2)*(x-4)^2+5/2*(x-4)*(x-2)^2sk =9/2*(-6+2*x)*(x-6)^2+27/4*(14-2*x)*(x-4)^2+5/2*(x-4)*(x-6)^2+2*(x-6)*(x-4)^2sk =27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+1/2*(x-6)*(x-10)^2+7/16*(x-10)*(x-6)^2m =5H =2 2 2 4lambda =0 0.5000 0.5000 0.3333mu =0.5000 0.5000 0.6667 0D =16.0000 27.0000 28.5000 25.0000 14.0000A =2.0000 0 0 0 00.5000 2.0000 0.5000 0 00 0.5000 2.0000 0.5000 00 0 0.6667 2.0000 0.33330 0 0 0 2.0000dY =8.00009.000010.00008.00007.0000输⼊程序>> x2=3,s2=2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+9/4*(x-2)*(x-4)^2+5/2*(x-4)*(x-2)^2x4=8,s4=27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+1/2*(x -6)*(x-10)^2+7/16*(x-10)*(x-6)^2运⾏后输出)3()3(2s f ≈和)8()8(4s f ≈的值如下3 25.7500 8 68.5000例2.5.2和例2.5.1计算的结果完全相同.计算⾃然样条的MATLAB 主程序function [m,H,lambda,mu,D,A,dY,sk]=NaturalSp(X,Y)m=length(X);A=zeros(m,m);n=m-1;H=zeros(1,n);lambda=zeros(1,n);mu=zeros(1,n);lambda(1)=1;A(1,1)=2;A(1,2)=lambda(1); D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1;D(1)=3*(Y(2)-Y(1)); for k=1:nhk=X(k+1)-X(k);H(k+1)=hk;endH=H(2:n+1);for k=1:n-1lambdak=H(k)/(H(k)+H(k+1)); lambda(k+1)=lambdak;muk=1-lambda(k+1);mu(k)=muk;dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)-Y(k+1))./H(k+1)));D(k+1)=dk;endD(m)= 3*(Y(m)-Y(m-1))/H(m-1);mu(n)=1; n;H;lambda;mu;D;for i=1:m-1A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i);enddY=A\D';syms xm=length(X);S=zeros(m-1,1);for k=2:msk=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+Y(k)*((H(k-1)+2*X(k)-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H (k-1)^2)end例2.5.3 ⽤计算⾃然样条的MATLAB 主程序NaturalSp.m 求例6.7.1的问题.解保存ClampedSp.m 为M ⽂件，输⼊程序>> X=[0:2:6,10],Y=[0,16,36,54,82],dyo=8,dyn=7,[m,H,lambda,mu,D,A,dY,sk]=NaturalSp(X,Y)运⾏后输出⾃然样条函数sk ， X 的维数m ，由1-i h ),,2,1(n i =、i λ、i µ和i d)1,,2,1(-=n i 的坐标组成的向量H 、lambda 和mu 、D ，（6.88）式的系数矩阵A ，线性⽅程组（6.88）的解向量dY 如下sk =2*(6-2*x)*x^2+915/172*x*(x-2)^2+117/86*(x-2)*x^2sk =2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+117/86*(x-2)*(x-4)^2+471/172*(x-4)*(x-2)^2sk =9/2*(-6+2*x)*(x-6)^2+27/4*(14-2*x)*(x-4)^2+471/172*(x-4)*(x-6)^2+333/172*(x-6)*(x-4)^2sk =27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+333/688*(x-6)*(x-10)^2+285/688*(x-10)*(x-6)^2m = 5H = 2 2 2 4lambda = 1 1/2 1/2 1/3mu = 1/2 1/2 2/3 1D = 48 27 57/2 25 21A =1/2 2 1/2 0 00 1/2 2 1/2 00 0 2/3 2 1/30 0 0 1 2dY =915/43234/43471/43333/43285/43输⼊程序>> x2=3,s2=2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+117/86*(x-2)*(x-4)^2+471/172*(x-4)*(x-2)^2x4=8,s4=s4=27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+333/688*(x-6)*(x-10)^2+285/688*(x-10)*(x-6)^2运⾏后输出)3()3(2s f ≈和)8()8(4s f ≈的值如下x2 = s2 = x 4= s4 =3 24.6221 8 68.5581例2.5.3和例2.5.2⽤的⽅法不同，所以计算的结果不同，但是两种⽅法计算的结果相近.2.5.6 ⽤MATLAB 计算三次样条分别求在下列条件下在插值点1，2处的压紧三次样条插值，并显⽰该样条函数的有关信息：（1）端点约束条件为5)1(=-'n S ，16.29)50.3(='n S ；（2）端点约束条件为0)1(=-'nS ，0)50.3(='n S . 解（1）输⼊MATLAB 程序>> X=[-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.823.50];Y=[-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31];XI=[-0.02 2.56];YI= spline (X, [5,Y,29.16],XI), PP = spline (X, [5,Y,29.16])运⾏后屏幕显⽰压紧样条分别在02.01-=x ，56.22=x 处的插值和该样条函数的有关信息如下YI =-3.1058 4.7834PP =form: 'pp'breaks: [-1 -0.5400 0.1300 1.1200 1.8900 2.0600 2.5400 2.8200 3.5000]coefs: [8x4 double]pieces: 8order: 4dim: 1（2）因为端点约束条件为0)1('=-n S ，0)50.3('=n S ，所以输⼊MATLAB 程序>> YI= spline (X, [0,Y,0],XI), PP= spline (X, [0,Y,0])运⾏后屏幕显⽰压紧三次样条分别在02.01-=x ，56.22=x 的插值和该样条函数的有关信息如下YI =-3.0192 4.7501PP =form: 'pp'breaks: [-1 -0.5400 0.1300 1.1200 1.8900 2.0600 2.5400 2.8200 3.5000]coefs: [8x4 double]pieces: 8order: 4dim: 1例 2.5.7 在默认端点约束条件下，⽤两种⽅法计算在下列插值点处的三次样条插值.（1）给出节点的数据与例2.5.6相同,插值点XI =2.56；（2）节点(X (i ),Y (i ))的横坐标向量X 从5-到5的整数，纵坐标向量Y =（-2.36，0.85，1.12，2.36，2.36，3，4，1.46，0.49，0.06, 0），插值点向量XI 是从4-到4的11个等分点.解（1）输⼊MATLAB 程序>>X=[-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.823.50];Y=[-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31];XI=2.56;Y1= spline (X, Y,XI),Y2=interp1(X,Y,XI,'spline')运⾏后屏幕显⽰三次样条插值的两种结果如下Y1 = 4.730 2，Y2 =4.730 2.（2）输⼊MATLAB 程序>> X= -5:5; Y= [-2.36 .85 1.12 2.36 2.36 3 4 1.46 .49 .06 0]; XI = linspace(-4,4,11); Y1= spline (X, Y,XI),Y2=interp1(X,Y,XI,'spline ')运⾏后屏幕显⽰Y1 =0.8500 1.0203 1.8857 2.4779 2.3683 3.00004.0656 2.5935 0.8247 0.4074 0.0600Y2 =0.8500 1.0203 1.8857 2.4779 2.3683 3.00004.0656 2.5935 0.8247 0.4074 0.0600例 2.5.8 设函数22511)(x x f +=定义在区间]1,1[-上，取10=n ，按等距节点构造分段三次样条函数)(x S n .试⽤MA TLAB 程序计算在各⼩区间中点处分段三次样条插值)(2/1+i n x S 及其相对误差.解 .在MATLAB ⼯作窗⼝输⼊程序>> h=0.2;x0=-1:h:1;y0=1./(1+25.*x0.^2);xi=-0.9:h:0.9;fi=1./(1+25.*xi.^2); yi=spline (x0,y0,xi);Ri=abs((fi-yi)./fi); xi,fi,yi,Ri,i=[xi',fi',yi',Ri']运⾏后屏幕显⽰各⼩区间中点x i 处的函数值f i ，插值s i ，相对误差值R i (略).2.5.7 ⼏种作三次样条有关图像的MATLAB程序求有关分段三次样条图形的MATLAB主程序（⼀）限定端点约束条件的作图程序function S=splinetx(x0,y0,xｊ,x,y,dy1,dyn)S = spline(x0,[dy1,y0,dyn],xｊ);Sn = spline(x0,[dy1,y0,dyn],x);plot(x0,y0,'o',x,Sn,'-',xｊ,S,'*',x,y,'-.')legend('节点(xi,yi)', '分段三次样条函数','插值点(x,S)','被插值函数y')（⼆）不限定端点约束条件的作图程序function S=splinetx1(x0,y0,xi,x,y)S= interp1(x0,y0,xi, 'spline'); Sn= interp1(x0,y0,x,'spline');plot(x0,y0,'o',x,Sn,'-',xi,S,'*',x,y,'-.')legend('节点(xi,yi)', '分段三次样条函数','插值点(x,S)','被插值函数y')。
计算方法-第2章-1、插值法(拉格朗日插值)
2019/1/15
26
证明：假设在区间[a,b]上f(x)的插值多项式为 Ln ( x) 令
Rn ( x) f ( x) Ln ( x)
显然在插值节点为 xi (i 0,1,, n)上 Rn ( xi ) f ( xi ) Ln ( xi ) 0 , i 0,1,, n 因此Rn ( x)在[a, b]上至少有n 1个零点
(k 0,1,2,, n)
且
n1 ( x) Ln ( x) yk ' ( x x ) k 0 k n 1 ( xk )
n
2019/1/15
18
总 结
于是, y f ( x)在节点xi (i 0 ,1, , n)上, 以l j ( x) (i 0 ,1, , n) 为插值基函数的插值多 项式(记为Ln ( x))为
本章只讨论多项式插值与分段插值
2019/1/15 7
§ 2.2
拉格朗日插值
• 此插值问题可表述为如下： • 问题 求作次数 n 多项式 Ln ( x) ，使满足条件
Ln x yi , (i 0,1,, n)
• 这就是所谓的拉格朗日（Lagrange）插值。
2019/1/15
8
§ 2.2.1
线性插值的局限性
2019/1/15
12
三、抛物插值
问题 求作二次式 L2 ( x) ，使满足条件
L2 ( x j ) y j
( j k 1, k , k 1)
二次插值的几何解释是用通过三个点
的抛物线来近似考察曲线，故称为拋物插值。类似于线性 插值，构造基函数，要求满足下式：
L2(x) yk 1lk 1 ( x) yklk ( x) yk 1lk 1 ( x)
2插值法
18
Ex：已知特殊角300，450，600 的正弦值，分别用一 次插值，二次插值函数近似sin500 的值。
19
3 代数多项式插值的存在唯一性 §2. 2.3
线性插值和二次插值都属于代数多项式插值。 对于一般的代数插值问题,就是寻求一个不高于n次 的代数多项式 Pn(x)=a0+a1x+a2x2+ … +anxn 使其在 给定的n+1个互异的插值基点上满足插值原则 …,n Pn(xi)=yi, i=0,1, i=0,1,…
� �
23
我们称Rn(x)为插值多项式Pn(x)的余项。 显然有 …,n Rn(xi)=0, i=0,1,2, i=0,1,2,… 下面给出插值多项式Pn(x)余项的表达式。 定理：设函数 f(x) 在区间［ a,b ］上具有 n+1 阶导 数， Pn(x)为次数不高于n的多项式,且 Pn(x0)=y0 Pn(x1)=y1 … Pn(xn)=yn
28
利用公式可以给出用多项式 Pn(x) 近似代替 f(x)的误差估计。这里还得说明几点: (1) 插值多项式本身只与插值基点及 f(x)在这些基 点上的函数值有关,而与函数f(x)本身并没有关系。 但余项Rn(x)却与f(x)联系很紧。
�
29
�
(2)若f(x)为次数不超过 n 的多项式 , 那么以 n+1个点为 ≡f(x)。这 基点的插值多项式就一定是其本身, 即Pn(x) (x)≡ 是因为此时 Rn(x)=0。 (3) 从余项 Rn(x) 中的 ω (n+1)(x) 知 , 当点 x 位于 x0,x1, … ,xn ω n+1(x)|比较小 , 精度要高一些, 而位于两端 的中部时 ,| ,|ω 时,精度要差一些;若x位于x0,x1,…,xn的外部,一般称外 插(或外推),此时精度一般不理想,使用时须注意。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第二章函数的插值学习目标：掌握多项式插值的Lagrange插值公式、牛顿插值公式等，等距节点插值、差分、差商、重节点差商与埃米特插值。
重点是多项式插值方法。
§2.1 多项式插值一、插值问题的基本概念：设有函数 ，只知道它在n+1个不同的点上的函数值，y 是另外一点。
不知道，如何求它的近似值？插值就是一种办法，它的思想是：找一个简单的已知解析表达式的函数 ，使得（1）并且 容易计算，我们就用 来代替。
)(x f 121,,,+n x x x )(y f )(x p ,1,,2,1),()(+==n i x f x p i i )(y p )(y p )(y f 称为插值函数， 称为被插值函数， 称为节点或结点。
如果限制 为n 次多项式，那么上述问题称为多项式插值， 称为 的n 次插值多项式。
)(x p )(x f 121,,,+n x x x )(x p )(x p )(x f 本节主要介绍插值问题的基本概念、方法和理论。
对于多项式插值，我们主要讨论以下几个问题： 插值问题是否可解，如果有解，解是否唯一？ 插值多项式的常用构造方法有哪些？插值函数逼近的误差如何估计，即截断误差的估计； 当插值节点无限加密时，插值函数是否收敛于 。
本节主要讨论前三个问题。
二、多项式插值的方法令是n 次多项式：(1)考察（1）式，就有方程组： (2)nn nk k k x a x a a x a x p +++==∑= 100)(1,,2,1,0+==∑=n i f x a i nk ki k )(x p )(x f刚巧是一个具有n+1个未知量 ,n+1个方程的线性方程组，它的系数矩阵为：n a a a ,,,10 ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=+++n n n n n nx x x x x x x x x A 121122221211111 而 恰为范德蒙达(Vandermonde) 行列式。
由线性代数知：)det(A 0)()det(11≠-=∏+≤<≤n j i i jx xA 因而方程组（2）有唯一解存在，也即由插值条件（1）可以唯一确定一个n 次多项式。
n a a a ,,,10定理1 由n+1个不同的点可以唯一确定一个n 次多项式 ，且满足条件（1）。
121,,,+n x x x )(x p 注意：定理1解决了n 次插值多项式的存在、唯一性。
以下我们主要讨论 的具体求法： )(x p （3）首先，从方程组（2），由克莱姆(Cram)法则我们知道：)det(A d a k k =其中是将系数矩阵A 的第k 列换为方程组（2）的右端向量形成的矩阵行列式 。
k d 从（3）式求，从而求得 ，从数学上来说是很清楚的，但从计算来看工作量太大了。
因为计算一个行列式是不容易的，花的代价较大。
以下我们介绍常用的求的方法 k a )(x p )(x p§1.1 拉格朗日途径考虑特殊的n 次多项式:)())(())(()())(())(()(1112111121++-++-----------=n k k k k k k k n k k k x x x x x x x x x x x x x x x x x x x x x q 它满足：⎩⎨⎧≠==ki ki x q i k ,0,1)(则： 若记 ∏+=-=11)()(n i i x x x w )()()()(k k k x w x x x w x q '-=由此看出多项式 ： ∑+=11)(n k k k x q f 是一个n 次多项式，在 处恰为 ix if ，因此满足条件（1），由定理1知的这种表示称为拉格朗日（Lagrange ）插值多项式，其中， (k=1,2,…,n+1)称为拉格朗日基本插值多项式。
∑+==11)()(n k k k x q f x P （4）)(x q k 例15)(,1)(,8)(,4,2,1321321======x f x f x f x x x 求二次插值多项式。
)(x P 211635)24)(14()2)(1(1)42)(12()4)(1(8)41)(21()4)(2()(2+-=----+----+----=x x x x x x x x x P 解 按拉格朗日方法，有：§1.2 内维尔途径内维尔途径是一种由两个n-1次插值多项式构造一个n 次插值多项式的方法。
记 是 在 上的一次插值多项式，同样 是 在 上的一次插值多项式，那么 上的二次插值多项式为)(2,1x P )(x f 21,x x )(3,2x P )(x f 32,x x 321,,x x x )()(1)()()(3,232,11133,21312,11333,2,1x P x x x P x x x x x P x x x x x P x x x x x P ---=--+--=这是很容易证明的。
因为 显然是二次多项式，且 ，由定理1就知道 是唯一的 在上的二次插值多项式。
)(3,2,1x P )(3,2,1x P )3,2,1)(()(3,2,1==i x f x P ii )(x f 321,,x x x一般地，若记 是 在 上的n-1次插值多项式， 是 在 上的n-1次插值多项式，那么 在 上 的n 次插值多项式为)(,,2,1x P n)(x f n x x x ,,,21 )(1,,,2x P n n + )(x f 121,,,,+n n x x x x )(x f 121,,,,+n n x x x x )()(1)()()(1..,3,21,,2,11111..,3,2111,,2,11111,,3,2,1x P x x x P x x x x x P x x x x x P x x x x x P n n n n n n n n n n n n ++++++++---=--+--= （5）这也是很容易验证的。
（5）有一个很重要的特点，就是，如果 ，那么，这意味着 是 和 的带权平均，且权系数是正的，这样得到的传播误差不会超过 和 两个误差中大的那一个，这在计算上是有好处11+≤≤n x x x ,0111≥--++x x x x n n 0111≥--+x x x x n )(1,,,2,1x P n n + )(,,2,1x P n )(1,,,2x P n n + )(,,2,1x P n )(1,,,2x P n n +例如，已知， 用内维尔方法求 的计算过程列表如下64)(,8)(,0)(,8)(,8,6,4,243214321===-=====x f x f x f x f x x x x )5(f ix )(i x f ix -54183241=--44143261=--12343281=---48101461=--220341481-=---2064381681-=----364 8 -18 6 10 4 3-8 2因此 1)5()5(=≈p f§1.3 牛顿途径对于n+1个不同的节点，考虑n 次多项式 121,,,+n x x x )())(())(()()(21212110n n x x x x x x c x x x x c x x c c x Q ---++--+-+=（6） 如果满足： ，那么它就是n+1个点上的n 次插值多项式，对于这样的，有 1,,,3,2,1,)()(+===n n i f x f x Q i i i )(x Q ⎪⎪⎪⎩⎪⎪⎪⎨⎧=---++-+==---++-+==-+===++++++--112111111011211110212102101)())(()()()())(()()()()()(n n n n n n n n n n n n n n n n f x x x x x x c x x c c x Q f x x x x x x c x x c c x Q f x x c c x Q f c x Q从可以求出 ，再从 可以求出 ，依次从 可以求出 ，最后从可以求出 。
)(1x Q 0c )(2x Q 1-n c )(1+n x Q 1c )(1+n x Q n c 例如，已知5)(,1)(,8)(,4,2,1321321======x f x f x f x x x 从，得 ， 101)(f c x Q ==80=c 从，得 212102)()(f x x c c x Q =-+=7)/()(12021-=--=x x c f c 再从 ，得，即： 32313213103))(()()(f x x x x c x x c c x Q =--+-+=32=c )2)(1(3)1(78)(--+--=x x x x Q 这样逐步获得，计算很方便。
n c c c ,,,10从的表达式（6）知道它的首项系数为 ，从上述计算过程还知道， 恰是两点 上一次插值的首项系数， 恰是3点上二次插值的首项系数，依次 恰是个i +1点 上i 次插值多项式 的首项系数。
的值取决于 )(x Q n c 2c 21,x x 1c 321,,x x x i c 121,,,+i x x x i c ,)),(,()),(,(2211 x f x x f x ))(,()),(,(11++i i i i x f x x f x 这i+1个平面上的点。
定义1在 上有定义，在其上确定的k 次多项式的首项系数称为在 上的k 阶差商，记为 )(x f 121,,,,+k k y y y y )(x f 121,,,,+k k y y y y ),,,,(121+k k y y y y f 推论1在 上的k 阶差商与 的次序无关。
即对任意一个的排列 ，都有 121,,,,+k k y y y y )(x f 121,,,,+k k y y y y 1,,,2,1+k k 121,,,+k i i i ),,,,(),,,,(121=+i i i i k k y y y y f y y y y f证明 由于在 上的k 次插值多项式与在上的k 次插值多项式是相同的，因而首项系数也是相同的。
)(x f 121,,,,+k k y y y y 121,,,,+k k i i i i y y y y 推论2在 上的k 阶差商为 )(x f 121,,,,+k k y y y y ∑+=++-+----=111111121)())(()()(),,,,(k i k i i i i i i i k k y y y y y y y y y f y y y y f (7) 证明 由拉格朗日插值公式知k 次插值多项式为：∑+=++-++---------=1111111111)())(()()())(()()()(k i k i i i i i i k i i i y y y y y y y y y x y x y x y x y f x P 它的首项系数就是（7）式的右边，因而（7）式成立。