CODEWITHSUNDEEP
pythonnumpymatplotlibscipy
pappu
pappu

Reputation: 129

Fitting Differential Equations: curve_fit converges to local minima

I am trying to fit the differential equation ay' + by''=0 to a curve by varying a and b The following code does not work. The problem with curve_fit seems to be that lack of initial guess results in failure in fitting. I also tried leastsq. Can anyone suggest me other ways to fit such differential equation? If I don't have good guess curve_fit fails!

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!

Upvotes: 0

Views: 1526

Answers (2)

pv.
pv.

Reputation: 35125

If your problem is that the default initial guesses, read the documentation curve_fit to find out how to specify them manually by giving it the p0 parameter. For instance, curve_fit(deriv, time, zn, p0=(12, 0.23)) if you want a=12 and b=0.23 be the initial guess.

Upvotes: 0

ovgolovin
ovgolovin

Reputation: 13410

Let's rewrite the equation:

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

So this equation may be represented in this way

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

The code in 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()

You may compare the resultant plots with the plots in WolframAlpha

Upvotes: 2

Related Questions