def FF_Yang_Dou_residual(vbyu, *args):
"""
The Yang_Dou residual function; to be used by numerical root finder
"""
(Re, rough) = args
Rstar = Re / (2 * vbyu * rough)
theta = np.pi * np.log( Rstar / 1.25) / np.log(100 / 1.25)
alpha = (1 - np.cos(theta)) / 2
beta = 1 - (1 - 0.107) * (alpha + theta/np.pi) / 2
R = Re / (2 * vbyu)
rt = 1.
for i in range(1,5):
rt = rt - 1. / np.e * ( i / factorial(i) * (67.8 / R) ** (2 * i))
return vbyu - (1 - rt) * R / 4. - rt * (2.5 * np.log(R) - 66.69 * R**-0.72 + 1.8 - (2.5 * np.log(
(1 + alpha * Rstar / 5) / (1 + alpha * beta * Rstar / 5)) + (5.8 + 1.25) * (alpha * Rstar / (
5 + alpha * Rstar)) ** 2 + 2.5 * (alpha * Rstar / (5 + alpha * Rstar)) - (5.8 + 1.25)
* (alpha * beta * Rstar / (5 + alpha * beta * Rstar)) ** 2 - 2.5 * (alpha * beta * Rstar / (
5 + alpha * beta * Rstar))))
python类factorial()的实例源码
def approx_expm(self,M,exp_t, scaling_terms):
#approximate the exp at the beginning to estimate the number of taylor terms and scaling and squaring needed
U=np.identity(len(M),dtype=M.dtype)
Mt=np.identity(len(M),dtype=M.dtype)
factorial=1.0 #for factorials
for ii in xrange(1,exp_t):
factorial*=ii
Mt=np.dot(Mt,M)
U+=Mt/((2.**float(ii*scaling_terms))*factorial) #scaling by 2**scaling_terms
for ii in xrange(scaling_terms):
U=np.dot(U,U) #squaring scaling times
return U
def approx_exp(self,M,exp_t, scaling_terms):
# the scaling and squaring of matrix exponential with taylor expansions
U=1.0
Mt=1.0
factorial=1.0 #for factorials
for ii in xrange(1,exp_t):
factorial*=ii
Mt=M*Mt
U+=Mt/((2.**float(ii*scaling_terms))*factorial) #scaling by 2**scaling_terms
for ii in xrange(scaling_terms):
U=np.dot(U,U) #squaring scaling times
return U
def sph_harm(m, n, az, el, type='complex'):
'''Compute sphercial harmonics
Parameters
----------
m : (int)
Order of the spherical harmonic. abs(m) <= n
n : (int)
Degree of the harmonic, sometimes called l. n >= 0
az: (float)
Azimuthal (longitudinal) coordinate [0, 2pi], also called Theta.
el : (float)
Elevation (colatitudinal) coordinate [0, pi], also called Phi.
Returns
-------
y_mn : (complex float)
Complex spherical harmonic of order m and degree n,
sampled at theta = az, phi = el
'''
if type == 'legacy':
return scy.sph_harm(m, n, az, el)
elif type == 'real':
Lnm = scy.lpmv(_np.abs(m), n, _np.cos(el))
factor_1 = (2 * n + 1) / (4 * _np.pi)
factor_2 = scy.factorial(n - _np.abs(m)) / scy.factorial(n + abs(m))
if m != 0:
factor_1 = 2 * factor_1
if m < 0:
return (-1) ** m * _np.sqrt(factor_1 * factor_2) * Lnm * _np.sin(m * az)
else:
return (-1) ** m * _np.sqrt(factor_1 * factor_2) * Lnm * _np.cos(m * az)
else:
# For the correct Condon–Shortley phase, all m>0 need to be increased by 1
return (-1) ** _np.float_(m - (m < 0) * (m % 2)) * scy.sph_harm(m, n, az, el)
def test_factorial():
if scipy_factorial is None:
return
for i in range(20):
assert_equal(factorial(i), scipy_factorial(i))
def factorial(N):
"""Compute the factorial of N.
If N <= 16, use a fast lookup table; otherwise use scipy.special.factorial
"""
if N < len(FACTORIALS):
return FACTORIALS[N]
elif scipy_special is None:
raise ValueError("need scipy for computing larger factorials")
else:
return int(scipy_special.factorial(N))
def fun(self, x, *args):
self.nfev += 1
return (prod(x) - factorial(self.N)) ** 2.0
def weights(dim, degree):
# 1D sigma-points (x) and weights (w)
x, w = hermegauss(degree)
# hermegauss() provides weights that cause posdef errors
w = factorial(degree) / (degree ** 2 * hermeval(x, [0] * (degree - 1) + [1]) ** 2)
return np.prod(cartesian([w] * dim), axis=1)
ReactionSystem.py 文件源码
项目:Chemical-Reaction-System-Simulator
作者: axlevisu
项目源码
文件源码
阅读 20
收藏 0
点赞 0
评论 0
def run(self, t,seed=None,plot=False):
"""
Gillespie Implementation
"""
l = self.reactions[:,0]
r = self.reactions[:,1]
change=r-l
time =0
times =[]
output =[]
times.append(time)
output.append(self.population)
while time<t:
lam = np.prod(comb(self.population,l)*factorial(l),axis=-1)*self.rates
lam_sum = np.sum(lam)
dt = np.log(1.0/np.random.uniform(0,1))*(1.0/(lam_sum))
react = np.where(np.random.multinomial(1,lam/lam_sum)==1)[0][0]
self.population = self.population + change[react]
time += dt
output.append(self.population)
times.append(time)
output = np.array(output)
times = np.array(times)
if plot:
for i in xrange(output.shape[1]):
label = self.species[i]
plt.plot(times, output[:, i], label=label)
plt.legend(loc='best')
plt.xlabel('t')
plt.title('Population Profile')
plt.grid()
plt.show()
return times,output
# def main():
# reactions = [[[1,0],[0,1]], [[0,1],[1,0]]]
# rates = [1,1]
# system = MassActionSystem(reactions,rates)
# system.set_concentrations([2.,4.])
# system.run()
# print system.current_concentrations()
# print system.dydt()
# system = StochasticSystem(reactions,rates)
# population = [10,0]
# system.set_population(population)
# t,o = system.run(10,plot=True)
# return