数值分析实验报告2

实验名称 插值法

实验目的

（1）学习并熟练掌握MA TLAB 语言的编程；

（2）通过课程实习能够应用MATLAB 软件来计算函数的插值，了解函数插值方法。 实验原理

牛顿差商形式多项式

P(x)=f(x0)+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+…+f[x0,x1,x2…xn](x-x0)…(x-xn-1) 牛顿插值多项式的余项 Rn(x)=f[x0,x1,x2…xn]wn+1(x) 实验题目

{1}已知函数在下列各点的值为

i x 0.2 0.4 0.6 0.8 1.0 ()i f x 0.98 0.92 0.81 0.64 0.38

试用4次牛顿插值多项式()4P x 及三次样条函数()Q x （自然边界条件）对数据进行插 值。用图给出{（,i i x y ），i x =0.2+0.08i ，i=0，1,11,10}，()4P x 及()Q x 。 ①实验过程

x1=[0.2 0.4 0.6 0.8 1.0];

y1=[0.98 0.92 0.81 0.64 0.38]; n=length(y1); c=y1(:);

for j=2:n %求差商 for i=n:-1:j

c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1)); end end

syms x df d ;

df(1)=1;d(1)=y1(1);

for i=2:n %求牛顿差值多项式 df(i)=df(i-1)*(x-x1(i-1)); d(i)=c(i-1)*df(i); end

P4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数 pp=csape(x1,y1, 'variational');%调用三次样条函数 q=pp.coefs;

q1=q(1,:)*[(x-.2)^3;(x-.2)^2;(x-.2);1]; q1=vpa(collect(q1),5)

q2=q(1,:)*[(x-.4)^3;(x-.4)^2;(x-.4);1]; q2=vpa(collect(q2),5)

q3=q(1,:)*[(x-.6)^3;(x-.6)^2;(x-.6);1];

q3=vpa(collect(q3),5)

q4=q(1,:)*[(x-.8)^3;(x-.8)^2;(x-.8);1];

q4=vpa(collect(q4),5)%求解并化简多项式

②运行结果

P4 =

- 0.20833333333212067373096942901611*x^4 - 0.20833333333212067373096942901611*x^3 + 0.15833*x^2 + 0.96833*x + 0.782

q1 =

- 1.3392857142898719757795333862305*x^3 + 0.80357*x^2 - 0.40714285714057041332125663757324*x + 1.04

q2 =

- 1.3392857142898719757795333862305*x^3 + 1.6071*x^2 - 0.8892857142927823588252067565918*x + 1.1643

q3 =

- 1.3392857142898719757795333862305*x^3 + 2.4107*x^2 - 1.6928571428579743951559066772461*x + 1.4171

q4 =

- 1.3392857142898719757795333862305*x^3 + 3.2143*x^2 - 2.8178571428288705646991729736328*x + 1.8629

③问题结果

4次牛顿差值多项式()4P x = 0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784

三次样条差值多项式()Q x =

{2}2.在区间[-1，1]上分别取n=10.20用两组等距节点对龙格函数f(x)=2

2511

x

+)做多项式及三次样条差值，对每个n 值，分别画出差值函数及f(x)的函数. 龙格图函数做多项式

function L=lagrange(a,b,a0) syms x n=length(a) L=0.0 for i=1:n l=b(i); for j=1:i-1

l=l.*(x-a(j))/(a(i)-a(j)); end

for j=i+1:n

l=l.*(x-a(j))/(a(i)-a(j)); end

L=L+l;

simplify(L);

end

L=collect(L)

L=vpa(L,6)

L=subs(L,'x',a0);

end

clear all

subplot(1,2,1);

a=linspace(-1,1,10);

b=1./(1+25.*a.^2);

L=Lagrange(a,b)

b0=subs(L,'x',a);

plot(a,b0,'rs')

hold on

plot(a,b)

title('n=10ʱµÄ²åÖµº¯Êý'); subplot(1,2,2);

a=linspace(-1,1,20);

b=1./(1+25.*a.^2);

L=Lagrange(a,b)

b0=subs(L,'x',a);

plot(a,b0,'rs');

hold on;

plot(a,b)

title('n=20ʱµÄ²åÖµº¯Êý');图像：