我们从Python开源项目中，提取了以下45个代码示例，用于说明如何使用scipy.optimize.root()。
def find_first_best(self): ''' Find the first best allocation ''' model = self.model S, Theta, G = self.S, self.Theta, self.G Uc, Un = model.Uc, model.Un def res(z): c = z[:S] n = z[S:] return np.hstack([Theta * Uc(c, n) + Un(c, n), Theta * n - c - G]) res = root(res, 0.5 * np.ones(2 * S)) if not res.success: raise Exception('Could not find first best') self.cFB = res.x[:S] self.nFB = res.x[S:] # Multiplier on the resource constraint self.XiFB = Uc(self.cFB, self.nFB) self.zFB = np.hstack([self.cFB, self.nFB, self.XiFB])
def time0_allocation(self, B_, s_0): ''' Finds the optimal allocation given initial government debt B_ and state s_0 ''' model, pi, Theta, G, beta = self.model, self.pi, self.Theta, self.G, self.beta Uc, Ucc, Un, Unn = model.Uc, model.Ucc, model.Un, model.Unn # First order conditions of planner's problem def FOC(z): mu, c, n, Xi = z xprime = self.time1_allocation(mu)[2] return np.hstack([Uc(c, n) * (c - B_) + Un(c, n) * n + beta * pi[s_0] @ xprime, Uc(c, n) - mu * (Ucc(c, n) * (c - B_) + Uc(c, n)) - Xi, Un(c, n) - mu * (Unn(c, n) * n + Un(c, n)) + Theta[s_0] * Xi, (Theta * n - c - G)[s_0]]) # Find root res = root(FOC, np.array( [0, self.cFB[s_0], self.nFB[s_0], self.XiFB[s_0]])) if not res.success: raise Exception('Could not find time 0 LS allocation.') return res.x
def find_first_best(self): ''' Find the first best allocation ''' Para = self.Para S,Theta,Uc,Un,G = self.S,self.Theta,Para.Uc,Para.Un,self.G def res(z): c = z[:S] n = z[S:] return np.hstack( [Theta*Uc(c,n)+Un(c,n), Theta*n - c - G] ) res = root(res,0.5*np.ones(2*S)) if not res.success: raise Exception('Could not find first best') self.cFB = res.x[:S] self.nFB = res.x[S:] self.XiFB = Uc(self.cFB,self.nFB) #multiplier on the resource constraint. self.zFB = np.hstack([self.cFB,self.nFB,self.XiFB])
def find_first_best(self): ''' Find the first best allocation ''' model = self.model S, Theta, Uc, Un, G = self.S, self.Theta, model.Uc, model.Un, self.G def res(z): c = z[:S] n = z[S:] return np.hstack( [Theta * Uc(c, n) + Un(c, n), Theta * n - c - G] ) res = root(res, 0.5 * np.ones(2 * S)) if not res.success: raise Exception('Could not find first best') self.cFB = res.x[:S] self.nFB = res.x[S:] # multiplier on the resource constraint. self.XiFB = Uc(self.cFB, self.nFB) self.zFB = np.hstack([self.cFB, self.nFB, self.XiFB])
def time0_allocation(self,B_,s_0): ''' Finds the optimal allocation given initial government debt B_ and state s_0 ''' Para,Pi,Theta,G,beta = self.Para,self.Pi,self.Theta,self.G,self.beta Uc,Ucc,Un,Unn = Para.Uc,Para.Ucc,Para.Un,Para.Unn #first order conditions of planner's problem def FOC(z): mu,c,n,Xi = z xprime = self.time1_allocation(mu)[2] return np.hstack([ Uc(c,n)*(c-B_) + Un(c,n)*n + beta*Pi[s_0].dot(xprime), Uc(c,n) - mu*(Ucc(c,n)*(c-B_) + Uc(c,n)) - Xi, Un(c,n) - mu*(Unn(c,n)*n+Un(c,n)) + Theta[s_0]*Xi, (Theta*n - c - G)[s_0] ]) #find root res = root(FOC,np.array([0.,self.cFB[s_0],self.nFB[s_0],self.XiFB[s_0]])) if not res.success: raise Exception('Could not find time 0 LS allocation.') return res.x
def _nodes_LGL_old(self, n): """Return Legendre-Gauss-Lobatto nodes. Gauss-Lobatto nodes are roots of P'_{n-1}(x) and -1, 1. ref. http://keisan.casio.jp/exec/system/1360718708 """ x0 = np.zeros(0) for i in range(2, n): xi = (1-3.0*(n-2)/8.0/(n-1)**3)*np.cos((4.0*i-3)/(4.0*(n-1)+1)*np.pi) x0 = np.append(x0, xi) x0 = np.sort(x0) roots = np.zeros(0) for x in x0: optResult = optimize.root(self._LegendreDerivative, x, args=(n-1,)) roots = np.append(roots, optResult.x) nodes = np.hstack((-1, roots, 1)) return nodes
def stat_estimator(samples, cutoff=200, confidence=0.99): '''Max Likelihood Estimator for censored exponential distribution. See "Estimation of Parameters of Truncated or Censored Exponential Distributions", Walter L. Deemer and David F. Votaw''' samples = samples.astype(float) n = (samples<cutoff).sum() N = len(samples) estimate = n/samples.sum() y_conf = stats.norm.ppf((1+confidence)/2) y = lambda c: N**0.5*(estimate/c-1)*(1-np.exp(-c*cutoff))**0.5 low = optimize.root(lambda c: y(c)-y_conf, estimate) high = optimize.root(lambda c: y(c)+y_conf, estimate) if not (low.success and high.success): raise RuntimeError('Could not find confidence interval for the given samples!') return np.array([1/estimate, 1/high.x, 1/low.x])
def inv_cdf_dlf(p, A, m1, a1, m2, a2, start=-26, end=-15): ''' Inverse Cumulative Schechter function. Second LF is set to be 2*A of first LF. @param p: probability @param A: Scale factor @param m1: Knee of distribution 1 @param a1: Faint-end turnover of first lf @param m2: Knee of distribution 2 @param a2: Faint-end turnover of second lf @param start: Brightest magnitude @param end: Faintest possible magnitude @return Magnitude associated with cdf probability p ''' def get_root(p): return root(lambda x: cdf_dlf(x,A,m1,a1,m2,a2,start)-p, (start + end)/2).x[0] if np.isscalar(p): return get_root(p) else: return np.fromiter(map(get_root,p),np.float,count=len(p))
def get_Mistar(self): """ get (Mi*)_max that is maximum inlet velocity ratio for supersonic starting""" Mistar0 = 1.5 def func(Mistar): a = self.get_Q()/(1-self.get_C()) b = Mistar**(2*self.gamma/(self.gamma -1)) * ((1-((self.gamma - 1)/(self.gamma + 1))*Mistar**2)/(Mistar**2 - ((self.gamma -1)/(self.gamma + 1))))**(1/(self.gamma -1)) return b - a sol = optimize.root(func,Mistar0) Mistar_max = sol.x[0] return Mistar_max
def get_all_stabilities(self, l0=np.array([1.0, 1.0]), isFullOutput=False): funcs = self.get_stability_functions() l = [] for f in funcs: l_full = optimize.root(f, l0, tol=1e-14, method='lm') l.append(l_full) if isFullOutput: return l else: l_num = [] for el in l: l_tmp = el.x[0] + 1j * el.x[1] l_num.append(l_tmp) l_num = np.array(l_num) return l_num
def solve_matrix_free(self): """Finds the stationary state using matrix free methods like broyden, krylov, etc.""" solmethod = self.funcp.solmethod # phi0_init = self.funcp.phi0_init if phi0_init is None: self.funcp.print_warning(0, "WARNING: For mfreeq=True no phi0_init is specified. "+ "Using phi0_init[0]=1.0 as a default. "+ "This warning will not be shown again.") phi0_init = self.set_phi0_init() # solmethod = solmethod if solmethod != 'n' else 'krylov' try: self.sol0 = optimize.root(self.generate_vec, phi0_init, args=(self), method=solmethod) self.phi0 = self.sol0.x self.success = self.sol0.success except Exception as exept: self.funcp.print_error(exept) self.phi0 = 0*phi0_init self.success = False self.funcp.solmethod = solmethod
def updateLocation(self,r0): lb = np.r_[-6.,-6.,-5.] ub = np.r_[6.,6.,-0.1] # Sol = op.minimize(self.computeMisfit,r0,method='Powell',options={'xtol':1e-5}) Sol = op.root(self.computeVecFcn,r0,method='lm',options={'xtol':1e-5}) r0 = Sol.x return r0,Sol ############################################################################## # DEFINE TEMTADSproblem class ##############################################################################
def updateLocation(self,r0): # Sol = op.minimize(self.computeMisfit,r0,method='dogleg') Sol = op.root(self.computeVecFcn,r0,method='lm',options={'xtol':1e-5}) r0 = Sol.x return r0,Sol ############################################################################## # DEFINE MPVproblem class ##############################################################################
def updateLocation(self,r0): # Sol = op.minimize(self.computeMisfit,r0,method='dogleg') Sol = op.root(self.computeVecFcn,r0,method='lm',options={'xtol':1e-5}) r0 = Sol.x return r0,Sol
def time1_allocation(self, mu): ''' Computes optimal allocation for time t\geq 1 for a given \mu ''' model = self.model S, Theta, G = self.S, self.Theta, self.G Uc, Ucc, Un, Unn = model.Uc, model.Ucc, model.Un, model.Unn def FOC(z): c = z[:S] n = z[S:2 * S] Xi = z[2 * S:] return np.hstack([Uc(c, n) - mu * (Ucc(c, n) * c + Uc(c, n)) - Xi, # FOC of c Un(c, n) - mu * (Unn(c, n) * n + Un(c, n)) + \ Theta * Xi, # FOC of n Theta * n - c - G]) # Find the root of the first order condition res = root(FOC, self.zFB) if not res.success: raise Exception('Could not find LS allocation.') z = res.x c, n, Xi = z[:S], z[S:2 * S], z[2 * S:] # Compute x I = Uc(c, n) * c + Un(c, n) * n x = np.linalg.solve(np.eye(S) - self.beta * self.pi, I) return c, n, x, Xi
def time1_allocation(self,mu): ''' Computes optimal allocation for time t\geq 1 for a given \mu ''' Para = self.Para S,Theta,G,Uc,Ucc,Un,Unn = self.S,self.Theta,self.G,Para.Uc,Para.Ucc,Para.Un,Para.Unn def FOC(z): c = z[:S] n = z[S:2*S] Xi = z[2*S:] return np.hstack([ Uc(c,n) - mu*(Ucc(c,n)*c+Uc(c,n)) -Xi, #foc c Un(c,n) - mu*(Unn(c,n)*n+Un(c,n)) + Theta*Xi, #foc n Theta*n - c - G #resource constraint ]) #find the root of the FOC res = root(FOC,self.zFB) if not res.success: raise Exception('Could not find LS allocation.') z = res.x c,n,Xi = z[:S],z[S:2*S],z[2*S:] #now compute x I = Uc(c,n)*c + Un(c,n)*n x = np.linalg.solve(np.eye(S) - self.beta*self.Pi, I ) return c,n,x,Xi
def find_first_best(self): ''' Find the first best allocation ''' Para = self.Para S,Theta,Uc,Un,G = self.S,self.Theta,Para.Uc,Para.Un,self.G def res(z): c = z[:S] n = z[S:] return np.hstack( [Theta*Uc(c,n)+Un(c,n), Theta*n - c - G] ) res = root(res,0.5*np.ones(2*S)) if not res.success: raise Exception('Could not find first best') self.cFB = res.x[:S] self.nFB = res.x[S:] IFB = Uc(self.cFB,self.nFB)*self.cFB + Un(self.cFB,self.nFB)*self.nFB self.xFB = np.linalg.solve(np.eye(S) - self.beta*self.Pi, IFB) self.zFB = {} for s in range(S): self.zFB[s] = np.hstack([self.cFB[s],self.nFB[s],self.xFB])
def time1_allocation(self, mu): ''' Computes optimal allocation for time t\geq 1 for a given \mu ''' model = self.model S, Theta, G, Uc, Ucc, Un, Unn = self.S, self.Theta, self.G, model.Uc, model.Ucc, model.Un, model.Unn def FOC(z): c = z[:S] n = z[S:2 * S] Xi = z[2 * S:] return np.hstack([ Uc(c, n) - mu * (Ucc(c, n) * c + Uc(c, n)) - Xi, # foc c Un(c, n) - mu * (Unn(c, n) * n + Un(c, n)) + Theta * Xi, # foc n Theta * n - c - G # resource constraint ]) # find the root of the FOC res = root(FOC, self.zFB) if not res.success: raise Exception('Could not find LS allocation.') z = res.x c, n, Xi = z[:S], z[S:2 * S], z[2 * S:] # now compute x I = Uc(c, n) * c + Un(c, n) * n x = np.linalg.solve(np.eye(S) - self.beta * self.pi, I) return c, n, x, Xi
def time0_allocation(self, B_, s_0): ''' Finds the optimal allocation given initial government debt B_ and state s_0 ''' model, pi, Theta, G, beta = self.model, self.pi, self.Theta, self.G, self.beta Uc, Ucc, Un, Unn = model.Uc, model.Ucc, model.Un, model.Unn # first order conditions of planner's problem def FOC(z): mu, c, n, Xi = z xprime = self.time1_allocation(mu)[2] return np.hstack([ Uc(c, n) * (c - B_) + Un(c, n) * n + beta * pi[s_0].dot(xprime), Uc(c, n) - mu * (Ucc(c, n) * (c - B_) + Uc(c, n)) - Xi, Un(c, n) - mu * (Unn(c, n) * n + Un(c, n)) + Theta[s_0] * Xi, (Theta * n - c - G)[s_0] ]) # find root res = root(FOC, np.array( [0., self.cFB[s_0], self.nFB[s_0], self.XiFB[s_0]])) if not res.success: raise Exception('Could not find time 0 LS allocation.') return res.x
def find_first_best(self): ''' Find the first best allocation ''' Para = self.Para S,Theta,Uc,Un,G = self.S,self.Theta,Para.Uc,Para.Un,self.G def res(z): c = z[:S] n = z[S:] return np.hstack( [Theta*Uc(c,n)+Un(c,n), Theta*n - c - G] ) res = root(res,0.5*np.ones(2*S)) if not res.success: raise Exception('Could not find first best') self.cFB = res.x[:S] self.nFB = res.x[S:] IFB = Uc(self.cFB,self.nFB)*self.cFB + Un(self.cFB,self.nFB)*self.nFB self.xFB = np.linalg.solve(np.eye(S) - self.beta*self.Pi, IFB) self.zFB = {} for s in range(S): self.zFB[s] = np.hstack([self.cFB[s],self.nFB[s],self.Pi[s].dot(self.xFB),0.])
def make_Ive_array(kperp,nb,Tpar_Tperp): global ive_arr ive_arr = np.empty((2*nb+2),dtype=np.float64) b=0.5*kperp**2/Tpar_Tperp js=np.arange(0,2*nb+2) j_shift=js-nb-1 for j in js: ive_arr[j] = ive(j_shift[j],b) #Generate plasma dispersion function array -- this has to be done any time omega or kpar #changes, i.e. also when the root finder moves
def DR(*arg): return dispersion_relation_analytical(*arg) #dispersion relation evaluation used by the root finding routine
def DR_RF(x): #For very negative imaginary parts, sometimes NaNs in the plasma dispersion function #occur. We just return 0,0 for gamma<-20, stopping the root solver. If all goes well, #the outer routines should discover the discontinuity of that solution and try again #at a different point, hopefully circumventing this problem. if x[1]<-20: return [1e9,1e9] omega=x[0]+1j*x[1] out=DR(omega,beta,tau,Tpar_Tperp,kperp,kpar,gam,eta,nb,theta,k) return [out.real,out.imag]
def DR_point(params,x): init_params(params,True) omega=x[0]+1j*x[1] return abs(DR(omega,beta,tau,Tpar_Tperp,kperp,kpar,gam,eta,nb,theta,k)) #Runs the root finder from a given start value and returns the result
def inverse_nfw_cumulative(self, p): ''' inverse of radial nfw cumulative distribution @param p: Probability @return float: Radius corresponding to probability p ''' R0 = self.__R0 def get_root(p): return root(lambda z: self.nfw_cumulative(z) - p,R0/2).x[0] if np.isscalar(p): get_root(p) else: return np.fromiter(map(get_root,p),np.float, count=len(p))
def get_mach_from_prandtle_meyer(self, v1): mach0 = 1.0 def func(mach, v1): return self.get_Pr(mach) - v1 sol = optimize.root(func, mach0, args=(v1)) mach = sol.x[0] return mach
def get_n_s(cvec, w, sigma, chi_n_vec, l_tilde, b, upsilon): n_args = (cvec, w, sigma, chi_n_vec, l_tilde, b, upsilon) n_guess = 0.5 * l_tilde * np.ones_like(cvec) results_n = opt.root(get_n_errors, n_guess, args=(n_args), method='lm') nvec = results_n.x return nvec
def get_n_s(cvec, w, sigma, chi_n_vec, l_tilde, b, upsilon, evec): n_args = (cvec, w, sigma, chi_n_vec, l_tilde, b, upsilon, evec) n_guess = 0.5 * l_tilde * np.ones_like(cvec) results_n = opt.root(get_n_errors, n_guess, args=(n_args), method='lm') nvec = results_n.x return nvec
def get_n_s(n_args): (cvec, bvec, r, w, mtrxparams, factor, sigma, chi_n_vec, l_tilde, b_ellip, upsilon) = n_args n_guess = 0.5 * l_tilde * np.ones_like(cvec) results_n = opt.root(get_n_errors, n_guess, args=(n_args), method='lm') nvec = results_n.x return nvec
def get_cnbvecs(c1, args): ''' bvec = (S,) vector, b_2, b3,...b_{S+1} ''' (r, w, X, factor, beta, sigma, chi_n_vec, l_tilde, b_ellip, upsilon, S, etrparams, mtrxparams, mtryparams) = args cvec = np.zeros(S) cvec[0] = c1 nvec = np.zeros(S) bvec = np.zeros(S) n1_args = (c1, 0.0, r, w, mtrxparams, factor, sigma, chi_n_vec[0], l_tilde, b_ellip, upsilon) n1 = get_n_s(n1_args) nvec[0] = n1 bs_args = (r, w, X, factor, etrparams) b2 = get_savings(c1, n1, 0.0, bs_args) bvec[0] = b2 for s in range(1, S): cs_init = cvec[s - 1] ns_init = nvec[s - 1] cn_init = np.array([cs_init, ns_init]) c_sm1 = cvec[s - 1] b_s = bvec[s - 1] chi_n_s = chi_n_vec[s] cn_args = (c_sm1, b_s, r, w, factor, beta, sigma, chi_n_s, l_tilde, b_ellip, upsilon, mtrxparams, mtryparams) results_cn = opt.root(get_cn, cn_init, args=(cn_args)) c_s, n_s = results_cn.x cvec[s] = c_s nvec[s] = n_s b_s = bvec[s - 1] b_sp1 = get_savings(c_s, n_s, b_s, bs_args) bvec[s] = b_sp1 return cvec, nvec, bvec
def get_stability(self, l0=np.array([1.0, 1.0])): funcs = self.get_stability_functions() l = [] for f in funcs: l_full = optimize.root(f, l0, tol=1e-14) l_num = l_full.x[0] + 1j * l_full.x[1] l.append(l_num) l = np.array(l) i_max = np.argmax(np.real(l)) return np.real(l[i_max]), np.imag(l[i_max])
def get_stability2(w, k, h, wc, tau, omega, topology, twist_number): d = topology.get_couling_derivate_matrix(h, twist_number, omega * tau) em, vm = np.linalg.eig(d) b = 1.0 / wc d_sum = np.sum(d[0, :]) # assumes that the sum of each row is the same for all rows lambda_zeta = [] for i_eigen in range(len(em)): zeta = em[i_eigen] print('\n\nzeta:', zeta, '\n\n') def func(l): alpha = np.real(zeta) beta = np.imag(zeta) x = np.zeros(2) x[0] = l[0] + b * l[0]**2 - b * l[1]**2 + k * d_sum - k * alpha * np.exp(-l[0] * tau) * np.cos(l[1] * tau) - k * beta * np.exp(-l[0] * tau) * np.sin(l[1] * tau) x[1] = l[1] + 2 * b * l[0] * l[1] + k * alpha * np.exp(-l[0] * tau) * np.sin(l[1] * tau) - k * beta * np.exp(-l[0] * tau) * np.cos(l[1] * tau) return x l_opt = optimize.root(func, np.array([0.1, 0.1]), tol=1e-14) l_tmp = l_opt.x[0] + 1j * l_opt.x[1] # Ignore solution for the eigenvector (a, a, a, ...) if np.max(np.abs(np.diff(vm[:, i_eigen]))) >= 1e-9: lambda_zeta.append(l_tmp) lambda_zeta = np.array(lambda_zeta) return lambda_zeta[np.argmax(np.real(lambda_zeta))]
def get_stability(self, l0=np.array([1.0, 1.0])): funcs = self.get_stability_functions() l = [] for f in funcs: l_full = optimize.root(f, l0, tol=1e-14) l_num = l_full.x[0] + 1j * l_full.x[1] l.append(l_num) l = np.array(l) i_max = np.argmax(np.real(l)) return l[i_max]
def solve(self, circuit, initial_conditions, **kwargs): ndarray_initial_conditions = np.array(initial_conditions) self._solution = root(self._get_equations_error, ndarray_initial_conditions, args=circuit) return circuit
def get_reduced_volume(self, Tr, pr, fluid, **kwargs): """ Returns the reduced volume :math:`v_r`. This is done by solving the following nonlinear equation: .. math:: \\frac{p_r v_r}{T_r} = 1 + \\frac{B(T_r)}{v_r} + \\frac{C(T_r)}{v_r^2} + \\frac{D(T_r)}{v_r^5} + \\frac{c_4}{T_r^3 v_r^2} (\\beta + \\frac{\gamma}{v_r^2}) \exp(-\\frac{\gamma}{v_r^2}) :param Tr: reduced temperature (:math:`T_r = T/T_c`) :type Tr: float :param pr: reduced pressure (:math:`p_r = p/p_c`) :type pr: float :param fluid: specifies if the fluid is *simple* or *reference* :type fluid: str :keyword init_vol: (optional kwarg) initial guess for reduced volume to be used in the nonlinear solver :type init_vol: float :return: reduced volume :rtype: float """ if fluid=='reference': i = 1 elif fluid=='simple': i = 0 else: return None B = self.get_B(Tr, fluid) C = self.get_C(Tr, fluid) D = self.get_D(Tr, fluid) def objfun(x): return -pr*x/Tr + 1 + B/x + C/x**2 + D/x**5 + self.c4[i]/Tr**3/x**2 * (self.beta[i] + self.gamma[i]/x**2) * exp(-self.gamma[i]/x**2) # initial gues for nonlinear solver if 'init_vol' in kwargs: init_vol = kwargs['init_vol'] else: init_vol = Tr/pr #if pr<1: # init_vol = 10. #else: # init_vol = 0.1 #if Tr/pr > 1: # init_vol = 10. #else: # init_vol = Tr/pr #vol, info, ier, mesg = fsolve(objfun, init_vol) result = spo.root(objfun, init_vol, method='lm') #print(result) if len(result.x) == 1: return result.x[0] else: return result.x
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 transition_arc(self, Rstar0, vmin, vmax, angle, delta_v): """ lower transition arc function Args: Rstar0 (float) : initial R* vmin (float) : small v(Prandtle-Meyer angle) [deg] vmax (float) : large v(Prandtle-Meyer angle) [deg] angle (float) : rotate angle from Y* to y* [deg] delta_v (float) : interval Prandtle-Meyer angle [deg] """ kmin = 1 kmax = int((vmax - vmin)/delta_v) + 1 k = np.arange(kmin, kmax) xstar_l = np.zeros(k.size) ystar_l = np.zeros(k.size) xstar_l[-1] = 0 ystar_l[-1] = Rstar0 func_Rstar = 2 * np.radians(vmax) - np.pi / 2 * (np.sqrt((self.gamma + 1) / (self.gamma - 1)) - 1) - np.radians(2 * (k - 1) * delta_v) def func(Rstar, func_Rstar_i): return 2 * self.chara_line(Rstar) - func_Rstar_i Rstar = np.zeros(func_Rstar.size) for (i, f_Rstar) in enumerate(func_Rstar): sol = optimize.root(func, Rstar0, args=(f_Rstar)) Rstar[i] = sol.x[0] fai_k = np.radians(vmax - vmin - (k - 1) * delta_v) xstar_k = - Rstar * np.sin(fai_k) ystar_k = Rstar * np.cos(fai_k) myu_k = - np.arcsin(np.sqrt(((self.gamma + 1) / 2) * Rstar**2 - (self.gamma - 1) / 2)) m_k =[] for (i, j) in enumerate(k[:-1]): m_k += [np.tan((fai_k[i] + fai_k[j])/2 + (myu_k[i] + myu_k[j])/2)] mbar_k = np.tan(fai_k[1:]) Xstar_l, Ystar_l = self.rotate(xstar_l[-1], ystar_l[-1], angle) x, y = [Xstar_l], [Ystar_l] xprev = 0 for i in range(kmax - kmin - 2, 0, -1): a = ystar_l[i+1] - mbar_k[i] * xstar_l[i+1] b = ystar_k[i] - m_k[i] * xstar_k[i] c = m_k[i] - mbar_k[i] xstar_l[i] = (a - b) / c ystar_l[i] = (m_k[i]*a -mbar_k[i]*b)/c if (angle > 0): xtemp = xstar_l[i] assert xtemp < xprev, "initial vu or vl must be changed" else: xtemp = - xstar_l[i] assert xtemp > xprev, "initial vu or v0 must be changed" Xstar_l, Ystar_l = self.rotate(xtemp, ystar_l[i], angle) x += [(Xstar_l)] y += [(Ystar_l)] xprev = xtemp return x, y
def calc_leading_edge(self, delta=8, offset=0.1): """ if you don't want to make the end small angle, adjust the end by change the angle and the offset Args: delta (float) : Increment of edge angle offset (float) : offset ratio by G*(shift) from row lower-surface """ t = self.shift * offset x1 = self.upper_convex_in_x_end y1 = self.upper_convex_in_y_end x2 = self.lower_concave_in_x_end y2 = self.lower_concave_in_y_end def func1(x): func1 = np.tan(np.deg2rad(self.beta_in)) * x + y1 - np.tan(np.deg2rad(self.beta_in)) * x1 - t func2 = np.tan(np.deg2rad(self.beta_in + delta)) * x + y2 - np.tan(np.deg2rad(self.beta_in + delta)) * x2 + self.shift return func1 - func2 x = optimize.root(func1,-1) y = np.tan(np.deg2rad(self.beta_in + delta)) * x.x[0] + y2 - np.tan(np.deg2rad(self.beta_in + delta)) * x2 + self.shift self.upper_straight_in_x = [x1, x.x[0]] self.upper_straight_in_y = [y1, y] self.edge_straight_in_x = [x2,x.x[0]] self.edge_straight_in_y = [y2 + self.shift + t,y] x1 = self.upper_convex_out_x_end y1 = self.upper_convex_out_y_end x2 = self.lower_concave_out_x_end y2 = self.lower_concave_out_y_end def func2(x): func1 = np.tan(np.deg2rad(self.beta_out)) * x + y1 - np.tan(np.deg2rad(self.beta_out)) * x1 - t func2 = np.tan(np.deg2rad(self.beta_out - delta)) * x + y2 - np.tan(np.deg2rad(self.beta_out - delta)) * x2 + self.shift return func1 - func2 x = optimize.root(func2,1) y = np.tan(np.deg2rad(self.beta_out - delta)) * x.x[0] + y2 - np.tan(np.deg2rad(self.beta_out - delta)) * x2 + self.shift self.upper_straight_out_x = [x1, x.x[0]] self.upper_straight_out_y = [y1, y] self.edge_straight_out_x = [x2,x.x[0]] self.edge_straight_out_y = [y2 + self.shift + t,y] self.shift = self.shift + t
def get_stability(n, w, k, h, m, tau, omega, wc): # Dependent parameters b = 1.0 / wc dhdt = h.get_derivative() phi_m = (2 * np.pi * m) / n # Compute stability for root-based frequencies alpha_plus = k * dhdt(-omega * tau - phi_m) alpha_minus = k * dhdt(-omega * tau + phi_m) e_mat = np.zeros((n, n)) e_mat[0, -1] = alpha_plus e_mat[0, 1] = alpha_minus for ik in range(1, n - 1): e_mat[ik, ik - 1] = alpha_plus e_mat[ik, ik + 1] = alpha_minus e_mat[-1, 0] = alpha_minus e_mat[-1, -2] = alpha_plus em, vm = np.linalg.eig(e_mat) lambda_nu = [] for inu in range(len(em)): nu = em[inu] #print 'tau = %.3f: %.6f + 1j * %.6f' % (tau, np.real(nu), np.imag(nu)) def func(l): mu = np.real(nu) gamma = np.imag(nu) x = np.zeros(2) x[0] = b * l[0]**2 - b * l[1]**2 + l[0] + 0.5 * (alpha_plus + alpha_minus) - 0.5 * mu * np.exp(-l[0] * tau) * np.cos(l[1] * tau) - 0.5 * gamma * np.exp(-l[0] * tau) * np.sin(l[1] * tau) x[1] = 2 * b * l[0] * l[1] + l[1] + 0.5 * mu * np.exp(-l[0] * tau) * np.sin(l[1] * tau) - 0.5 * gamma * np.exp(-l[0] * tau) * np.cos(l[1] * tau) return x l_opt = optimize.root(func, np.array([1.0, 1.0]), tol=1e-15) l_tmp = l_opt.x[0] + 1j * l_opt.x[1] if np.max(np.abs(np.diff(vm[:, inu]))) >= 1e-9: lambda_nu.append(l_tmp) lambda_nu = np.array(lambda_nu) return lambda_nu[np.argmax(np.real(lambda_nu))] # ##############################################################################
def get_stability(n, w, k, h, m, tau, omega, wc): # Dependent parameters b = 1.0 / wc dhdt = h.get_derivative() phi_m = (2 * np.pi * m) / n # Compute stability for root-based frequencies alpha_plus = k * dhdt(-omega * tau - phi_m) alpha_minus = k * dhdt(-omega * tau + phi_m) e_mat = np.zeros((n, n)) e_mat[0, -1] = alpha_plus e_mat[0, 1] = alpha_minus for ik in range(1, n - 1): e_mat[ik, ik - 1] = alpha_plus e_mat[ik, ik + 1] = alpha_minus e_mat[-1, 0] = alpha_minus e_mat[-1, -2] = alpha_plus em, vm = np.linalg.eig(e_mat) lambda_nu = [] for inu in range(len(em)): nu = em[inu] #print 'tau = %.3f: %.6f + 1j * %.6f' % (tau, np.real(nu), np.imag(nu)) def func(l): mu = np.real(nu) gamma = np.imag(nu) x = np.zeros(2) x[0] = b * l[0]**2 - b * l[1]**2 + l[0] + 0.5 * (alpha_plus + alpha_minus) - 0.5 * mu * np.exp(-l[0] * tau) * np.cos(l[1] * tau) - 0.5 * gamma * np.exp(-l[0] * tau) * np.sin(l[1] * tau) x[1] = 2 * b * l[0] * l[1] + l[1] + 0.5 * mu * np.exp(-l[0] * tau) * np.sin(l[1] * tau) - 0.5 * gamma * np.exp(-l[0] * tau) * np.cos(l[1] * tau) return x l_opt = optimize.root(func, np.array([1.0, 1.0]), tol=1e-15) l_tmp = l_opt.x[0] + 1j * l_opt.x[1] if np.max(np.abs(np.diff(vm[:, inu]))) >= 1e-9: lambda_nu.append(l_tmp) lambda_nu = np.array(lambda_nu) return lambda_nu[np.argmax(np.real(lambda_nu))] # ############################################################################## # 1. Define triangle wave # ##############################################################################
