matlab找多项式最值,[转载]Matlab Polyfit高阶多项式插值取最小值

昨晚在看Lammps前辈 Michael的博客，看到他《MATLAB计算平衡晶格常数》的博文，今早拿着Lammps的脚本算出来，后用他给的脚本在Matlab中运行，结果提示：

Warning: Polynomial is badly conditioned. Add points with

distinct X values, reduce the degree of the

polynomial, or try centering and scaling as described in HELP

POLYFIT.

他博文中采用了8阶多项式插值，我们可以看到Lammps循环输出的数据点>9，因此多项式插值的方程数条件是满足的，那么问题出在哪呢？随后我将数据导入到cftool中，在没用smooth的前提下，发现高阶多项式插值效果都还不错。也就是说cftool的polyN和polyfit的算法还是不同的，同时高阶多项式插值是可以得到插值系数的。

多番尝试后，我将x轴的数据做了线性处理：x=

(x-mean)./std(x)后，发现polyfit居然可以工作了。。。究其原因应该是在这套数据体系下，x的差别特别小，使得求解插值系数的矩阵病态。

这是修改后的代码：

function cal_lat_const(N,inFileName)

% calcuate the lattice constant according to "lat_const

cohesive_energy"

% Input:

% N:

order of the polynomial fitting.

%

inFileNmae: name of the file storing "lat_const

cohesive_energy"

% Example:

%

cal_lat_const(8,'Au')

%

Here, 'Au' is a file stores "lat_const cohesive_energy"

% Powered by Xianbao Duan

% Email: xianbao.d@gmail.com

% Website: http://www.52souji.net/

%

----------------------------------------------------------------------

% Add centering and scaling with X values by mxio.

% Website: http://blog.sina.com.cn/mxio

%

______________________________________________________________________

% read in data from the file

data = load(inFileName,'-ascii');

x = data(:,1);

y = data(:,2);

[Ecoh,S,mu]

= polyfit(x,y,N);

% polynomial fitting and scale X values

dEcoh = polyder(Ecoh);

% derivation of the polynomial equation

zero_points = roots(dEcoh)*mu(2)+mu(1);

% solve the zero points and rescale

% calculate the lattice constant and coresponding cohesive

energy

for i = 1:

length(zero_points)

if isreal(zero_points(i))

if zero_points(i)

> x(1)

if zero_points(i)

< x(end)

lat_const = zero_points(i)

coh_energy = spline(x,y,lat_const)

end

end

end

end

% plot the curve

xx = x(1):0.01:x(end);

yy = polyval(Ecoh,xx,[],mu);

figure;

plot(x,y,'o',xx,yy);

xlabel('lattice

constant / A');

ylabel('cohesive

energy / eV');

legend('Original

points','Fitted

curve');

mxio

2015.3.24