%举例:
%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33; %f=InterpLagrange(x,y) %f=InterpLagrange(x,y,x0)
if length(x)==length(y) n=length(x); else
disp('节点个数和函数值个数不同!') f=' '; return; end p=0; for i=1:n l=y(i); for j=1:n if j==i continue; end
%利用卷积计算Lagrange基函数 l=conv(l,[1 -x(j)]./(x(i)-x(j))); end
%p是一向量,表示插值多项式的系数 p=p+l;
end
if nargin==3
f=polyval(p,x0);%计算插值多项式在x0处的值 else
f=poly2str(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式 end
function f=InterpLagrange2(x,y,x0) %构造Lagrange插值多项式
%此函数中借助符号运算来求Lagrange基函数,运算速度较慢,不推荐此种方法 %f=InterpLagrange2(x,y,x0)
%f:插值多项式或者是插值多项式在x0处的值 %x:节点 %y:函数值 %x0:某一测试点 %
%调用格式:
%f=InterpLagrange2(x,y) 返回插值多项式
%f=InterpLagrange2(x,y,x0) 返回插值多项式在点x0处的值 %举例:
%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33; %f=InterpLagrange2(x,y)
%f=InterpLagrange2(x,y,x0)
if length(x)==length(y) n=length(x); else
disp('节点个数和函数值个数不同!') f=' '; return; end syms t; f=0; for i=1:n l=y(i); for j=1:n if j==i continue; end
l=l*(t-x(j))/(x(i)-x(j));%借助符号运算,计算Lagrange基函数 end f=f+l;
simplify(f);%化简多项式 if i==n
if nargin==3
f=subs(f,'t',x0);%计算插值多项式f在点x0处的值 else
f=collect(f);%计算插值多项式,展开并合并同类项 f=vpa(f,6);%设置多项式系数的有效数字 end end end
4、Newdon插值法
function f=InterpNewdon(x,y,x0) %Newdon插值多项式 %f=InterpNewdon(x,y,x0)
%f:插值多项式或者是插值多项式在x0处的值 %x:节点 %y:函数值 %x0:某一测试点 %
%调用格式
%f=InterpNewdon(x,y) 返回插值多项式
%f=InterpNewdon(x,y,x0) 返回插值多项式在x0点的值 %应用举例:
%x=[1 2 3 4 5];y=[1 4 7 8 6];x0=6;
%f=InterpNewdon(x,y) %f=InterpNewdon(x,y,x0)
if length(x)==length(y) n=length(x); else
disp('节点个数和函数值个数不同!') f=' '; return; end
A=zeros(n);%初始化差商矩阵 for i=1:n
A(i,1)=y(i);%差商矩阵的第一列是函数值 end
%计算差商矩阵
%差商矩阵中对角线上的元素为Newdon插值多项式的系数 for j=2:n for i=j:n
A(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1)); end end
%求Newdon插值多项式 p=zeros(1,n);