数值分析实验(2)

实验二 插值法 P50

专业班级：信计131班 姓名：段雨博 学号：2013014907 一、实验目的

1、熟悉MATLAB 编程；

2、学习插值方法及程序设计算法。 二、实验题目

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 及三次样条函数()S x （自然边界条件）对数据进行插值用图给出

(){},,0.20.08,0,1,11,10i

i

i

x y x i i =+=，()4

P x 及()S x 。

2、在区间[]1,1-上分别取10,20n =用两组等距节点对龙格函数()2

1

125f x x

=

+作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及()f x 的图形。 3、下列数据点的插值 x 0 1 4 9 16 25 36 49 64 y 0 1

2

3

4

5

6

7

8

可以得到平方根函数的近似，在区间[]0,64上作图 （1）用这9个点作8次多项式插值()8L x （2）用三次样条（第一边界条件）程序求()S x

从得到结果看在[]0,64上，哪个插值更精确；在区间[]0,1上，两种插值哪个更精确？ 三、实验原理与理论基础 1、拉格朗日差值公式

)()(111k k

k k

k k x x x x y y y x L ---+

=++ 点斜式

k

k k

k k k k k x x x x y x x x x y x L --+--=++++11

111)( 两点式

2、n 次插值基函数 ....,2,1,0,)()(0

n j y x l y x L i

j

n

k k

k j n ===

∑=

n k x x x x x x x x x x x x x l n k n k k k k k ,...,1,0,)

()

(...

)

()

(...

)

()

()(1100=------=

--

3、牛顿插值多项式

...))(](,,[)](,[)()(102100100+--+++=x x x x x x x f x x x x f x f x P n ))...(](,...,[100---+n n x x x x x x f

)(],...,,[)()()(10x x x x f x P x f x R n n n n +=-=ω

4、三次样条函数

若函数],,[)(2

b a C x S ∈且在每个小区间],[1+j j x x 上是三次多项式，其中，

b x x x a n =<<<=...10是给定节点，则称)(x S 是节点n x x x ,...,,10上的三次样条函数。

若在节点j x 上给定函数值),,...,2,1,0)((n j x f y j i ==并成立,,...,2,1,0,)(n j y x S i j ==则称)(x S 为三次样条插值函数。 5、三次样条函数的边界条件

（1）0)()('

'''''00''====n n f x S f x S

（2）'

''00')(,)(n n f x S f x S == 四、实验内容 1、M 文件：

function [p]=Newton_Polyfit(X,Y)

format long g r=size(X); n=r(2); M=ones(n,n); M(:,1)=Y'; for i=2:n for j=i:n

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

p0=[zeros(1,n-1) M(1,1)];p=p0; for i=1:n-1

p1=M(i+1,i+1).*poly(X(1:i)); p0=[zeros(1,n-i-1) p1]; p=p+p0; end

3、M 文件：

function f=Language(~,~,~)

%%求已知数据点的拉格朗日插值多项式 %%已知数据点的x 坐标向量：x %%已知数据点的y 坐标向量：y %%插值点的x 坐标：x0 %%解得的拉格朗日插值多项式f

x=[0.0 0.4 0.8 1.2 1.6];

y=[0 0.428392 0.722101 0.910314 0.970348]; x0=[0.3 0.5];

syms t l;

if(length(x)==length(y))

n=length(x);

else

disp('x，y维数不一样');

return;

end

p=sym(0);

for i=1:n

l=sym(y(i));

for k=1:i-1

l=l*(t-x(k))/(x(i)-x(k));

end

for k=i+1:n

l=l*(t-x(k))/(x(i)-x(k));

end

p=p+1;

end

simplify(p);

f=subs(p,'t',x0);

f=vpa(f,6);

end

五、实验结果

1、

>> X=[0.2 0.4 0.6 0.8 1.0];

>> Y=[0.98 0.92 0.81 0.64 0.38];

>> [p]=Newton_Polyfit(X,Y);

>> Y2=polyval(p,X);

>> X1=0:0.01;1;

>> Y3=interp1(X,Y,X1,'spline');

>> plot(X,Y,'o',X,Y2,'r',X1,Y3,'g') 图像：