CODEWITHSUNDEEP
pythonscipyscipy-optimize
horlust
horlust

Reputation: 111

Why is curve_fit giving me different results in different computers?

I have the following model:

-dc/dt = kc(t) - k'n(t)

I compute c(t) and n(t) from a simulation trajectory with CuPy. I then compute k and k' with curve_fit, through the following code:

# compute dc/dt
dc_dt = np.gradient(c_t, times)
    
# define the model function
def model(t, k, k_prime):
    return -k * np.array(c_t) + k_prime * np.array(n_t)
    
initial_guess = [0.5, 0.5]
bounds = ([0.0, 0.0], [np.inf, np.inf])
    
# fit the model to the data
params, covariance = curve_fit(model, times, dc_dt)
k1_optimized, k2_optimized = params

However, for all temperatures, k' = k2 = 1.0, which is incorrect.

I tried saving the c(t) and n(t) results to an Excel sheet and computing k1 and k2 in another computer; in this case, k2 wasn't 1.0. I was wondering what is going wrong?

Upvotes: 1

Views: 87

Answers (1)

lastchance
lastchance

Reputation: 6745

You have a standard least-squares problem, similar to that which would give you linear regression.

Let c, dc/dt and n be arrays. Form the squared error

enter image description here

Partial-differentiate with respect to k and k’ and set these equal to 0 (aiming to minimize S2).

enter image description here

These give two simultaneous linear equations:

enter image description here

which can be solved by elimination for k and k’.

Code:

import numpy as np

# Generate some data
N = 1000
times = np.linspace( 0, 1, N )
c_0 = 1
k0, k0_prime = 2, 3
n_t = 4 * np.ones( N )
c_t = c_0 * np.exp( -k0 * times ) + ( k0_prime * n_t / k0 ) * ( 1 - np.exp( -k0 * times ) )

dc_dt = np.gradient(c_t, times[1] - times[0] )

# Form relevant sums
Scc = ( c_t * c_t ).sum()
Snn = ( n_t * n_t ).sum()
Scn = ( c_t * n_t ).sum()
Scd = ( c_t * dc_dt ).sum()
Snd = ( n_t * dc_dt ).sum()
k = -( Snn * Scd  - Scn * Snd ) / ( Snn * Scc - Scn * Scn )
kprime = ( Scc * k + Scd ) / Scn
print( f'k = {k:6.3f}, kprime = {kprime:6.3f}' )

Output:

k =  2.000, kprime =  3.000

Upvotes: 2

Related Questions