分段线性插值与三次样条插值法.doc

实验报告
实验项目
插值法
实验日期
2016/9/30
理论内容
分段线性插值与三次样条插值
授课日期
2016/9/0
实验室名称
文理管203
微机编号
E1
实验目的及要求：
1、 学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；
2、 明确插值多项式和分段插值多项式各自的优缺点；
3、 熟悉插值方法的程序编制；如果绘出插值函数的曲线，观察其光滑性。
实验内容：
编写分段线性插值法及三次样条插值法通用子程序，依据数据表

2.0
2.1
2.2
2.3
2.4

1.414214
1.449138
1.483340
1.516575
1.549193
构造相应的插值多项式，并计算函数在的近似值。

实验步骤及程序：
1、分段线性插值法流程图
是
否
是
输入xi,yi,i=1,⋯,n及x
ia⟵1,ib⟵n
ib-ia>1?
im⟵ia+ib2的整数部分
xm2.4];
y=[1.414214 1.449138 1.483340 1.516575 1.549193];
y_1=0.5*y.^(-0.5);
x0=2.15;
f = 0.0;
if(length(x) == length(y))
if(length(y) == length(y_1))
n = length(x);
else
disp('y和y的导数的维数不相等！');
return;
end
else
disp('x和y的维数不相等！');
return;
end
for i=1:n
if(x(i)<=x0)&& (x(i+1)>=x0)
index = i;
break;
end
end
h = x(index+1) - x(index);
fl = y(index)*(1+2*(x0-x(index))/h)*(x0-x(index+1))^2/h/h + ...
y(index+1)*(1-2*(x0-x(index+1))/h)*(x0-x(index))^2/h/h;
f = fl;

3、三次样条插值法流程图

是
输入数据flag
flag=1？
flag=2？
flag=3？
输出失败信息
执行子程序“第一类边界下的样条三次插值”
执行子程序“第二类边界下的样条三次插值”
执行子程序“第三类边界下的样条三次插值”
是
是
否
否

否

4、三次样条插值法源程序
function[yy, b, c, d]=spline3(~,~,~,~,~,~)
%三次样条插值函数 %(x,y)为插值节点，xx为插值点；
%flag表端点边界条件类型：
%flag=0:自然样条(端点二阶导数为0)；
%flag=1:第一类边界条件(端点一阶导数给定)；
%flag=2:第二类边界条件(端点二阶导数给定)；
%vl,vr表左右端点处的在边界条件值。
%样条函数为：Si(x)=yi+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3
%b,c,d分别为各子区间上的系数值
%yy表插值点处的函数值.
x=[2.0 2.1 2.2 2.3 2.4];
y=[1.414214 1.449138 1.483340 1.516575 1.549193];
xx=2.15;
flag=0;
vl=2.0;
vr=2.4;
if length(x)length(y)
n=length(x);
a=zeros(n-1,1);
b=a;
d=a;
dx=a;
dy=a;
A=zeros(n);
B=zeros(n,1);
end
for i=1:n-1
a(i)=y(i);
dx(i)=x(i+1)-x(i);
dy(i)=y(i+1)-y(i);
end
for i=2:n-1
A(i,i-1)=dx(i-1);
A(i,i)=2*(dx(i-1)+dx(i));
A(i,i+1)=dx(i);
B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));
end
%自然样条端点条件(端点二阶导数为零)%
if flag==0;
A(1,1)=1;
A(n,n)=1;
end
%---------------------------------%
%端点一阶导数条件%
if flag==1
A(1,1)=2*dx(1);
A(1,2)=dx(1);
A(n,n-1)=dx(n-1);
A(n,n)=2*dx(n-1);
B(1,1)=3*(dy(1)/dx(1)-vl);
B(n,1)=3*(vr-dy(n-1)/dx(n-1))
end
%---------------%
%端点二阶导数条件%
if flag==2
A(1,1)=2;
A(n,n)=2;
B(1,1)=vl;
B(n,1)=vr;
end
%---------------%
c=A\B;
for i=1:n-1
d(i)=(c(i+1)-c(i))/(3*dx(i));
b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3;
end
[mm,nn]=size(xx);
yy=zeros(mm,nn);
for i=1:mm*nn
for ii=1:n-1
if xx(i)>=x(ii) && xx(i)