数值分析第一次实习题报告

clear;
x=linspace(-5,5,10);
y=1./(1+x.*x);
t=linspace(-5,5,500);
Y=1./(1+t.^(2));
xx = linspace(-5,5,10);    
y1 = interp1(x,y,xx,'linear');  
plot(x,y,'ro',t,Y,'g',xx,y1,'b');grid;

clear;
clc;
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)*df(i);
end
disp('4次牛顿插值多项式');
P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数
pp=csape(x1,y1, 'variational');%调用三次样条函数
q=pp.coefs;
disp('三次样条函数');
for i=1:4
    S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1];
    S=vpa(collect(S),5)
end
x2=0.2:0.08:1.08;
dot=[1 2 11 12]; 
figure
ezplot(P4,[0.2,1.08]);
hold on;
y2=fnval(pp,x2);
x=x2(dot);
y3=eval(P4);
y4=fnval(pp,x2(dot));
plot(x2,y2,'r',x2(dot),y4,'co');
title('4次牛顿插值及三次样条');

clear;
x=linspace(-1,1,10);
y=1./(1+25*x.*x);
t=linspace(-1,1,200);
Y=1./(1+25*t.^(2));
f=Language(x,y)
f =21.6248*t.^8 - 44.9155*t.^6 + 30.7285*t.^4 - 8.26092*t.^2 + 0.861538
S=interp1(x,y,t,'spline');
plot(x,y,'ro',t,Y,'g',t,f,'b',t,S,'g:');grid;

clear;
x=linspace(-1,1,20);
y=1./(1+25*x.*x);
t=linspace(-1,1,200);
Y=1./(1+25*t.^(2));
f=Language(x,y)
f=- 25671.2*t.^18 + 95604.8*t.^16 - 146791.0*t.^14 + 121019.0*t.^12 - 58585.0*t.^10 + 17172.5*t.^8 - 3055.32*t.^6 + 327.726*t.^4 - 21.6235*t.^2 + 0.992681
S=interp1(x,y,t,'spline');
plot(x,y,'ro',t,Y,'g',t,f,'b',t,S,'g:');grid;

（1） 8 次多项式插值：
首先建立新的 M-file:
输入如下代码（此为拉格朗日插值的功能函数） 并保存：
function f=Language(x,y,x0)
%求已知数据点的拉格朗日插值多项式
%已知数据点的 x 坐标向量： x
%已知数据点的 y 坐标向量： y
%插值的 x 坐标： x0
%求得的拉格朗日插值多项式或在 x0 处的插值： f
syms t;
if(length(x)==length(y))
n=length(x);
else
disp('x 和 y 的维数不相等!');
return;
end %检错
f=0.0;
for(i=1:n)
l=y(i);
for(j=1:i-1)
l=l*(t-x(j))/(x(i)-x(j));
end;
for(j=i+1:n)
l=l*(t-x(j))/(x(i)-x(j)); %计算拉格朗日基函数
end;
f=f+l; %计算拉格朗日插值函数
simplify(f); %化简
if(i==n)
if(nargin==3)
f=subs(f,'t',x0); %计算插值点的函数值
else
f=collect(f); %将插值多项式展开
f=vpa(f,6); %将插值多项式的系数化成 6 位精度的小数
end
end
end
再建立新的 M-file：
输入：
clear;
x=[0 1 4 9 16 25 36 49 64];
y=[0:1:8];
f=Language(x,y)
运行得到
f =
1.32574*t-.381410*t^2+.604294e-1*t^3+.222972e-3*t^5-.542921e-5*t^6+.671268e-7*t^7-.328063e-9*t^8-.498071e-2*t^4
这就是 8 次多项式插值
L_8(x)=1.32574*t-.381410*t^2+.604294e-1*t^3+.222972e-3*t^5-.542921e-5*t^6+.671268e-7*t^7-.328063e-9*t^8-.498071e-2*t^4.
(2)三次样条插值：
建立新的 M-file：
输入：
clear;
x=[0 1 4 9 16 25 36 49 64];
y=[0:8];
t=[0:0.1:64];
Y=t.^(0.5);
f=Language(x,y)
f=1.32574*t-.381410*t.^2+.604294e-1*t.^3+.222972e-3*t.^5-.542921e-5*t.^6+.671268e-7*t.^7-.328063e-9*t.^8-.498071e-2*t.^4;
S=interp1(x,y,t,'spline');
plot(x,y,'ro',t,Y,'r',t,f,'b',t,S,'g:');grid;
运行程序得到如下图：