昨晚在看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) <</SPAN> 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