今天给大家带来的是插值思想最经典的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('输入的点必须是互异的.')endL=zeros(n); %存储插值基函数for i=1:npx=poly(xdata([1:i-1,i+1:n]));%构造以x_j为根的多项式(j=1:i-1,i_1:n)L(i,:)=px/polyval(px,xdata(i));%求插值基函数并存储endy=sum(bsxfun(@times,L,ydata(:)));%求插值多项式if nargin==3 %若输入参数为3y=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 onplot(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 longs=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(x-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endf = 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