scipy.curve_fit 返回的不合逻辑的参数

Question

我正在用 Python 模拟一个球在流体中落下，并使用阻尼系数（a 和b）和流体密度将模型函数拟合到一组数据点，但拟合值为流体密度返回负值，我不知道代码有什么问题。我的代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import curve_fit


##%%Parameters and constants 
m       = 0.1                               #mass of object in kg
g       = 9.81                              #acceleration due to gravity in m/s**2
rho     = 700                               #density of object in kg/m**3
v0      = 0                                 #velocity at t=0
y0      = 0                                 #position at t=0
V       = m / rho                           #volume in cubic meters
r       = ((3/4)*(V/np.pi))**(1/3)          #radius of the sphere
asample = 0.0001                            #sample value for a 
bsample = 0.0001                            #sample value for b

#%%Defining integrating function 
##fuction => y'' = g*(1-(rhof/rho))-((a/m)y'+(b/m)y'**2)
##   y' = v
##   v' = g*(1-rhof/rho)-((a/m)v+(b/m)v**2)

def sinkingball(Y, time, a, b, rhof):
    return [Y[1],(1/m)*(V*g*(rho-rhof)-a*Y[1]-b*(Y[1]**2))]

def balldepth(time, a, b, rhof):
    solutions = odeint(sinkingball, [y0,v0], time, args=(a, b, rhof))
    return solutions[:,0] 

time    = np.linspace(0,15,151)

# imported some experimental values and named the array data

a, b, rhof = curve_fit(balldepth, time, data, p0=(asample, bsample, 100))[0]
print(a,b,rhof) 

Answer 1

提供你实际得到的输出会很有帮助，关于time 未被sinkingball() 使用的评论值得关注。

您可能会发现 lmfit (https://lmfit.github.io/lmfit-py) 很有用。这为曲线拟合提供了更高级别的接口，除其他外，它允许对参数设置界限，以便它们可以保持物理上的合理性。我认为您的问题会从 curve_fit 转换为 lmfit 为：

from lmfit import Model
def balldepth(time, a, b, rhof):
    solutions = odeint(sinkingball, [y0,v0], time, args=(a, b, rhof))
    return solutions[:,0] 

# create a model based on the model function "balldepth"
ballmodel = Model(balldepth)

# create parameters, which will be named using the names of the
# function arguments, and provide initial values
params = bollmodel.make_params(a=0.001, b=0.001, rhof=100)

# you wanted rhof to **not** vary in the fit:
params['rhof'].vary = False

# set min/max values on `a` and `b`:
params['a'].min = 0
params['b'].min = 0

# run the fit
result = ballmodel.fit(data, params, time=time)

# print out full report of results
print(result.fit_report())

# get / print out best-fit parameters:
for parname, param in result.params.items():
    print("%s = %f +/- %f" % (parname, param.value, param.stderr))