0

I am currently using "scipy.optimize.curve_fit"

In the equation (piecewise) I fit with, "theta" is calculated by one of the params ("D") first, then the choosing of equations is dependent on "theta" (therefore, dependent on "D"). However, it seems that scipy input the params as a list or something, so it does not choose the right equation for each "D" individually. Instead, it asks me to use "theta.any()" or "theta.all()" which is not what I want.

Is there any module that can do what I want?

Many thanks,

Christopher

import numpy as np
from scipy.optimize import curve_fit
def MahonOldham(t,D:'m^2/s',c:'mM'):
    n=1
    F=96485     #A s / mol
    pi=3.14
    a=12.5*10**-6           # m
    D*=10**-10
    theta=D*t/a**2
    if theta <=1.281: #HERE!!!
        factor=1/(np.sqrt(pi*theta))+1+np.sqrt(theta/(4*pi))-3*theta/25+(3*theta**(3/2))/226
        I=n*pi*F*c*D*a*factor #nA
    else:
        factor=4/pi+8/np.sqrt(pi**5*theta)+25*theta**(-3/2)/2792-theta**(-5/2)/3880-theta**(-7/2)/4500
        I=n*pi*F*c*D*a*factor #nA
        
    return I*10**9
sylvanas
  • 1
  • 2
  • You can try [lmfit][1]. It provides an interface to express parameters as a function of other parameters. [1]: https://lmfit.github.io/lmfit-py/index.html – mss Jul 04 '22 at 10:25
  • Check also this: https://stackoverflow.com/questions/41641880/using-scipy-curve-fit-with-piecewise-function?rq=1 – mss Jul 04 '22 at 10:47

1 Answers1

0

This is the module that you want, you just need to adapt to it's API. The function curve_fit passes to your function arrays for D and c rather than single values, so you need to vectorize MahonOldham (and, please, rename to mahon_oldham, PEP8).

In particular replace if statement by np.where:

I = np.where(
    theta <= 1.281, 
    1/(np.sqrt..., # first case
    4/pi+8/np(..., # second case
)
psarka
  • 1,562
  • 1
  • 13
  • 25
  • Thank you, Now i have another problem, lol. my equation is like f(t)/(t+dt), where t is the independent var and dt is a param to fit. my dataset included t=0 and dt=0 is also allowed. Therefore when t=dt=0, the result is not finite, and I think that's what raised the error: "ValueError: Residuals are not finite in the initial point." is there anyway to solve this? Thank you!!! – sylvanas Jul 07 '22 at 12:48
  • @sylvanas Not sure, but try to provide initial point that is not t=dt=0. – psarka Jul 07 '22 at 15:11