问题是你最初的猜测与实际的解决方案相去甚远。如果在chisqfunc()中添加一个print语句,比如print (a,b),然后重新运行代码,您将得到如下结果:(0, 0)
(1.4901161193847656e-08, 0.0)
(0.0, 1.4901161193847656e-08)
这意味着minimize只在这些点计算函数。
如果您现在尝试在这3对值上计算chisqfunc(),您将看到它们完全匹配,例如print chisqfunc((0,0))==chisqfunc((1.4901161193847656e-08,0))
True
这是因为浮点运算的舍入。换句话说,在计算stress - model时,varstress的数量级太多,大于model,结果被截断。
然后可以尝试对其进行粗处理,增加浮点精度,在用loadtxt加载数据之后写入data=data.astype(np.float128)。minimize失败,有result.success=False,但有一条有用的消息Desired error not necessarily achieved due to precision loss.
一种可能是提供一个更好的初始猜测,以便在减法中stress - model的model部分具有相同的数量级,另一种可能是重新调整数据,以便解决方案更接近您的初始猜测(0,0)。
如果你只是重新缩放数据,例如使某个应力值无量纲(比如这种材料的嘶嘶声/破裂声),会好得多
这是拟合的一个例子,使用最大测量应力作为应力标度。代码中很少有更改:import numpy
import scipy.optimize as opt
filename = 'data.csv'
data = numpy.loadtxt(open(filename,"r"),delimiter=",")
stress = data[:,0]
strain = data[:,1]
err_stress = data[:,2]
smax = stress.max()
stress = stress/smax
#I am assuming the errors err_stress are in the same units of stress.
err_stress = err_stress/smax
def chisqfunc((a, b)):
model = a + b*strain
chisq = numpy.sum(((stress - model)/err_stress)**2)
return chisq
x0 = numpy.array([0,0])
result = opt.minimize(chisqfunc, x0)
print result
assert result.success==True
a,b=result.x*smax
plot(strain,stress*smax)
plot(strain,a+b*strain)
你的线性模型很好,也就是说,你的材料在这个变形范围内有一个非常线性的行为(它到底是什么材料?):