I want to apply the operation curve_fit (to find coefficients c1-c7) to the integral function
The problem is that the subintegral function depends on several variables(lists defor,stress,init_def) and the integration limits are variable too (lists init_def,defor).
So, the integral function is time=f(defor,stress)
The integration should be over the parameter (defor)
I tried following code, but errors occur. I would be very appreciate if you could help me. Thank you very much in advance.
import numpy as np
from scipy import integrate
from scipy.optimize import curve_fit
from scipy.integrate import quad
time =np.array([0, 18, 24, 42, 48, 66, 72, 0, 4, 22, 28, 46, 52, 70])
defor =np.array([0.11, 0.62, 0.73, 0.91, 1.0, 1.17, 1.22, 0.15, 0.26, 0.4, 0.43, 0.51, 0.51, 0.58])
init_def =np.array([0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15])
stress =np.array([10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0,7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5])
temp=943
e=2.71828
arg=zip(defor,stress,init_def)
# subint_funct- enter subintegral function with с1-с7 uncertain coefficients
def subint_funct(arg, c1,c2, c3,c4,c5,c6,c7):
return ((c1**(-1))*(temp**c2)*(stress**(-c3))*e**((c5-c6*stress)/temp)*
((defor+1)**(-c3))*(defor**c4)*e**(-((c6*stress+c7)*defor/temp)) )
# integration- integration of subint_funct with lower and upper integration bounds init_def and defor
def integration(arg,c1,c2,c3,c4,c5,c6,c7):
return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]
# curve fitting, dependent data-time
parameter = sp.curve_fit(integration,
arg,
time)
parameter=parameter[0]
print (parameter)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [11], in <cell line: 25>()
22 def integration(arg,c1,c2,c3,c4,c5,c6,c7):
23 return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]
---> 25 parameter = sp.curve_fit(integration,
26 arg,
27 time)
28 parameter=parameter[0]
29 print (parameter)
File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:789, in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
787 # Remove full_output from kwargs, otherwise we're passing it in twice.
788 return_full = kwargs.pop('full_output', False)
--> 789 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
790 popt, pcov, infodict, errmsg, ier = res
791 ysize = len(infodict['fvec'])
File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:410, in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
408 if not isinstance(args, tuple):
409 args = (args,)
--> 410 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
411 m = shape[0]
413 if n > m:
File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:24, in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
22 def _check_func(checker, argname, thefunc, x0, args, numinputs,
23 output_shape=None):
---> 24 res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
25 if (output_shape is not None) and (shape(res) != output_shape):
26 if (output_shape[0] != 1):
File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:485, in _wrap_func.<locals>.func_wrapped(params)
484 def func_wrapped(params):
--> 485 return func(xdata, *params) - ydata
Input In [11], in integration(arg, c1, c2, c3, c4, c5, c6, c7)
22 def integration(arg,c1,c2,c3,c4,c5,c6,c7):
---> 23 return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]
Input In [11], in <listcomp>(.0)
22 def integration(arg,c1,c2,c3,c4,c5,c6,c7):
---> 23 return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]
File ~\anaconda3\lib\site-packages\scipy\integrate\quadpack.py:351, in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
348 flip, a, b = b < a, min(a, b), max(a, b)
350 if weight is None:
--> 351 retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
352 points)
353 else:
354 if points is not None:
File ~\anaconda3\lib\site-packages\scipy\integrate\quadpack.py:463, in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
461 if points is None:
462 if infbounds == 0:
--> 463 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
464 else:
465 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
TypeError: only size-1 arrays can be converted to Python scalars
As I am not 100% sure about the details of the OP, here an example that does something similar. Simple modifications should do the trick then.
import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit
### not clever, but whatever...
def integ( arg, c1, c2, c3 ):
x, k1 = arg
return k1**c1 * x**c2 + c3
def integ_first_wrapper( x, k, c):
return(integ( ( x, k ), *c ) )
def fitwrapper( xlist, c1, c2, c3 ):
if len( xlist.shape ) == 2: ##array
out = np.fromiter(
( fitwrapper( x, c1, c2, c3 ) for x in xlist ),
np.float
)
else:
l, u, s = xlist
out = quad(
integ_first_wrapper, l, u, args=( s, ( c1, c2, c3 ) )
)[0]
return out
c10 = 2.33
c20 = 3.11
c30 = 4.11
lolimits = 5 * np.random.random( size=10 )
uplimits = 10 + 5 * np.random.random( size=10 )
stress = 28 * np.random.random( size=10 )
xin = np.transpose( ( lolimits, uplimits, stress ) )
### create test data
data = list()
for l,u, s in zip( lolimits, uplimits, stress ):
# ~data.append( quad( integ_first_wrapper, l, u, args=( s, (c10, c20, c30 ) ) )[0] ) ### without error
data.append(
( 1 + 0.01 * ( 1 - 2 * np.random.random() ) )
* quad(
integ_first_wrapper, l, u, args=( s, (c10, c20, c30 ) )
)[0]
) ### wit relative error
# ~print(xin)
res, err = curve_fit( fitwrapper, xin, data, p0=( 1, 1, 1 ) )
print( res )
print( err )