今天给大家带来的是插值思想最经典的Lagrange插值,理论部分详见各种教科书,木木只负责提供代码,所演示的代码可直接拿走食用。
例:列车从直道进入到半径为R的圆弧弯道时,为安全行驶,往往需要进行一段缓冲区(缓和曲线哦,木木本科就是道路方向的哈哈哈),如图所示:
我国铁路大部分采用三次抛物线型缓和曲线。现测得某缓和曲线上4点不同的坐标如下表所示:
x坐标(m) | 0 | 60 | 120 | 180 |
y坐标(m) | 0 | 0.45 | 3.6 | 12.15 |
是利用Lagrange多项式插值法求缓和曲线的方程,并求x=115m时y的坐标
代码1:
function varargout = lagrange_interp(xdata,ydata,xi)
%拉格朗日插值
n=length(xdata); %求向量xdata的长度
if length(unique(xdata))<n %若输入的点出现相同时,给出错误提示
error('输入的点必须是互异的.')
end
L=zeros(n); %存储插值基函数
for i=1:n
px=poly(xdata([1:i-1,i+1:n]));%构造以x_j为根的多项式(j=1:i-1,i_1:n)
L(i,:)=px/polyval(px,xdata(i));%求插值基函数并存储
end
y=sum(bsxfun(@times,L,ydata(:)));%求插值多项式
if nargin==3 %若输入参数为3
y=polyval(y,xi); %根据插值多项式求指定点处的值
end
[varargout{1:2}]=deal(y,... %第一个输出参数为插值多项式或其在某点处的值
L); %第2个输出参数为插值基函数
命令行输入:
xdata=[0 60 120 180];
ydata=[0 0.45 3.6 12.15];
y=Lagrange_interp(xdata,ydata)
x=linspace(0,180);
f=polyval(y,x); %返回多项式的值
plot(x,f)
hold on
plot(xdata,ydata,'ko')
title(['缓和曲线方程:f(x)=0',strtrim(poly2str(y,'x'))])
绘图:
代码2:
function f=arui_lagrange(x0,y0,x)
%x0为节点变量,y0为节点上的函数值,x为插值点,f为返回插值
n=length(x0);m=length(x);
format long
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(x-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
f = s;
end
命令行输入:
xdata=[0 60 120 180];
ydata=[0 0.45 3.6 12.15];
y=Lagrange_interp(xdata,ydata,115)
f=arui_lagrange(xdata,ydata,115)
得出结果:
y= 3.168489583333333
f= 3.168489583333333