第二章 函数的插值
- 格式:pps
- 大小:327.50 KB
- 文档页数:15
数值分析--第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)。
第⼆章函数的插值⽅法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,xj,x,y,dy1,dyn)S = spline(x0,[dy1,y0,dyn],xj);Sn = spline(x0,[dy1,y0,dyn],x);plot(x0,y0,'o',x,Sn,'-',xj,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')。
第二章函数的插值学习目标:掌握多项式插值的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)式成立。