当前位置：文档之家 › 第六章 函数的插值方法

第六章 函数的插值方法

  • 格式：doc
  • 大小：802.50 KB
  • 文档页数：29

下载文档原格式

  / 29

合集下载

函数的插值法

函数的插值法

函数的插值法
函数插值法是一种数值分析中的有效算法，它可以用来近似一个函数的值，既可以在给定的若干个点上，也可以在连续的一段区间上。

在插值法中，首先我们需要找到一个合适的插值函数，它可以比较好的拟合给定的点。

有许多常见的插值函数，比如拉格朗日插值函数、牛顿插值函数、指数插值函数以及Hermit插值函数等等。

当确定插值函数之后，就可以使用插值算法，将这些函数应用于给定的点，然后得到函数的极值、极点以及拐点等特征。

例如，如果给定的是三个点，则可以使用牛顿插值算法，将这三个点拟合起来，然后计算该函数的极值、极点和拐点，以此得到函数的一般性表达式。

如果要对函数在连续的一段区间上进行插值，可以使用多项式拟合法来求解，首先在这一段区间内构造一个有限多项式，然后使用这个多项式来拟合给定的点，充分利用多项式的特性，从而得到函数的一般性表达式。

总之，函数插值法是一种有效的数值计算方法，它能够帮助我们求解函数的理论表达式，而不必为此去解方程，大大提高了计算的效率，也使得函数的理论研究变得更加容易。

第6讲(1)插值

第6讲(1)插值

8
其中 Ai 为待定系数，利用li ( xi ) = 1可解得：
Ai = ( xi − x0 )
1 ( xi − xi−1 )( xi − xi+1 )
(xi − xn )
从而
∏ li ( x) =
j≠i
x − xj xi − x j
基本插值函数 (插值基函数)
9
2-2 Lagrange 插值多项式
数类，插值函数 P( x) 满足 P( xi ) = yi (i = 0, , n) ， 即
a0ϕ0 ( xi ) + a1ϕ1 ( xi ) + + anϕn ( xi ) = yi , i = 0, , n
6
若插值基函数{ϕ
i
(
x
)}n i=
0
线性无关，则上述方程组
有唯一的解{ai
}n i=
0
(3)
P(xk ) = I1,
,n ( xk
)
xk x1
− xn+1 − xn+1
+
I2,
,n+1 ( xk
)
xk − x1 xn+1 − x1
=
f
( xk
)
xk x1
− −
xn +1 xn +1
+
f
( xk
)
xk − x1 xn+1 − x1
=
f (xk )
(2 ≤ k ≤ n)
23
注 由上述性质可知， P(x) 是 f (x) 的关于节点 x1, , xn+1的 n 次插值多项式. 它实质上是对两个 n −1次的插值多项式，再经过线性插值求出的.

第六章 函数插值

第六章 函数插值

牛顿前差公式
设 x x0 t h ，则 N n ( x ) N n ( x0 t h) n
Rn ( x)
f (n1) ( ) t(t 1)...( t n)hn1 ,
(n 1)!
k0
( x0, xn )
t k
k f ( x0 )
牛顿后差公式
将节点顺序倒置：
Nn( x) f ( xn ) f [xn , xn1]( x xn ) ... f [xn , ... , x0]( x xn )...( x x1)
Nn(x)
ai = f [ x0, …, xi ]
Rn(x)
牛顿插值公式的优点是：当增加一个节点时，只要再增加 一项就行了，即有递推式: Nk1( x) Nk ( x) ( x x0 )( x x1 )( x xk ) f [ x0 ,, xk , xk1 ]
由插值的唯一性可知 Nn(x) Ln(x)， 故其余项也相同，即
也可以将P2 (x)写成：
P2 (x) a0 a1(x x0 ) a2 (x x0 )( x x1)
令x = x0，由插值条件（6.2），有
P2 ( x0 ) y0 a0
令x = x1，由（6.2），有
f ( x1 ) P2 ( x1) y0 a1( x1 x0 )
a1
y1 y0 x1 x0
这样的分段插值函数在分段上要求多项式次数低而在节点上不仅连续还存在连续的低阶导数我们把满足这样条件的插值函数称为样条插值函数它所对应的曲线称为样条曲线其节点称为样点这种插值方法称为样条插值
第六章 函数插值
6.1 代数插值
设已知某个函数关系y = f (x)在某些离散点上的函数值：
x x0 x1 x2 xn

函数插值法

函数插值法

