Search code examples
pythonnumpyscipycurve-fittingscipy-optimize

Scipy.Curve fit struggling with exponential function


I'm trying to fit a curve of the equation:

y = ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2

where  A = (-np.log((k1/v)*Co))/k2

given to me by a supervisor to a dataset that looks like a rough exponential that flattens to a straight horizontal line at its top. When I fit the equation i am receiving only a straight line from the curve fit and a corresponding Warning:

<ipython-input-24-7e57039f2862>:36: RuntimeWarning: overflow encountered in exp
 return ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2

the code I am using looks like:

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

xData_forfit = [1.07683e+13, 1.16162e+13, 1.24611e+13, 1.31921e+13, 1.40400e+13, 2.65830e+13,
 2.79396e+13, 2.86676e+13, 2.95155e+13, 3.03605e+13, 3.12055e+13, 3.20534e+13,
 3.27814e+13, 3.36293e+13, 3.44772e+13, 3.53251e+13, 3.61730e+13, 3.77459e+13,
 3.85909e+13, 3.94388e+13, 4.02838e+13, 4.11317e+13, 4.19767e+13, 4.27076e+13,
 5.52477e+13, 5.64143e+13, 5.72622e+13, 5.81071e+13, 5.89550e+13, 5.98000e+13,
 6.05280e+13, 6.13759e+13, 6.22209e+13, 6.30658e+13, 6.39137e+13, 6.46418e+13,
 6.55101e+13, 6.63551e+13, 6.72030e+13, 6.80480e+13, 6.88929e+13, 6.97408e+13,
 7.04688e+13, 7.13167e+13, 7.21617e+13, 8.50497e+13, 8.58947e+13, 8.67426e+13,
 8.75876e+13, 8.83185e+13, 9.00114e+13, 9.08563e+13, 9.17013e+13]
yData_forfit = [1375.409524, 1378.095238, 1412.552381, 1382.904762, 1495.2, 1352.4,
 1907.971429, 1953.52381,  1857.352381, 1873.990476, 1925.114286, 1957.085714,
 2030.52381,  1989.8,      2042.733333, 2060.095238, 2134.361905, 2200.742857,
 2342.72381,  2456.047619, 2604.542857, 2707.971429 ,2759.87619,  2880.52381,
 3009.590476, 3118.771429, 3051.52381,  3019.771429, 3003.561905, 3083.0,
 3082.885714, 2799.866667, 3012.419048, 3013.266667, 3106.714286, 3090.47619,
 3216.638095, 3108.447619, 3199.304762, 3154.257143, 3112.419048, 3284.066667,
 3185.942857, 3157.380952, 3158.47619,  3464.257143, 3434.67619,  3291.457143,
 2851.371429, 3251.904762, 3056.152381, 3455.07619,  3386.942857]

def fnct_to_opt(t, k2, k1):
    #EXPERIMENTAL CONSTANTS
    v = 105
    Co = 1500

    A = (-np.log((k1/v)*Co))/k2
    return ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2


initial_k2k1 = [100, 1*10**-3]
constants = curve_fit(fnct_to_opt, xData_forfit, yData_forfit, p0=initial_k2k1)
k2_fit = constants[0][0]
k1_fit = constants[0][1]


fit = []
for i in xData_forfit:
    fit.append(fnct_to_opt(i,k2_fit,k1_fit))


plt.plot(xData_forfit, yData_forfit, 'or', ms='2')
plt.plot(xData_forfit, fit)

this is giving me this plot as a result: enter image description here

As far as i can tell, the code isn't producing a useful output due to a too large value for the np.exp term, but i don't know how to go about diagnosing where this overflow is coming from or how to fix the issue. any help would be appreciated, thanks.


Solution

  • The overflow is happening exactly where the error message tells you: in the return expression of fnct_to_opt. I asked you to print the offending values just before the error point; this would show you the problem.

    At the point of error, the values in A are in the range e+13 to e+14. t is insignificant; k2 is a bit under -10000.0

    Thus, the values in your argument to np.exp are well out of the domain that the function can handle. Just add a line to you function and watch the results:

    def fnct_to_opt(t, k2, k1):
        #EXPERIMENTAL CONSTANTS
        v = 105
        Co = 1500
    
        A = (-np.log((k1/v)*Co))/k2
        print("TRACE", "\nk2", k2, "\nt", t, "\nA", A, "\nother", k1, v, Co)
        return ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2