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