常州大学数值分析第三章

第三章作业

1.设节点x 0=0,x 1=π/8,x2=π/4,x3=3π/8,x4=π/2,试适当选取上述节点，用拉格朗日插值法分别构造cosx 在区间[0，π/2]上的一次、二次、四次差值多项式P 1（x ）,P 2(x)和P 4（x)，并分别计算P 1（π/3）,P 2(π/3)和P 4（π/3). 解： x0 x1 x2 x3 x4 x

π/8 π/4 3π/8 π/2 y=cosx 1

0.923879

0.707106

0.382683

（1）选择x0=0,x4=π/2的节点

y0=cosx0=1,y4=cosx4=0,可得

)

()

()()()(0101

1010

1x x x x y x x x x y x P --+--=，即 333333

.0)3/(1636620.0)(11≈+-≈πP x x P

（2）选择x0=0,x2=π/4,x4=π/2的节点

y0=cosx0=1,y2=cosx2=0.707106,y4=cosx4=0,可得

)

)(())(())(()

)(())(())(()(1202102

2101201

2010210

1x x x x x x x x y x x x x x x x x y x x x x x x x x y x P ----+----+----=,即

145968

.1)3/(1511124.5482067.1)(222≈++-≈πP x x x P

（3）选择x0=0,,x1=π/8,x2=π/4,x3=3π/8，x4=π/2的节点y0=cosx0=1,y1=cosx1=0.923879,y2=cosx2

=0.707106,y3=cosx3=0.382683,y4=cosx4=0可得

)

(

)(4

,04

4∏

∑≠==--=i

j j j

i j i i x x x x y x P ， 得

P3(x)=1+0.0031x-0.51542x +0.02423

x +0.02

844

x

4(3) 0.5001P π=/

7.解：

选取0123=0=1=2=3x x x x ，，，为节点 >> T0=[0.0 0.5];x=[1 2 3]';y=[1.25 2.75 3.5]';x0=2.8;T=aitken(x,y,x0,T0) T =

0.0000 0.5000 0.0000 0.0000 0.0000 1.0000 1.2500 2.6000 0.0000 0.0000 2.0000 2.7500 3.6500 4.4900 0.0000 3.0000 3.5000 3.3000 3.2300 3.4820 16

1）拉格朗日差值

．选取

函数

],[),sin()cos(ππ-∈+=x x x y

x0=-pi:0.5*pi:pi; y0=cos(x0);

x=-pi:0.05*pi:pi;

if length(x0)~=length(y0)

error('The length of x0 must be equal to it of y0'); end

w=length(x0); n=w-1;

L=zeros(w,w); for k=1:n+1 V=1;

for j=1:n+1 if k~=j if

abs(x0(k)-x0(j))

error('Divided by Zero,there are two nodes are the same'); end

V=conv(V,poly(x0(j)))/(x0(k)-x0(j));

end end L(k,:)=V;

end

C=y0*L;

Y=polyval(C,x); y=cos(x)+sin(x); r=y-Y;

plot(x,Y,'r--',x,y,'b-',x,r,'k-.','LineWidth',2);

legend('Lagrange polynomial','The

original f(x) ','Error',0)

2）牛顿差值

．选取函数]5,5[,11

2

-∈+=

x x y >> x0=-5:1:5;

>>y0=1./(1+x0.*x0); >>x=-5:0.1:5;

>>if length(x0)~=length(y0)

>>error('The length of x0 must be equal to it of y0'); end

>>n=length(x0); >>D=zeros(n,n); >>D(:,1)=y0'; >>for j=2:n for k=j:n

if abs(x0(k)-x0(k-j+1))

error('Divided by Zero,there are two nodes are the same'); end

>>D(k,j)=(D(k,j-1)-D(k-1,j-1))/(x0(k )-x0(k-j+1)); end end

C=D(n,n);

for k=(n-1):-1:1

C=conv(C,poly(x0(k))); m=length(C);

C(m)=C(m)+D(k,k); end

Y=polyval(C,x); y=1./(1+x.*x+20*x); r=y-Y;

>plot(x,Y ,'r--',x,y,'b-',x,r,'k-.','LineWidth',2); >legend('Newton polynomial','The original f(x)

','Error',0)