数值分析实验和程序

(2)、试构造Lagrange多项式L6(x)。 结果为f(1.8)?0.164762, f(6.15)?0.001266 在主显示区，输入程序： >> X=[1 2 3 4 5 6 7];

>> Y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001]; >> [A1,LN,L1,B1]=lagrange1(X,Y) >> plot(X,A1); >> F=poly2sym(A1)

运行后，输出结果的Lagrange多项式L6(x)的结果： A1 =

0.0001 -0.0016 0.0186 -0.1175 0.4419 -0.9683 0.9950 F =

(4304240283865561*x^6)/73786976294838206464 (7417128346304051*x^5)/4611686018427387904 199/200

+

(223*x^4)/12000

- -

(2821*x^3)/24000 + (994976512675275*x^2)/2251799813685248 - (19367*x)/20000 +

Lagrange多项式L6(x)的图如下：

二、计算函数值

计算函数值的主程序：lagrangezhi.m function [y,R]=lagrangezhi(X,Y,x,M) n=length(X); m=length(x); for i=1:m z=x(i);s=0.0; for k=1:n

p=1.0; q1=1.0; c1=1.0; for j=1:n if j~=k p=p*(z-X(j))/(X(k)-X(j)); end

q1=abs(q1*(z-X(j)));c1=c1*j; end s=p*Y(k)+s; end y(i)=s; end R=M*q1/c1;

（1）、计算f(0.596)、f(0.99) 的值。 在主显示区，输入程序：

>> x=0.596; M=1; X=[0.4,0.55,0.65,0.80,0.95,1.05]; >> Y=[0.41075,0.57815,0.69675,0.90,1.00,1.25382]; >> [y,R]=lagrangezhi(X,Y,x,M) 运行结果： y = 0.6257 R =

2.2170e-008

在主显示区，输入程序：

>> x=0.99; M=1; X=[0.4,0.55,0.65,0.80,0.95,1.05]; >> Y=[0.41075,0.57815,0.69675,0.90,1.00,1.25382]; >> [y,R]=lagrangezhi(X,Y,x,M) 运行结果： y = 1.0542 R =

5.5901e-008

（2）、计算f(1.8)、f(6.15)的值 在主显示区，输入程序： >> x=1.8; M=1; X=[1,2,3,4,5,6,7];

>> Y=[0.368,0.135,0.050,0.018,0.007,0.002,0.001]; >> [y,R]=lagrangezhi(X,Y,x,M) 运行结果： y = 0.1648 R = 0.0059

在主显示区，输入程序： >> x=6.15; M=1; X=[1,2,3,4,5,6,7];

>> Y=[0.368,0.135,0.050,0.018,0.007,0.002,0.001];

>> [y,R]=lagrangezhi(X,Y,x,M) 运行结果： y = 0.0013 R = 0.0042

三、Newton插值多项式

Newton插值多项式主程序M文件：Newton.m function [A,C,L,wcgs,Cw]= Newton(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:n for i=j:n

A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end

b=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1; end

C=A(n,n); b=poly(X(n)); q1=conv(q1,b); for k=(n-1):-1:1

C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k); end

L(k,:)=poly2sym(C); Q=poly2sym(q1); syms M

wcgs=M*Q/c1; Cw=q1/c1; 在主显示区，输入L5(x)的程序: >> x=[0.4 0.55 0.65 0.80 0.95 1.05];

>> y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382]; >> [A,C,L,wcgs,Cw]= Newton(x,y) >> syms x; >> ezplot(L,[0 1.1]); 运行结果如下，得到L5(x) A =