拟合微分方程：curve_fit收敛到局部最小值

Question

我正在尝试通过调整参数a和b，把一个微分方程ay' + by''=0拟合到一条曲线上。但是我写的代码没有成功。问题出在curve_fit这个函数上，因为没有一个好的初始猜测，它就无法进行拟合。我也试过leastsq这个方法。有没有人能给我推荐其他的方法来拟合这样的微分方程？如果我没有一个好的初始猜测，curve_fit就会失败！

from scipy.integrate import odeint
from scipy.optimize import curve_fit
from numpy import linspace, random, array

time = linspace(0.0,10.0,100)
def deriv(time,a,b): 
    dy=lambda y,t : array([ y[1], a*y[0]+b*y[1] ])
    yinit = array([0.0005,0.2]) # initial values
    Y=odeint(dy,yinit,time)
    return Y[:,0]

z = deriv(time, 2, 0.1)
zn = z + 0.1*random.normal(size=len(time))

popt, pcov = curve_fit(deriv, time, zn)
print popt  # it only outputs the initial values of a, b!

Answer 1

如果你的问题是默认的初始猜测不合适，可以查看文档 curve_fit，了解如何手动指定这些初始值，方法是给它传递 p0 参数。例如，如果你希望初始猜测为 a=12 和 b=0.23，可以这样写：curve_fit(deriv, time, zn, p0=(12, 0.23))。

Answer 2

我们来重新写一下这个方程：

ay' + by''=0
y'' = -a/b*y'

所以这个方程可以这样表示：

dy/dt = y'
d(y')/dt = -a/b*y'

下面是用Python 2.7写的代码：

from scipy.integrate import odeint
from pylab import *

a = -2
b = -0.1

def deriv(Y,t):
    '''Get the derivatives of Y at the time moment t
Y = [y, y' ]'''
    return array([ Y[1], -a/b*Y[1] ])

time = linspace(0.0,1.0,1000)
yinit = array([0.0005,0.2]) # initial values
y = odeint(deriv,yinit,time)
figure()
plot(time,y[:,0])
xlabel('t')
ylabel('y')
show()

你可以把得到的图和WolframAlpha上的图进行比较。