多项式插值和三次样条插值.docx

已知某产品从1900年到2010年每隔10年的产量,用多项式插值和三次样条插值的方法，画出每隔一年的插值曲线的图形,试计算并比较在不同方法下的2005年以及2015年的产量。
年份190019101920193019401950产量75.99591.972105.711123.203131.699150.697年份196019701980199020002010产量179.323203.212226.505251.525291.854325.433思想算法：多项式插值采用牛顿多项式插值法，该算法可以克服多项式插值和拉格朗日插值法的缺点，即：当用已知的n+1个数据点求出插值多项式后，又获得了新的数据点，要用它连同原有的n+1个数据点一起求出插值多项式，从原已计算出的n次插值多项式计算出新的n+1次插值多项式是很困难的，必须全部重新计算。而牛顿插值法可以克服这一缺点。
三次样条插值不仅在节点增多使子区间减少时，误差随之减少，也使曲线具有足够的光滑性。
Matlab程序如下
程序一：牛顿插值法
源程序名称Newton.m
clearall;
x=0:10:110;
y=[75.995,91.972,105.711,123.203,131.699,150.697,179.323,203.212,226.505,251.525,291.854,325.433];
n=length(x);symst;
fork=2:n
fori=n:-1:k
y(i)=(y(i)-y(i-1))/(x(i)-x(i-k+1));
end;
end;
N=y(1);
fori=2:n
W=1;
forj=1:(i-1)
W=W*(t-x(j));
end;
N=N+y(i)*W;
end;
N=expand(N);
ezplot(N,[0,120]);
holdon;
formatshort;
Q=[];
fori=0:120
Q(i+1)=subs(N,t,i);
end;
T=0:120;
plot(T,Q,'^');
title('产量随时间变化曲线');
xlabel('T/时间');
ylabel('Q/产量');
N105=subs(N,t,105);
N115=subs(N,t,115);
程序二：三次样条插值
源程序名称SPLINEM.m
clearall;
x=0:10:110;
y=[75.995,91.972,105.711,123.203,131.699,150.697,179.323,203.212,226.505,251.525,291.854,325.433];
n=length(x);symst;
fori=1:n
p(i)=y(i);
end;
fork=2:3
fori=n:-1:k
p(i)=(p(i)-p(i-1))/(x(i)-x(i+1-k));
end;
end;
h(2)=x(2)-x(1);
fori=2:(n-1)
h(i+1)=x(i+1)-x(i);
c(i)=h(i+1)/(h(i)+h(i+1));
a(i)=1-c(i);
b(i)=2;
p(i)=6*p(i+1);
end;
p(1)=0;p(n)=0;c(1)=0;b(1)=2;b(n)=2;a(n)=0;
u(1)=b(1);z(1)=p(1);
fork=2:n
l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);
z(k)=p(k)-l(k)*z(k-1);
end;
M(n)=z(n)/u(n);
fork=(n-1):-1:1
M(k)=(z(k)-c(k)*M(k+1))/u(k);
end;
fori=2:n
S(i)=M(i-1)*(x(i)-t)^3/(6*h(i))+M(i)*(t-x(i-1))^3/(6*h(i))+(y(i-1)-M(i-1)*h(i)^2/6)*(x(i)-t)/h(i)+(y(i)-M(i)*h(i)^2/6)*(t-x(i-1))/h(i);
p=0;
fork=((i-2)*10+1):((i-1)*10)
Q(k)=subs(S(i),t,x(i-1)+p);p=p+1;
end;

ezplot(S(i),[x(i-1),x(i)]);
holdon;
end;
ezplot(S(12),[110,120]);
holdon;
fork=111:121
Q(k)=subs(S(12),t,k-1);
end;
T=0:120;
plot(T,Q,'^');
axis([0,120,0,400]);
title('产量随时间变化曲线');
xlabel('T/时间');
ylabel('Q/产量');
S105=subs(S(12),t,105);
S115=