(xk ,
yk
),(xk+1 ,
yk+1
)
,求一条
由直线的点斜式公式可知：
P
1
( x ) =
y
k
+
y x
k + 1 k + 1
− y − x
k k
( x − x
k
)
把此式按照
yk
和
yk+1
写成两项：
P
1
(x) =
x − xk+1 x − xk yk + yk+1 xk − xk+1 xk+1 − xk
。
yk = f ( xk )
求 一 个 次 数 不 超 过 二 次 的 多 项 式
P2 ( x )
, 使 其 满 足 ,
P2 ( x k − 1 ) = y k − 1
P2 ( x k + 1 ) = y k + 1
。
,
P2 ( x k ) = y k
其几何意义为:已知平面上的三个点
(x k−1
, y
f (x)
在 区 间
[ xk , xk+1 ]
的 端 点 上 的 函 数 值
y k = f ( x k ), y k + 1 = f ( x k + 1 )
，求一个一次函数 。
y = P1 ( x )
使得
y k = P1 ( x k ), y k + 1 = P1 ( x k + 1 )
其几何意义是已知平面上两点 直线过该已知两点。 2.插值函数和插值基函数
要证明插值多项式存在唯一，只要证明参数 其系数行列式不为零即可。 系数行列式为：

插值法数学计算方法

插值法数学计算方法

插值法数学计算方法插值法是一种数学计算方法，用于在已知数据点的基础上，通过构建一条插值曲线来估计未知数据点的值。

插值法可以应用于各种数学问题中，例如逼近函数、插值多项式、差值等。

本文将详细介绍插值法的原理和常见的插值方法。

一、插值法的原理插值法的基本思想是通过已知数据点的函数值来构建一个函数表达式，该函数可以通过插值曲线来估计任意点的函数值。

根据已知数据点的数量和分布，插值法可以采用不同的插值方法来构建插值函数。

插值法的原理可以用以下几个步骤来描述：1.收集已知数据点：首先，需要收集一组已知的数据点。

这些数据点可以是实际测量得到的，也可以是其他方式获得的。

2.选择插值方法：根据问题的特性和数据点的分布，选择适合的插值方法。

常见的插值方法包括拉格朗日插值法、牛顿插值法、埃尔米特插值法等。

3.构建插值函数：通过已知数据点，利用选择的插值方法构建插值函数。

这个函数可以拟合已知数据点，并通过插值曲线来估计未知数据点。

4.估计未知数据点：利用构建的插值函数，可以估计任意点的函数值。

通过插值曲线，可以对未知数据点进行预测，获得相应的数值结果。

二、常见的插值方法1.拉格朗日插值法：拉格朗日插值法基于拉格朗日多项式，通过构建一个具有多项式形式的插值函数来逼近已知数据点。

插值函数可以通过拉格朗日基函数计算得到，式子如下：P(x) = ∑[f(xi) * l(x)], i=0 to n其中，P(x)表示插值函数，f(xi)表示已知数据点的函数值，l(x)表示拉格朗日基函数。

2.牛顿插值法：牛顿插值法基于牛顿差商公式，通过构建一个递归的差商表来逼近已知数据点。

插值函数可以通过牛顿插值多项式计算得到，式子如下：P(x) = f(x0) + ∑[(f[x0, x1, ..., xi] * (x - x0) * (x - x1)* ... * (x - xi-1)] , i=1 to n其中，P(x)表示插值函数，f[x0, x1, ..., xi]表示xi对应的差商。

数值分析第六章插值法PPT学习教案

数值分析第六章插值法PPT学习教案

l0 (x) 9l1(x) 23l2 (x) 3l3 (x)
11 x3 45 x 2 1 x 1
4
4
2
为便于上机计算,常将拉格朗日插值多项式(6.8)改写成
Ln (x)
n k 0
yk
n
i0 ik
x xi xk xi
第16页/共81页
例6.5 已知f(x)的观测数据 x 1234
1 3
3 1
2
f (1.5) p(1.5) 1.25
第13页/共81页
例6.3 已知x=1, 4, 9 的平方根值,
7
用抛物插值公式, (x–x1)(x–x2)
(x–x0)(x–x2)
p2(x) =
y0
+
y1
求
(x0–x1)(x0–x2) (x–x0)(x–x1)
(x1–x0)(x1–x2)
+
第4页/共81页
这惟是一一性个说关明于，待不定论参用数何种a0方, a法1, 来构, a造n 的，n也+不1阶论线用性何方种 程组形,式其来系表数示矩插阵值行多列项式式为,只要满足插值条件(6.1)其结
果都是相互恒等的。
1 x0 x02 x0n
1 V
x1
x12
x1n
n i 1
i 1
(xi x j )
P(x) l0 (x) y0 l1(x) y1 ln (x) yn 事实上，由于每个插值基函数 lk (x)(k 0,1,, n) 都是n次值多项式,所以他们的线性组合
n
P(x) lk (x) yk k 0
（6.8）
是次数不超过n次的多项式 , 称形如（6.8）式的插
值多项式为n次拉格朗日插值多项式。并记为 Ln (x)

数值分析第六章_数值插值方法


M n1 (n 1)!
n1 ( x)
说明：
n=1时，
R1 ( x)

1 2
f
( )2 (x)

1 2
f
( )(x
x0 )(x
x1)
n=2时，
( [x0 , x1])
R2 (x)

1 6
f
( )(x
x0 )(x
x1)(x
x2 )
( [x0 , x2 ])
,
x1,
Hale Waihona Puke xn)1
x1
x12

x1n

n
( xi
ni j1

xj)
1 xn xn2 xnn
因 xi x j (i j) 故上式不为0。
据Cramer法则，方程组解存在且唯一。 故Pn (x)存在且唯一。虽然直接求解上述方程组 可求得插值多项式，但繁琐复杂，一般不用。
得关于a0,a1,…,an的n+1阶线性方程组
a0 a1x0 a0 a1x1
an x0n an x1n

y0 y1
a0 a1xn an xnn yn
其系数行列式是Vandermonde行列式
1 x0 x02 x0n
V
( x0
jk jk
(j,k=0,1)
称l0 (x)及l1 (x)为线性插值基函数。
2. 抛物插值：n=2情形
假定插值节点为x0, x1, x2 ，求二次插值多项式 L2 (x)，使 L2(xj)=yj (j=0,1,2) y= L2 (x)的几何意义就是过 (x0, y0)，(x1, y1) ， (x2, y2)三点的抛物线。 采用基函数方法，设

插值法内插与外插

插值法内插与外插插值法(内插与外插)Interpolation and Extrapolation数学问题探讨当我们在表现资料时，常常会有需要比实际量测点上的值更细密的情况，或者是有需要在范围外预测其值。

比方说天气图的绘制，不论是气压或是雨量，都不可能做到处处都有测量站，又例如我们关心一天之中温度随时间的变化，但是实际上记录气温的动作可能只是每小时一次，则我们要作一个连续的图时，就会用到插值法。

插值法的中心议题是：在我们己具备一组表列数(tabulated value)的情况下，如何得出没被定到之区域的值。

什么样的函数才能被插值，这是数学上讨论的间问题，参见课文p.99第四段。

然而，我们会要用到插值法的场合往往都不知道描述对象背后的函数是什么形式(但相信其有连续的本质)，因此我们也只能尽力求真实。

使用插值法所建立的函数，在表列点上一定要重现原本给定的表列值，否则就不是插值法而是函数近似或曲线拟合的间问题了，它们是不一样的。

插值的作法，很直观地来讲，就是，(1)先从表列值来获得函数f(x)，再(2)用函数f(x)求出我们所要的任何特定x之f(x)函数值。

然而，比较精密且系统化的数值方法却不是用这两个步骤来进行插值，原因是前述两阶段方法对于插值的精密度并没有控制，效率较差，也比较会有进位误差。

一般在做插值法，是从欲插值点x附近的几个表列点xi开始，建立插值函数f(x)，并且也试著网罗更多表列点来插值，看随著项数变多误差会不会变小，如此找出最适合的函数f(x)。

我们会比较希望演算法在从表列值建立插值用函数时，也能提供误差分析以供我们或程式来判断。

毕竟可用的插值函数f(x)并非唯一，而即便是己设定了采用一种方法，如多项式法，也会有该使用多少项才最恰当的问题。

建立插值函数所需之邻近表列值个数，我们称之为插值法的order(阶)，较高阶未必保证得到较合理的插值，这点在多项式插值法尤其如此，要小心注意。

详见课文中之例图上两图实线都是原现象背后的真正值，短虚线代表低阶多项式插值结果，长虚线代表高阶多项式插值结果。

第六章 插值法


的近似值． 7 的近似值．
解 插值条件为 y (4) = 2, y (9) = 3, y (16) = 4
(7 − 9)(7 −16) (7 − 4)(7 −16) (7 − 4)(7 − 9) 7 ≈ y2 (7) = ×2+ ×3 + ×4 (4 − 9)(4 −16) (9 − 4)(9 −16) (16 − 4)(16 − 9) 3 81 2 = + − = 2.6285714 5 35 7
定义6.1 设f(x)在[a,b]上有定义 相异的点 i, i=0,1,2,…,n, 都在 上有定义,相异的点 都在[a,b] 定义 在 上有定义 相异的点x 上，不妨设
a ≤ x0 < x1 < L< xn ≤ b
又设f(x 为 在这些点上的准确值， 又设 i)为f(x)在这些点上的准确值，若存在一个多项式 在这些点上的准确值 若存在一个多项式y(x)，使 ，
第六章 插值法
插值法在数值分析这门课程中是最基础，且应用最广泛的知识． 插值法在数值分析这门课程中是最基础，且应用最广泛的知识． 在工程应用中对于函数 y =f (x)常常不能得到一个具体的解析表达 常常不能得到一个具体的解析表达 它可能是通过实验、 式，它可能是通过实验、测量或者中间计算而得到的一组数据 (xi,f(xi)) i=0,1,2,…,n，或者虽然有函数 =f (x)的解析表达式，但其 的解析表达式， ，或者虽然有函数y 的解析表达式 关系式相当复杂，不便于计算和使用． 关系式相当复杂，不便于计算和使用．因此我们需要用一个比较 简单的函数 y = y(x) 来近似代替数据, ,或近似代替函数 =f (x) ,使 来近似代替数据 或近似代替函数y 使 或近似代替函数
y2(xi ) = f (xi )

函数的插值法

实验二、函数的插值法一、插值算法的基本思想：函数的变化规律往往是通过一组实验数据给出，为了研究此变化规律往往需要求出不在表上的一些函数值。

因此我们希望根据给定的函数表做一个既能反函数)(x f 的特性，又便于计算的简单函数)(x p ，用)(x p 近似)(x f 。

通常选一类较简单的函数作为)(x p ，并使成立。

对n i x f x p i i ,,1,0)()( ==这样确定的)(x p 就是我们希望得到的插值函数，主要有①Lagrange 插值；②Newton 插值；③Hermit 插值；④分段线性插值；⑤三次样条插值。

二、实验要求：1、利用Lagrange 插值公式编写插值多项式程序。

2、给出分段二次和三次多项式的一般表达式，并根据节点选取规则对上述两种插值法进行比较分析。

3、对此插值问题若用牛顿插值多项式结果如何？牛顿插值多项式公式如下：三、目的和意义：1、学会常用的插值计算方法以求函数的近似表达式， 以解决其它科学实验的计算问题；2、通过比较加深理解插值多项式与分段插值多项式的优劣问题；3、熟悉各种插值方法的程序设计；4、通过画出的各种插值曲线比较起插值的光滑度问题。

四、算法步骤：()()nn kkk L x lx y ==∑0()1,()()(0,1,2,,)0,()nj i i j j i j j ix x j il x l x i n j i x x =≠-=⎧===⎨≠-⎩∏1001000100()()[,,]()()[,,]()k nn n j k j j iki n ki ij j j iN x f x f x x x x x f x f x x x xx -==≠==≠=+-=-∑∏∑∏六、Newton 插值1）输入各个节点和节点的函数值ym y y xm x x ,,1,0,,1,0 和； 2）构造均差表，定义)(x f ，[]0,x x f ，[]n x x x f ,,,0 ；3）定义来得)(x f ，[]0,x x f ，[]n x x x f ,,,0 代入Newton 公式，计算N 次插值多项式流程图：七、MATLAB 程序源代码Newton 插值程序源代码function f=newton(x0,y0) %x0为入的节点值，y0相应节点的函数值 n=length(x0) syms x for i=1:nf(i,1)=y0(i) endhx=f(1,1) xx=(1.0) for k=2:n for i=k:nf(i,k)=(f(i,k-1)-f(i-1,k-1))/(x0(i)-x0(i-k+1)) %构造差商表 endxx=xx*(x-x0(k-1)) hx=hx+f(k,k)*xx endf=expand(hx) %多项式展开● 数值例题1）数值插入：x0=[0.40 0.55 0.65 0.80] y0=[0.41075 0.57815 0.69675 0.88811] 运行：newton(x0,y0)2）结果：f=-25542915886587131/11258999068426240000+4572354581681722811/4503599627370496000*x-5029019583899/140737488355328*x^2+444355163233905/2251799813685248*x^3● 数值例题：多项式2x f =在66≤≤-x 取n=6.按等距节点求分段线性插值函数L(x):在matlab 命令窗口中输入如下程序：a=-6;b=6;n=6; f='x^2';fenduan(a,b,n,f)L =[ -10*X-24, -6*X-8, -2*X, 2*X, 6*X-8,10*X-24]数值例题% p60clearx=[.25 .3 .39 .45 .53];y=[.5 .5477 .6245 .6708 .728];% 1%dy=[1 .6868];% x0=[.4 .47];[F,a]=scyt(x,y)当x属于区间[0.25,0.3]时这个函数为f =-6.8*(.30-1.*t)^3-4.9*(t-.25)^3+.26+.95*t当x属于区间[0.3,0.39]时这个函数为f =-2.7*(.39-1.*t)^3-1.9*(t-.30)^3+.30+.85*t当x属于区间[0.39,0.45]时这个函数为f =-2.9*(.45-1.*t)^3-2.2*(t-.39)^3+.33+.77*t当x属于区间[0.45,0.53]时这个函数为f =-1.7*(.53-1.*t)^3-1.4*(t-.45)^3+.35+.71*ta = 0.6325 0.6855当x属于区间[0.25,0.3]时这个函数为f =-6.8*(.30-1.*t)^3-4.9*(t-.25)^3+.26+.95*t当x属于区间[0.3,0.39]时这个函数为f =-2.7*(.39-1.*t)^3-1.9*(t-.30)^3+.30+.85*t当x属于区间[0.39,0.45]时这个函数为f =-2.9*(.45-1.*t)^3-2.2*(t-.39)^3+.33+.77*t当x属于区间[0.45,0.53]时这个函数为f =-1.7*(.53-1.*t)^3-1.4*(t-.45)^3+.35+.71*ta = 0.6325 0.6855clearx=[.25 .3 .39 .45 .53];y=[.5 .5477 .6245 .6708 .728];% 2% x0=[.4 .47];% ddy=[0 0];[F,a]=scyt(x,y)当x属于区间[0.25,0.3]时这个函数为f =-6.3*(t-.25)^3+.26+.97*t当x属于区间[0.3,0.39]时这个函数为f =-3.5*(.39-1.*t)^3-1.6*(t-.30)^3+.30+.84*t 当x属于区间[0.39,0.45]时这个函数为f =-2.4*(.45-1.*t)^3-2.9*(t-.39)^3+.32+.77*t 当x属于区间[0.45,0.53]时这个函数为f =-2.1*(.53-1.*t)^3+.36+.70*ta =0.6324 0.6855% p56clearx=[27.7 28 29 30];y=[4.1 4.3 4.1 3];[F,a]=scyt(x,y)% 1% dy=[3.0 -4];% x0=27.7:0.1:30;当x属于区间[27.7,28]时这个函数为f =-13.*(28.-1.*t)^3+.22*(t-28.)^3+19.-.53*t当x属于区间[28,29]时这个函数为f =.66e-1*(29.-1.*t)^3+.14*(t-28.)^3+12.-.27*t当x属于区间[29,30]时这个函数为f =.14*(30.-1.*t)^3-1.5*(t-29.)^3-12.+.56*ta =[4.1000 4.2956 4.3357 4.3000 4.2550 4.2144 4.1787 4.14824.1234 4.1047 4.0926 4.0875 4.0898 4.1000 4.1167 4.13184.1354 4.1173 4.0678 3.9769 3.8346 3.6310 3.3561 3.0000]clear n=7;x=linspace(-pi,pi,n)y=sin(x);[F,a]=scyt(x,y)% 3% x0=0当x属于区间[-3.1416,-2.0944]时这个函数为f =.15*(t+3.1)^3-3.1-.99*t当x属于区间[-2.0944,-1.0472]时这个函数为f =.15*(-1.0-1.*t)^3+.15*(t+2.1)^3-1.2当x属于区间[-1.0472,0]时这个函数为f =-.15*t^3+.32e-16*(t+1.0)^3+.99*t-.39e-16当x属于区间[0,1.0472]时这个函数为f =.32e-16*(1.0-1.*t)^3-.15*t^3-.39e-16+.99*t当x属于区间[1.0472,2.0944]时这个函数为f =-.15*(2.1-1.*t)^3-.15*(t-1.0)^3+1.2当x属于区间[2.0944,3.1416]时这个函数为f =-.15*(3.1-1.*t)^3+3.1-.99*t a=1.8489e-032或采用：function [S,y]=bkj(x0,y0,p,x,s1,s2)n=length(x0);w=length(x);for i=1:n-1h(i)=x0(i+1)-x0(i);g(i)=y0(i+1)-y0(i);endfor j=2:n-1d(j)=6*((g(j)/h(j))-(g(j-1)/h(j-1)))/(h(j-1)+h(j));u(j)=h(j-1)/(h(j-1)+h(j));r(j)=h(j)/(h(j-1)+h(j)); endif p==1r(1)=1;u(n)=1;d(1)=6/h(1)*(g(1)/h(1)-s1);d(n)=6/h(n-1)*(s2-g(n-1)/h(n-1));elseif p==2r(1)=0;u(n)=0;d(1)=2*s1;d(n)=2*s2;elser(n)=h(1)/(h(n-1)+h(1));u(n)=1-r(n);d(n)=6*((g(1)/h(1))-(g(n-1)/h(n-1)))/(h(n-1)+h(1)); endif p==3for i=1:n-2for j=1:n-2if i==jA(i,j)=2;A(i,j+1)=r(i);A(i+1,j)=u(i+1);endendendA(n-1,n-1)=2;A(1,n-1)=u(1);A(n-1,1)=r(n);for i=1:n-1b(i)=d(i+1);endL=A\b';for i=1:n-1M(i+1)=L(i);endM(1)=M(n);elsefor i=1:n-1for j=1:n-1if i==jA(i,j)=2;A(i,j+1)=r(i);A(i+1,j)=u(i+1);endendendA(n,n)=2;M=A\d';endfor a=1:wif ((x(a)>=x0(1))&(x(a)<=x0(n)))for j=0:n-1if x(a)<x0(j+1)s=sym('x');R=poly2str(sym2poly(M(j)*(x0(j+1)-s)^3/(6*h(j))+M(j+1)*(s-x0(j))^3/(6 *h(j))+(y0(j)-(M(j)*h(j)^2)/6)*(x0(j+1)-s)/h(j)+(y0(j+1)-(M(j+1)*h(j) ^2)/6)*(s-x0(j))/h(j)),'x');for t=1:length(R)S(a,t)=R(t);endy(a)=M(j)*(x0(j+1)-x(a))^3/(6*h(j))+M(j+1)*(x(a)-x0(j))^3/(6*h(j))+(y 0(j)-(M(j)*h(j)^2)/6)*(x0(j+1)-x(a))/h(j)+(y0(j+1)-(M(j+1)*h(j)^2)/6) *(x(a)-x0(j))/h(j);break;endendelsey(a)=0;disp( x(a)) ;disp('this number is not exist it this erea') end end.3301.4293.4281.47.27)(i i x f x.4)30(0.3)7.27(''-==S S (注释: 第一种边界条件)X=[27.75 28.85 29.95]Y= [4.2222 4.0877 3.1888]得到的S(x)= 13.293 x^3 - 1116.4136 x^2 + 31253.5671 x - 291635.1316)28,7.27(∈x ： 0.072277 x^3 - 5.8733 x^2 + 158.4366 x - 1413.9139 )29,28(∈x ： -1.6574 x^3 + 144.6109 x^2 - 4205.604 x + 40771.8119CPU 时间: 0.0250s结果分析：本程序主要特点是通用性比较长,三种边界条件都可运用,而且可一次性求出多个点的值.比拉格郎日插值,牛顿插值和分段低次插值它有光滑性,而且很近似原函数值;。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

6.1 插值问题及其误差6.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)例 6.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)6.2 拉格朗日（Lagrange ）插值及其MATLAB 程序6.2.1 线性插值及其MATLAB 程序例 6.2.1 已知函数)(x f 在]3,1[上具有二阶连续导数，5)("≤x f ，且满足条件2)3(,1)1(==f f .求线性插值多项式和函数值)5.1(f ，并估计其误差.解 输入程序>> 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 为-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 =1.8750例6.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.6.2.2 抛物线插值及其MATLAB 程序例 6.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) ,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=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) ,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),P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6; Y = polyval(P,x)运行后输出基函数l 01、l 11和l 21及其插值多项式的系数向量P 、插值多项式L 和插值Y 为l0 =228155022448185/281474976710656*x^2-2150310427208497/1125899906842624*x+1l1 =-228155022448185/140737488355328*x^2+5734161139222659/2251799813685248*xl2 =228155022448185/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.6.2.3 n 次拉格朗日（Lagrange ）插值及其MATLAB 程序例 6.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+1l2 =-1/3*x^3+4/3*x ，l3 =1/8*x^3+1/8*x^2-1/4*x输入程序>> 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(-∈ξ.6.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例6.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)运行后输出五次拉格朗日插值多项式L 及其系数向量C ，基函数l 及其系数矩阵L 1如下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(∈ξ.6.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:np=1.0; q1=1.0; c1=1.0;for j=1:nif j~=kp=p*(z-X(j))/(X(k)-X(j));endq1=abs(q1*(z-X(j)));c1=c1*j;ends=p*Y(k)+s;endy(i)=s;endR=M*q1/c1;例 6.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.6.3 牛顿（Newton ）插值及其MATLAB 程序6.3.3 牛顿插值多项式、差商和误差公式的MATLAB 程序求牛顿插值多项式和差商的MATLAB 主程序function [A,C,L,wcgs,Cw]= newploy(X,Y)n=length(X); A=zeros(n,n); A(:,1)=Y';s=0.0; p=1.0; q=1.0; c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); endb=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1;endC=A(n,n); b=poly(X(n)); q1=conv(q1,b);for k=(n-1):-1:1C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k);endL(k,:)=poly2sym(C); Q=poly2sym(q1);syms Mwcgs=M*Q/c1; Cw=q1/c1;例6.3.3 给出节点数据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 ,作五阶牛顿插值多项式和差商，并写出其估计误差的公式.解 （1）保存名为newpoly.m 的M 文件.（2）输入MATLAB 程序Y=[17.03 7.24 1.05 2.03 17.06 23.05];[A,C,L,wcgs,Cw]= newdcw (X,Y)运行后输出差商矩阵A ,五阶牛顿插值多项式L 及其系数向量C , 插值余项公式L 及其向量C w 如下A =17.0300 0 0 0 0 07.2400 -8.5130 0 0 0 01.0500 -6.1287 1.1039 0 0 02.0300 0.97033.5144 0.7604 0 017.0600 14.8812 6.8866 1.1129 0.0843 023.0500 4.9098 -4.4715 -3.5056 -1.0867 -0.2169C =-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =-7813219284746629/36028797018963968*x^5+583849564517807/9007199254740992*x^4+593245028711263/281474976710656*x^3+3823593773002357/1125899906842624*x^2-321902673270315/70368744177664*x+308328649211299/281474976710656wcgs =1/720*M*(x^6-79/25*x^5-14201/2500*x^4+4934097026900981/281474976710656*x^3+154500712237335/35184372088832*x^2-8170642380559269/562949953421312*x+5212760744134241/36028797018963968)Cw =0.0014 -0.0044 -0.0079 0.0243 0.0061 -0.0202 0.0002即L =1.0954-4.5745*x +3.3960*x ^2+2.1076*x ^3+0.0648*x ^4-0.2169*x ^5.估计其误差的公式为)(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(∈ξ.例6.3.4 求函数7)(-=x f e 5/x -在]6,2[上五阶牛顿插值多项式，估计其误差的公式和误差限公式.用它们计算)156.3(f ，并估计其误差.解 （1）输入MATLAB 程序>> X=2:4/5:6; Y=-7*exp(-X/5);[A,C,L,wcgs,Cw]= newploy(X,Y), x1=2:0.001:6;M=max(-7*exp(-x1/5)/(5^6)),运行后输出差商矩阵A , 五阶牛顿插值多项式L 及其系数向量C , 插值余项公式L 及其向量C w 如下A =-4.6922 0 0 0 0 0-3.9985 0.8672 0 0 0 0-3.4073 0.7390 -0.0801 0 0 0-2.9035 0.6297 -0.0683 0.0049 0 0-2.4742 0.5366 -0.0582 0.0042 -0.0002 0-2.1084 0.4573 -0.0496 0.0036 -0.0002 0.0000C =0.0000 -0.0004 0.0089 -0.1389 1.3985 -6.9991L =9721799720875/1152921504606846976*x^5-3503994098647815/9223372036854775808*x^4+160742008798419/18014398509481984*x^3-1251152213853501/9007199254740992*x^2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312wcgs =1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+7276634802928539/2199023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/2199023255552)Cw =0.0014 -0.0333 0.3256 -1.6533 4.5959 -6.6159 3.8438M =-1.3494e-004（2）输入MATLAB 程序>> syms xwcgs1=1/720*M*1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+7276634802928539/2199023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/2199023255552),运行后输出误差限公式wcgs1如下wcgs1 =5565367633581443/158456325028528675187087900672*x^6-16696102900744329/19807040628566084398385987584*x^5+1630652716639362799/198070406285660843983859875840*x^4-517579189923074199/12379400392853802748991242240*x^3+40497147813610772932297365501777/348449143727040986586495598010130648530944*x^2-29148396970323855270001164718917/174224571863520493293247799005065324265472*x+33870489914310283926665276375603/348449143727040986586495598010130648530944（3）输入MATLAB 程序>> x=3.156;y=9721799720875/1152921504606846976*x^5-3503994098647815/9223372036854775808*x^4+160742008798419/18014398509481984*x^3-1251152213853501/9007199254740992*x^2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312wcgs2=1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+7276634802928539/2199023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/2199023255552)运行后输出)156.3(f 的近似值y ，及其误差限wcgs 2如下y =-3.7237wcgs2 =-2.4764e-0076.3.4 牛顿插值及其误差估计的MATLAB 程序牛顿插值及其误差估计的MATLAB 主程序function [y,R]= newcz(X,Y,x,M)n=length(X); m=length(x);for t=1:mz=x(t); A=zeros(n,n);A(:,1)=Y';s=0.0; p=1.0; q1=1.0; c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C); C(d)=C(d)+A(k,k);endy(k)= polyval(C, z);endR=M*q1/c1;例6.3.5 已知5.030sin =，1707.045sin = ，0866.060sin = ，用牛顿插值法求40sin 的近似值，估计其误差，并与例 6.2.6的计算结果比较.解 方法1（牛顿插值及其误差估计的MATLAB 主程序）输入MATLAB 程序>> x=2*pi/9;M=1; X=[pi/6 ,pi/4, pi/3];Y=[0.5,0.7071,0.8660]; [y,R]= newcz(X,Y,x,M)运行后输出y = R =0.6434 8.8610e-004方法2 （求牛顿插值多项式和差商的MATLAB 主程序）输入MATLAB 程序>> x=2*pi/9; X=[pi/6 ,pi/4, pi/3];Y=[0.5,0.7071,0.8660]; M=1;[A,C,L,wcgs,Cw]= newploy(X,Y), y=polyval(C,x)运行后输出结果A =0.5000 0 00.7071 0.7911 00.8660 0.6070 -0.3516C =-0.3516 1.2513 -0.0588L =-1583578379808357/4503599627370496*x^2+1408883391907005/1125899906842624*x-132405829044691/2251799813685248wcgs =1/6*M*(x^3-3/4*x^2*pi+4012734077357799/2251799813685248*x-7757769783530263/18014398509481984)Cw =0.1667 -0.3927 0.2970 -0.0718y =0.6434上述两种方法计算y 的结果相同.6.3.5 牛顿插值法的MATLAB 综合程序求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序function [y,R,A,C,L]=newdscg(X,Y,x,M)n=length(X); m=length(x);for t=1:mz=x(t); A=zeros(n,n);A(:,1)=Y';s=0.0; p=1.0; q1=1.0; c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endy(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);例6.3.6 给出节点数据00.27)00.4(=-f ，00.1)00.0(=f ，00.2)00.1(=f ，00.17)00.2(=f ，作三阶牛顿插值多项式，计算)345.2(-f ，并估计其误差.解 首先将名为newdscg.m 的程序保存为M 文件，然后在MA TLAB 工作窗口输入程序>> syms M,X=[-4,0,1,2]; Y =[27,1,2,17]; x=-2.345;[y,R,A,C,P]=newdscg(X,Y,x,M)运行后输出插值y )345.2(-≈f 及其误差限公式R ，三阶牛顿插值多项式P 及其系数向量C ，差商的矩阵A 如下y =22.3211R =1323077530165133/562949953421312*M （即R =2.3503*M ）A=27.0000 0 0 01.0000 -6.5000 0 02.0000 1.0000 1.5000 017.0000 15.0000 7.0000 0.9167C =0.9167 4.2500 -4.1667 1.0000P =11/12*x^3+17/4*x^2-25/6*x+1例 6.3.7 将区间 [0,π/2] 分成n 等份)3,2(=n ，用x x f y sin )(==产生1+n 个节点，求二阶和三阶牛顿插值多项式)(2x P 和)(3x P .用它们分别计算sin (π/7) (取四位有效数字)，并估计其误差.解 首先将名为newdscg.m 的程序保存为M 文件，然后在MA TLAB 工作窗口输入程序>> X1=0:pi/4:pi/2; Y1 =sin(X1); M=1; x=pi/7;X2=0:pi/6:pi/2; Y2 =sin(X2);[y1,R1,A1,C1,P2]=newdscg(X1,Y1,x,M),[y2,R2,A2,C2,P3]=newdscg(X2,Y2,x,M)运行后输出插值y 1=)7/sin()7/(2π≈πP 和y 2=)7/sin()7/(3π≈πP 及其误差限R 1和R 2，二阶和三阶牛顿插值多项式P 2和P 3及其系数向量C 1和C 2，差商的矩阵A 1和A 2如下y1 =0.4548R1 =0.0282A1 =0 0 00.7071 0.9003 01.0000 0.3729 -0.3357C1 =-0.3357 1.1640 0P2 =-3024156947890437/9007199254740992*x^2+163820246322191/140737488355328*xy2 =0.4345R2 =9.3913e-004A2 =0 0 0 00.5000 0.9549 0 00.8660 0.6991 -0.2443 01.0000 0.2559 -0.4232 -0.1139C2 =-0.1139 -0.0655 1.0204 0P3 =-1025666884451963/9007199254740992*x^3-4717668559161127/72057594037927936*x^2+4595602396951547/4503599627370496*x6.4 埃尔米特（Hermite ）插值及其MATLAB 程序6.4.3 埃尔米特插值多项式和误差公式的MATLAB 程序求埃尔米特插值多项式和误差公式的MATLAB 主程序function [Hc, Hk,wcgs,Cw]= hermite (X,Y,Y1)m=length(X); n=M1;s=0; H=0;q=1;c1=1; L=ones(m,m);G=ones(1,m);for k=1:n+1V=1;for i=1:n+1if k~=is=s+(1/(X(k)-X(i)));V=conv(V,poly(X(i)))/(X(k)-X(i));endh=poly(X(k)); g=(1-2*h*s); G=g*Y(k)+h*Y1(k);endH=H+conv(G,conv(V,V)); b=poly(X(k));b2=conv(b,b);q=conv(q,b2); t=2*n+2;Hc=H;Hk=poly2sym (H); Q=poly2sym(q);endfor i=1:tc1=c1*i;endsyms M,wcgs=M*Q/c1; Cw=q/c1;例6.4.3 给定函数)(x f 在点2/,4/,6/210π=π=π=x x x 处的函数值5.0)(0=x f , 1707.0)(1=x f ，1)(2=x f 和导数值0866.0)(0'=x f ，1707.0)(1'=x f ，0000.0)(0'=x f ，且1)()6(≤x f ，求函数)(x f 在点210,,x x x 处的五阶埃尔米特插值多项式)(5x H 和误差公式，计算)1.567(f 并估计其误差. 解 （1）保存名为hermite.m 的M 文件.（2）在MATLAB 工作窗口输入程序>>X=[pi/6,pi/4,pi/2]; Y=[0.5,0.7071,1];Y1=[0.8660,0.7071,0]; [Hc, Hk,wcgs,Cw]= hermite (X,Y,Y1) 运行后输出五阶埃尔米特插值多项式H k 及其系数向量H c ，误差公式wcgs 及其系数向量C w 如下Hc =1.0e+003 *0.1912 -0.9273 1.6903 -1.4380 0.5751 -0.0866 Hk =6725828781679091/35184372088832*x^5-4078286086775209/4398046511104*x^4+7434035571017927/4398046511104*x^3-3162205449085973/2199023255552*x^2+5058863928652835/8796093022208*x-6094057839958843/70368744177664wcgs =1/720*M*(x^6-11/6*x^5*pi+7446708432019761/562949953421312*x^4-4363745503235773/281474976710656*x^3+21569239021155/2199023255552*x^2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)Cw =0.0014 -0.0080 0.0184 -0.0215 0.0136 -0.0044 0.0006 当M x f =≤1)()6(时的误差公式为R=0.001 4* x ^6-0.008 0*x ^5+0.018 4*x ^4-0.021 5*x ^3+0.013 6*x ^2-0.004 4*x +0.000 6（3）在MATLAB 工作窗口输入程序>> x=1.567;M=1;Hk=6725828781679091/35184372088832*x^5-4078286086775209/4398046511104*x^4+7434035571017927/4398046511104*x^3-3162205449085973/2199023255552*x^2+5058863928652835/8796093022208*x-6094057wcgs=1/720*M*(x^6-11/6*x^5*pi+7446708432019761/562949953421312*x^4-4363745503235773/281474976710656*x^3+21569239021155/2199023255552*x^2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)运行后输出)1.567(f 的近似值Hk 及其误差wcgs 如下 Hk =2.5265wcgs =1.3313e-0086.5 分段插值及其MATLAB 程序6.5.1 高次插值的振荡和MATLAB 程序例6.5.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)(x x 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;endy(i)=s;end（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;y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));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;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)'),pausen=3; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);y1=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 的图形(略).6.5.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;yi = interp1(x0,y0,xi);'o'':'legend('节点(xi,yi)','分段线性插值函数','被插值函数y=cosx ')title('y=cosx 及其分段线性插值函数和节点的图形')例6.5.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.048642533936656.5.5 用MATLAB 计算有关分段线性插值的误差例6.5.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;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((sqrt(3)+sin(2.*x))./(3+4*x.^2))); si=xxczhjt1(x0,y0,xi,x,y);Ri=abs((fi-si)./fi);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 工作窗口输入程序y=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.001698934907116.6 分段埃尔米特插值及其MATLAB 程序6.6.2 分段埃尔米特插值的MATLAB 程序调用格式一：YI=interp1(X,Y ,XI,'pchip')调用格式二：YI=interp1(X,XI,'pchip')例6.6.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 (略).6.6.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;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 及其分段埃尔米特插值函数和节点的图形')>> 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 及其分段埃尔米特插值函数和节点的图形')例6.6.6 设函数211)(x x 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 的图形(略).例6.6.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 的图形(略).6.6.4 用MATLAB 计算有关分段埃尔米特插值的误差例6.6.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 工作窗口输入程序>> syms h,x=-1:0.0001:1;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 工作窗口输入程序运行后输出误差限为R =0.06250000000000例6.6.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.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2)+16.*(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)).^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)).(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)+si n(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))./(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.*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。