我们从Python开源项目中，提取了以下22个代码示例，用于说明如何使用scipy.optimize.fsolve()。
def control(self, tau): """ I referred to robot modeling and control when making this program :param deltaf: :return: author: weiwei date: 20170417 """ func = lambda th: self.__linkI*(th-self.__thnp-self.__thnp+self.__thnpp)/(self.__deltat*self.__deltat)+\ self.__dampc*(th-self.__thnp)/self.__deltat+self.__mgl*math.sin(th)-tau # initial guess th_ig = 1 th_solution = fsolve(func, th_ig)[0] # update previous values self.__thnpp = self.__thnp self.__thnp = th_solution return th_solution
def vx_pert(mode, x, z, t, W, K, R1): vx_hat_vals = np.zeros((len(x), len(z)), dtype=complex) vx_vals = np.zeros((len(x), len(z)), dtype=complex) for i in range(len(x)): for j in range(len(z)): def func(r): return [r[0] - x[i] + xix(mode, r[0], r[1], t, W, K, R1), r[1] - z[j] + xiz(mode, r[0], r[1], t, W, K, R1)] sol = np.real(fsolve(func, [x[i],z[j]], xtol=1e-03)) if abs(sol[0]) <= K: vx_hat_vals[i,j] = (constB(mode, W, K, R1)*sc.cosh(m0(W)*sol[0]) + constC(mode, W, K, R1)*sc.sinh(m0(W)*sol[0])) elif sol[0] < -K: vx_hat_vals[i,j] = constA(mode, W, K, R1)*(sc.cosh(m1(W, R1)*sol[0]) + sc.sinh(m1(W, R1)*sol[0])) elif sol[0] > K: vx_hat_vals[i,j] = constD(mode, W, K, R1)*(sc.cosh(m2(W)*sol[0]) - sc.sinh(m2(W)*sol[0])) vx_vals[i,j] = vx_hat_vals[i,j] * np.exp(1j*(z[j]-t)) return vx_vals
def vz_pert(mode, x, z, t, W, K, R1): vz_hat_vals = np.zeros((len(x), len(z)), dtype=complex) vz_vals = np.zeros((len(x), len(z)), dtype=complex) for i in range(len(x)): for j in range(len(z)): def func(r): return [r[0] - x[i] + xix(mode, r[0], r[1], t, W, K, R1), r[1] - z[j] + xiz(mode, r[0], r[1], t, W, K, R1)] sol = np.real(fsolve(func, [x[i],z[j]], xtol=1e-03)) if abs(sol[0]) <= K: vz_hat_vals[i,j] = (1j * c0**2 / (c0**2 - W**2)) * m0(W)*(constB(mode, W, K, R1)*sc.sinh(m0(W)*sol[0]) + constC(mode, W, K, R1)*sc.cosh(m0(W)*sol[0])) elif sol[0] < -K: vz_hat_vals[i,j] = (1j * c1(R1)**2 / (c1(R1)**2 - W**2)) * m1(W, R1)*constA(mode, W, K, R1)*(sc.sinh(m1(W, R1)*sol[0]) + sc.cosh(m1(W, R1)*sol[0])) elif sol[0] > K: vz_hat_vals[i,j] = (1j * c2**2 / (c2**2 - W**2)) * m2(W)*constD(mode, W, K, R1)*(sc.sinh(m2(W)*sol[0]) - sc.cosh(m2(W)*sol[0])) vz_vals[i,j] = vz_hat_vals[i,j] * np.exp(1j*(z[j]-t)) return vz_vals # Displacement amplitude
def implied_volatility(self): impvol = lambda x: self.BS(self.S, self.K, self.T, x, self.r, self.q) - self.opt_price iv = fsolve(impvol, SOLVER_STARTING_VALUE, fprime=self._fprime, xtol=IMPLIED_VOLATILITY_TOLERANCE) return iv[0]
def euler_implicit(fun_yprime, t, y0, h): y = [y0] for i in range(len(t)-1): to_be_solved = lambda yp: y[i]+h*fun_yprime(yp, t[i+1])-yp y.append(fsolve(to_be_solved, y[i])) return y
def original_seeds_non_int(moved_seeds, mins, maxes, n, mode, x, z, t, W, K, R1): # CANT GET THIS TO WORK """ Find original seed points Parameters ---------- moved_seeds: list of tuples moved seed points, e.g. [(x,y,z),(x,y,z)] disp_x, disp_y, disp_z: 3d array x, y, z displacement mins, maxes: list mins and maxes of x, y, z coords. e.g. [xmin, ymin, zmin] """ def rescale_n(unscaled, axis): return unscaled * (maxes[axis]-mins[axis]) / n[axis] + mins[axis] def rescale_maxmin(unscaled, axis, minus_min=False): if minus_min == True: return (unscaled - mins[axis]) * n[axis] / (maxes[axis]-mins[axis]) else: return unscaled * n[axis] / (maxes[axis]-mins[axis]) seeds = [] for seed in moved_seeds: seed = list(seed) def function(orig_seed): return [orig_seed[0] - seed[0] + H * rescale_maxmin(np.real(sf.xix(mode, rescale_n(orig_seed[0], axis=0), rescale_n(orig_seed[1], axis=1), t, W, K, R1)), axis=0), orig_seed[1] - seed[1], # + H * rescale_maxmin(np.real(sf.xiz(mode, rescale_n(orig_seed[0], axis=0), rescale_n(orig_seed[1], axis=1), t, W, K, R1)), axis=2), orig_seed[2] - seed[2] + H * rescale_maxmin(np.real(sf.xiy(mode, rescale_n(orig_seed[0], axis=0), rescale_n(orig_seed[1], axis=1), t, W, K)), axis=1)] original_seed = list(np.real(fsolve(function, seed, xtol=1e-03))) seeds.append(tuple(original_seed)) return seeds
def rho_pert(mode, x, z, t, W, K, R1): rho_hat_vals = np.zeros((len(x), len(z)), dtype=complex) rho_vals = np.zeros((len(x), len(z)), dtype=complex) for i in range(len(x)): for j in range(len(z)): def func(r): return [r[0] - x[i] + xix(mode, r[0], r[1], t, W, K, R1), r[1] - z[j] + xiz(mode, r[0], r[1], t, W, K, R1)] sol = np.real(fsolve(func, [x[i],z[j]], xtol=1e-03)) if abs(sol[0]) <= K: rho_hat_vals[i,j] = m0(W)*(constB(mode, W, K, R1)*sc.sinh(m0(W)*sol[0]) + constC(mode, W, K, R1)*sc.cosh(m0(W)*sol[0])) * lamb00(W) / (c0**2 * m00(W)) # if mode in slow_surf_mode_options + slow_body_1_mode_options + slow_body_2_mode_options: # rho_hat_vals[i,j] = -rho_hat_vals[i,j] elif sol[0] < -K: rho_hat_vals[i,j] = constA(mode, W, K, R1)*(sc.sinh(m1(W, R1)*sol[0]) + sc.cosh(m1(W, R1)*sol[0])) * lamb1(W, R1) / c1(R1)**2 elif sol[0] > K: rho_hat_vals[i,j] = constD(mode, W, K, R1)*(sc.sinh(m2(W)*sol[0]) - sc.cosh(m2(W)*sol[0])) * lamb2(W) / c2**2 rho_vals[i,j] = rho_hat_vals[i,j] * np.exp(1j*(z[j]-t)) return rho_vals ############################################################################## # Alfven mixed driver functions # Frequency of the wave is dependant on x
def _schedule(self, rho, v): if rho < 1e-100: epsilon = 0 else: fun = lambda epsilon: pow(epsilon, 2) + v * pow(epsilon, 3 / 2) - pow(rho, 2) epsilon = optimize.fsolve(fun, rho / 2) return (epsilon)
def MultiDimNewtRaph(f,x0,dx=1e-6,args=(),ytol=1e-5,w=1.0,JustOneStep=False): """ A Newton-Raphson solver where the Jacobian is always re-evaluated rather than re-using the information as in the fsolve method of scipy.optimize """ #Convert to numpy array, force type conversion to float x=np.array(x0,dtype=np.float) error=999 J=np.zeros((len(x),len(x))) #If a float is passed in for dx, convert to a numpy-like list the same shape #as x if isinstance(dx,int): dx=dx*np.ones_like(float(x)) elif isinstance(dx,float): dx=dx*np.ones_like(x) r0=array(f(x,*args)) while abs(error)>ytol: #Build the Jacobian matrix by columns for i in range(len(x)): epsilon=np.zeros_like(x) epsilon[i]=dx[i] J[:,i]=(array(f(x+epsilon,*args))-r0)/epsilon[i] v=np.dot(-inv(J),r0) x=x+w*v #Calculate the residual vector at the new step r0=f(x,*args) error = np.max(np.abs(r0)) #Just do one step and stop if JustOneStep==True: return x return x
def FF_Yang_Dou(Re, rough): """ This function approximate the friction factor for the given Reynold's number and the relative roughness arrays with the Yang Dou approximation (see Yang, S. Dou, G. (2010). Turbulent drag reduction with polymer additive in rough pipes. Journal of Fluid Mechanics, 642 279-294). The function is implicit and utilize a numerical root finder Arguments: ReNum (ndarray-float): Reynold's number rough (ndarray-float): 1/relative roughness (w/roughness length scale) Returns: ndarray-float : frction factor """ ff_args = (Re, rough) sol_vbyu = fsolve(FF_Yang_Dou_residual, 15., ff_args) # sol_vbyu = brentq(FF_Yang_Dou_residual, 18., 100., ff_args) ff_Yang_Dou = 2 / sol_vbyu ** 2 Rplus = Re / (2 * sol_vbyu) ff_Man_Strkl = 0.143 / 4 / rough ** (1 / 3) if Rplus >= 100 * rough: ff = ff_Man_Strkl else: ff = ff_Yang_Dou if rough < 32 and ff > ff_Man_Strkl: ff = ff_Man_Strkl return ff #-----------------------------------------------------------------------------------------------------------------------
def solve(self,t=0): solution = fsolve(self.r,self.getSolutionVec(),band=self._band)
def conshist(c1, *HHparams): # unpack HHparams wbar = HHparams[0] rbar = HHparams[1] # initialize vectors for consumption, labor and savings chist = np.zeros(S) nhist = np.zeros(S) bhist = np.zeros(S+1) # set intitial consumption value chist[0] = c1 for t in range(0, S): # set up lambda function for fsolve with n as the input f1 = lambda n: LLEuler(n, chist[t], *HHparams) # solve LLEuler for value of n given c nhist[t] = opt.fsolve(f1, .99) if printcheck: # check that LLEuler is close to zero and report check = f1(nhist[t]) print("nhist", t, ": ", nhist[t], " check-n: ", check) # solve for b given c and n bhist[t+1] = wbar*nhist[t] + (1.+rbar)*bhist[t] - chist[t] # if not the final period solve for next period's c from interpemporal # Euler equation if t< S-1: chist[t+1] = chist[t]*(bet*(1.+rbar))**(1/sig) return chist, nhist, bhist # find IT Euler error at death, given a value for c(1), wbar, rbar and model # parameters
def findc1(c1, *HHparams): # This function is for fsolve it takes initial consumption and returns # final IT Euler error. chist, nhist, bhist = conshist(c1, *HHparams) return chist[S-1] - chib**(-1.0/sig)*bhist[S] # find updated wbar and rbar from initial guesses
def updatebar(bar, UPparams): # unpack the bar vector wbar = bar[0] rbar = bar[1] # set parameters to pass to findc1 in fsolve HHparams = (wbar, rbar, sig, chin, chib, b, ups, bet, S, printcheck) # solve findc1 for value of c1 c1opt = opt.fsolve(findc1, .1, args=HHparams) # check that final savings is close to zero and report # check = findc1(c1opt, *HHparams) # print "check-c: ", check, "wbar & rbar: ", wbar, rbar # get the consumer's full history chist, nhist, bhist = conshist(c1opt, *HHparams) # sum savings and labor to get aggregate capital annd labor inputs Kbar = np.sum(bhist) if Kbar < .01: Kbar = .01 Lbar = np.sum(nhist) # solve for the implied values of wages and interest rates wbarnew = (1-alf)*(Kbar/Lbar)**alf rbarnew = alf*(Lbar/Kbar)**(1-alf) - delta # put into barnew array barnew = np.array([wbarnew, rbarnew]) return barnew # Main program follows # set model parameter values
def LinApp_FindSS(funcname,param,guessXY,Zbar,nx,ny): ''' Finds the steady state for a DSGE model numerically Parameters ----------- funcname: function the name of the function which generates a column vector from ny+nx dynamic equations. The ny equations to be linearized into the form below in the first ny rows. A X(t) + B X(t-1) + C Y(t) + D Z(t) = 0 The function must be written so as to evaluate to zero for all rows in the steady state. param: array, dtype=float a vector of parameter values to be passed to funcname. guessXY: array, dtype=float guess for the steady state values of X and Y Zbar: array, dtype=float nz vector of Z steady state values nx: number, dtype=int number of X variables ny: number, dtype=int number of Y variables Returns -------- XYbar: array, dtype=float 1-by-(nx+ny) vector of X and Y steady state values, with the X values in positions 1 - nx and the Y values in nx+1 - nx+ny. ''' f = lambda XYbar: steady(XYbar, Zbar, funcname, param, nx, ny) XYbar = opt.fsolve(f, x0=guessXY) return XYbar
def LinApp_FindSS(funcname, param, guessXY, Zbar, nx, ny): # ''' # Finds the steady state for a DSGE model numerically # # Parameters # ----------- # funcname: function # the name of the function which generates a column vector # from ny+nx dynamic equations. # The ny equations to be linearized into the form below in the first # ny rows. # A X(t) + B X(t-1) + C Y(t) + D Z(t) = 0 # The function must be written so as to evaluate to zero for all rows # in the steady state. # # param: array, dtype=float # a vector of parameter values to be passed to funcname. # # guessXY: array, dtype=float # guess for the steady state values of X and Y # # Zbar: array, dtype=float # nz vector of Z steady state values # # nx: number, dtype=int # number of X variables # # ny: number, dtype=int # number of Y variables # # Returns # -------- # XYbar: array, dtype=float # 1-by-(nx+ny) vector of X and Y steady state values, with the X # values in positions 1 - nx and the Y values in nx+1 - nx+ny. # ''' f = lambda XYbar: steady(XYbar, Zbar, funcname, param, nx, ny) XYbar = opt.fsolve(f, guessXY) return XYbar
def fit(self, data, p = None, r = None): if p is None or r is None: av = np.average(data) va = np.var(data) r = (av*av)/(va-av) p = (va-av)/(va) sm = np.sum(data)/len(data) x = optimize.fsolve(self.mleFun, np.array([p, r]), args=(data, sm)) self.p = x[0] self.r = x[1]
def Thermosyphon(Fluid,L_a,L_c,L_e,D_o,D_i,F_r,h_eo,k_wall,h_co,T_eo,T_co): # Step 1: Design parameter specification L_eff=L_a+(L_c+L_e)/2 # thermosyphon effective length A_eo=np.pi*L_e*D_o # evaporator external area A_ei=np.pi*L_e*D_i # evaporator internal area A_co=np.pi*L_c*D_o # condenser external area A_ci=np.pi*L_c*D_i # Condenser internal cross-sectional area A_axial=np.pi/4*(D_o**2-D_i**2) # tube annuli cross-sectional area # Step 2: Determine thermal resistance R_eo =1/(A_eo*h_eo) # Evaporator external thermal resistance R_we=np.log(D_o/D_i)/(2*np.pi*L_e*k_wall)# evaporator wall thermal resistance R_wc=np.log(D_o/D_i)/(2*np.pi*L_c*k_wall)# condenser wall thermal resistance R_co = 1/(A_co*h_co)# condenser external thermal resistance R_axial = L_eff/(k_wall*A_axial) # tube axial thermal resistance # set of equation def f(x): # Qdot_rad: x[0] # T_we_o: x[1] # T_we_i: x[2] # Tsat_e: x[3] # T_wc_i: x[4] # T_wc_o: x[5] print('.') Psat_e = PropsSI('P','T',x[3],'Q',1,Fluid) DeltaP_v = ThermosyphonVapourPressureDrop(Fluid,x[3],x[0],D_i,L_eff) Psat_c = Psat_e-DeltaP_v Tsat_c = PropsSI('T','P',Psat_c,'Q',1,Fluid) R_ei = BoilingThermalResistanceESDU(Fluid,x[0],D_i,L_e,x[3],F_r) R_ci = CondThermalResistanceESDU(Fluid,x[0],Tsat_c,D_i,L_c) T_we = (x[1]+x[2])/2 T_wc = (x[4]+x[5])/2 Qdot_axial = (T_we-T_wc)/R_axial return [ x[0]+Qdot_axial-(T_eo-x[1])/R_eo, x[0]-(x[1]-x[2])/R_we, x[0]-(x[2]-x[3])/R_ei, x[0]-(Tsat_c-x[4])/R_ci, x[0]-(x[4]-x[5])/R_wc, x[0]+Qdot_axial-(x[5]-T_co)/R_co ] [Qdot_rad,T_we_o,T_we_i,Tsat_e,T_wc_i,T_wc_o]=fsolve(f,[1.0e3,500.0+273,250+273.0,180.0+273.15,170.0+273.15,150.0+273.15],ytol=1e-6) return Qdot_rad, T_we_o,T_we_i,Tsat_e,T_wc_i,T_wc_o
def DR_Solve(params,start): init_params(params) solt=root(DR_RF,x0=start,method='lm',tol=1e-12) if solt.get('success')==False: x, infodict, ier, mesg = fsolve(DR_RF, x0=start,maxfev=0, full_output=1,xtol=1e-6, epsfcn=1e-6) nfev=infodict['nfev'] else: x = solt.get('x') ier=1 mesg=solt.get('message') nfev=solt.get('nfev') if ier: print('Converged after', nfev, 'iterations',file=outfile) print('\n RESULT:',(x[0],x[1]),'\n',file=outfile) residual=DR(x[0]+1j*x[1],beta,tau,Tpar_Tperp,kperp,kpar,gam,eta,nb,theta,k) print('DR residual:', residual,file=outfile) if abs(residual)>1e-8: print('WARNING: Large residual, root may not be properly converged!',file=outfile) return None, None, None, None, 0,residual #convergence_issue+=1 omega_z = x[0]+1j*x[1] compeigvec=None #for complex components u, v = eig(mat) eps=2*amin(abs(u)) move=abs(x[0]-start[0])+abs(x[1]-start[1]) #will be 0 if it doesn't move # # Get the eigenvectors. # print('Smallest eigenvalue',eps,file=outfile) By=None dn=None for i in range(3): if abs(u[i]) < eps and eps < 1e-4: #print ("Smallest eigenvalue and corresponding eigenvector:",file=outfile) #print (abs(u[i]),file=outfile) compeigvec=np.array([0+1j*0,0+1j*0,0+1j*0]) for j in range(3): compeigvec[j]=v[j,i] By=compeigvec[0]/compeigvec[1] dn=compeigvec[2]/compeigvec[1] return x[0],x[1],By,dn,1,residual else: return None, None, None, None, 0,-1
def compute_rad_robin_eigenfrequencies(param, l, n_roots=10, show_plot=False): a2, a1, a0, alpha, beta = param eta = -a1 / 2. / a2 def characteristic_equation(om): if np.round(om, 200) != 0.: zero = (alpha + beta) * np.cos(om * l) + ((eta + beta) * (alpha - eta) / om - om) * np.sin(om * l) else: zero = (alpha + beta) * np.cos(om * l) + (eta + beta) * (alpha - eta) * l - om * np.sin(om * l) return zero def complex_characteristic_equation(om): if np.round(om, 200) != 0.: zero = (alpha + beta) * np.cosh(om * l) + ((eta + beta) * (alpha - eta) / om + om) * np.sinh(om * l) else: zero = (alpha + beta) * np.cosh(om * l) + (eta + beta) * (alpha - eta) * l + om * np.sinh(om * l) return zero # assume 1 root per pi/l (safety factor = 3) om_end = 3 * n_roots * np.pi / l start_values = np.arange(0, om_end, .1) om = ut.find_roots(characteristic_equation, 2 * n_roots, start_values, rtol=int(np.log10(l) - 6), show_plot=show_plot).tolist() # delete all around om = 0 om.reverse() for i in range(np.sum(np.array(om) < np.pi / l / 2e1)): om.pop() om.reverse() # if om = 0 is a root then add 0 to the list zero_limit = alpha + beta + (eta + beta) * (alpha - eta) * l if np.round(zero_limit, 6 + int(np.log10(l))) == 0.: om.insert(0, 0.) # regard complex roots om_squared = np.power(om, 2).tolist() complex_root = fsolve(complex_characteristic_equation, om_end) if np.round(complex_root, 6 + int(np.log10(l))) != 0.: om_squared.insert(0, -complex_root[0] ** 2) # basically complex eigenfrequencies om = np.sqrt(np.array(om_squared).astype(complex)) if len(om) < n_roots: raise ValueError("RadRobinEigenvalues.compute_eigen_frequencies()" "can not find enough roots") eig_frequencies = om[:n_roots] eig_values = a0 - a2 * eig_frequencies ** 2 - a1 ** 2 / 4. / a2 return eig_frequencies, eig_values
def arctan_fit(first_point, coef1, coef2, coef3, abil_deprec, init_guesses): ''' -------------------------------------------------------------------- This function fits an arctan function to the last 20 years of the ability levels of a particular ability group to extrapolate abilities by trying to match the slope in the 80th year and the ability depreciation rate between years 80 and 100. -------------------------------------------------------------------- INPUTS: first_point = scalar > 0, ability level at age 80 coef1 = scalar, coefficient in log ability equation on linear term in age coef2 = scalar, coefficient in log ability equation on quadratic term in age coef3 = scalar, coefficient in log ability equation on cubic term in age abil_deprec = scalar in (0, 1), ability depreciation rate between ages 80 and 100 init_guesses = (3,) vector, initial guesses OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: arc_error() arctan_func() OBJECTS CREATED WITHIN FUNCTION: a = scalar, scale parameter for arctan function b = scalar, curvature parameter for arctan function c = scalar, shift parameter for arctan function old_ages = (20,) vector, annual ages 81 to 100 abil_last = (20,) vector, extrapolated ability levels for ages 81 to 100 RETURNS: abil_last -------------------------------------------------------------------- ''' params = [first_point, coef1, coef2, coef3, abil_deprec] a, b, c = opt.fsolve(arc_error, init_guesses, params) old_ages = np.linspace(81, 100, 20) abil_last = arctan_func(old_ages, a, b, c) return abil_last
