第七章 matlab在数值分析中的应用

主要介绍有关插值法的一些基本概念, 及多项式插值的基础理论
一、插值问题
对函数f (x), 其函数形式可能很复杂, 且不利于在计算机上
运算, 假如可以通过实验或测量 , 可以获得f ( x )在区间[ a , b ] 上的一组n + 1个不同的点
a ≤ x0 < x1 < x2 < < xn ≤ b
n
n = 2 ,4 ,6 ,8 ,10
%lagrangen.m function y=lagrangen(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i);s=0; for k=1:n L=1; for j=1:n if j~=k L=L*(z-x0(j))/(x0(k)-x0(j)); end end s=s+L*y0(k); end y(i)=s; end y;
数值分析方法的主要任务： 数值分析方法的主要任务： 1.将计算机上不能执行的运算化为在计算机 将计算机上不能执行的运算化为在计算机 上可执行的运算 2.针对所求解的数值问题研究在计算机上可 针对所求解的数值问题研究在计算机上可 执行的且有效的计算公式 3.因为可能采用了近似等价运算，故要进行 因为可能采用了近似等价运算， 因为可能采用了近似等价运算 误差分析， 误差分析，即数值问题的性态及数值方法的 稳定性
C= conv(a,b)。 例： 求a=x3+3x2+2x1+1和b=4x3+3x2+9x1+10的乘积。 a=[1 3 2 1]; b=[4 3 9 10]; c=conv(a,b)
2. 多项式的除法dconv 多项式的除法的运算实质就是多项式系数的解卷积 运算。向量a对向量c进行解卷积得到的商向量q和余 向量r，满足： k c(k)-r(k)= Σ a(i)q(k+1-i) i=1 在Matlab，完成此功能的命令是dconv： [q,r]= dconv(c,a) 其中： c=conv(a,q)+r
整体误差的大小反映了插值函数的好坏 为了使插值函数更方便在计算机上运算，一般插值函 数都使用代数多项式和有理函数 讨论的就是代数插值多项式 二、代数插值多项式的存在唯一性
设函数 y = f ( x ) 在区间[ a , b ]上的代数插值多项式为
Pn ( x ) = a0 + a1 x + a2 x 2 + + an x n
其中
( x xi ) l j (x ) = ∏ i = 0 ( x j xi )
n i≠ j
称( Ln ( x) 为 y = f ( x) 的Lagrange插值多项式
称 l j ( x) (i = 0,1, , n) 为n次Lagrange 插值基函数
1 设函数 f ( x ) = , x ∈ [ 5 ,5] 例. 2 1+ x 10 将[ 5 ,5]n等份取n + 1个节点xi = 5 + ih , h = , i = 0 ,1, , n n
例： 求a=x3+3x2+2x1+1和b=4x3+3x2+9x1+10的乘和商 的导数。 a=[1 3 2 1]; b=[4 3 9 10]; polyder(a,b) [q,d]=polyder(b,a)
插值法
要点 用简单的函数(如多项式函数)作为一个 复杂函数的近似，最简单实用的方法就是 插值
上的函数值 yi = f ( xi ),
i = 0 ,1,2 , , n
能否存在一个性能优良、便于计算的函数
比如多项式函数 P( x )
P( xi ) = yi
i = 0 ,1 , 2 , , n
------(1)
并且用P( x )近似代替f ( x )
这就是插值问题, (1)式为插值条件,
称函数P( x )为函数f ( x )的插值函数
4. 多项式的求值polyval
在Matlab中，多项式求值有两种方式：按数组规 则和按矩阵规则计算多项式值。 y=polyval(p,x)：按数组规则运算。用来计算多 项式p在自变量x处的值，可以是数、向量和矩 阵。如果x是向量或矩阵，则该命令将对x的每 个元素计算p的值，结果大小和x的大小相同
0 .2
0 .1 0 0 0 .5 1 1 .. 5 1 5 2 2 2 .. 5 2 5 3 3 3 .. 5 3 5
xx x
对于被插函数f ( x )和插值函数P( x ) 在节点xi处的函数值必然相等
但在节点外P( x )的值可能就会偏离f ( x ) 因此P( x )近似代替f ( x )必然存在着误差
虽然线性方程组(4)推出的插值多项式存在且唯一 但通过解线性方程组(4)求插值多项式却不是好方法
y = f ( x)在节点xi (i = 0,1, , n)上, 以l j ( x) (i = 0,1, , n) 为插值基函数的插值多项式(记为Ln ( x))为
Ln ( x ) = y0l0 ( x ) + y1l1 ( x ) + + ynln ( x )
例： p=[2 3 3 0 6]; x=[1 -1]; y=polyval(p,x)
5. 多项式的部分分式展开
没有多重根： B(s) R(1) R(2) R(n) ---- = -------- + -------- + ... + -------- + K(s) A(s) s - P(1) s - P(2) s - P(n) 有多重根时，则包含有：
试就n = 2 ,4 ,6 ,8 ,10作f ( x )的n次Lagrange插值多项式
并作图比较。 解:
1 yi = f ( xi ) = 1 + xi2
作n次Lagrange插值多项式
n ( x xi ) 1 Ln ( x ) = ∑ 1 + x 2 ∏ ( x j xi ) j =0 i =0 j i≠ j
如果P(x )为多项式函数, 则称之为插值多项式
称点 xi , i = 0 ,1,2 , , n为插值节点
称区间[ a , b ]为插值区间 如函数y = sin x , 若给定[0 ,π ]上5个等分点
其插值函数的图象如图
s in x 的
1
的 的
y
y
0 .9
0 .8 0 .7 0 .6 0 .5 0 .4 0 .3
例：
A=[6 -8 6;1 0 0;0 2 0]; p=poly(A), roots(p), D=eig(A)
多项式求根 r=roots(p)：求多项式p的根，r是根组成 的向量。 例： p=[2 3 3 0 6]; roots(p)
多项式的基本运算
1. 多项式的乘法conv 多项式的乘法本质上是多项式系数向量的卷积运算。 卷积的定义为： k c(k)= Σ a(i)b(k+1-i) i=1 在Matlab，完成此功能的命令是conv：
R(j) R(j+1) R(j+m-1) -------- + ------------ + ... + -----------s - P(j) (s - P(j))^2 (s - P(j))^m
在Matlab中完成此任务的指令是residue： [r,p,k]=residue(b,a)：求b(s)/a(s)的部分分式展开，返 回值r是留数，p是极点，k是整式向量，a,b分别是 分母和分子； [b,a]=residue(r,p,k)：将展开式还原为两个多项式b(s) 和a(s)相除的形式。
n 1
n
1 xn
2 n xn xn
由Cramer法则，线性方程组(4)有唯一解 定理1. 若插值节点 xi ≠ x j (i ≠ j ),
Pn ( xi ) = yi i = 0 ,1,2 , , n
则满足插值条件 --------(3) --------(2)
的插值多项式 Pn ( x ) = a0 + a1 x + a2 x 2 + + an x n 存在且唯一。
2 n a0 + a1 x0 + a2 x0 + + an x0 = y0
--------(4)
上述方程组的系数行列式为n+1阶Vandermond行列式
1 x0
2 n x0 x0 2 1 n 1
xi ≠ x j 1 x1 x x = V= ∏ j∏1( x j xi ) ≠ 0 i =0 =i +
数值分析的应用： 解线性方程组 插值与拟合 数值积分与数值微分 非线性方程的根 常微分方程的数值解 ……
多项式运算 （62~68页）
多项式的表示和创建
在Matlab中，任意的多项式都是用一个行 向量表示，将多项式的系数按降幂排列存 放在行向量中。多项式p(x)=a0xn+a1xn-1+… +an-1x+an对应的行向量为： P=[a0 a1 … an]。
例: a=[1 3 2 1]; c=[4 15 26 47 51 29 10]; [q,r]=deconv(c,a)
3. 多项式的加减 对多项式所对应的向量利用数组加减运算 规则执行加减运算。
注：多项式所对应的向量如果大小不相等， 必须用0补齐。 例： 求a=x3+3x2+2x1+1和b=4x3+3x2+9x1+10的和与差。 a=[1 3 2 1]; b=[4 3 9 10]; c=a+b, c1=a-b,
2
不同次数的Lagrange插值多项式的比较图
f(x )= 1/(1+ x 2 ) n= 10
1.5
1