def smeft_evolve(C_in, scale_high, scale_in, scale_out, **kwargs):
def func(y, t0):
return beta.beta_array(C=beta.C_array2dict(y.view(complex)),
HIGHSCALE=scale_high).view(float)/(16*pi**2)
y0 = beta.C_dict2array(C_in).view(float)
sol = odeint(func=func, y0=y0, t=(log(scale_in), log(scale_out)), **kwargs)
return beta.C_array2dict(sol[1].view(complex))
python类odeint()的实例源码
def simulate(crn, par, initial_cond, start_t, end_t, incr):
"""Simulate the deterministic dynamics."""
# time
times = np.arange(start_t, end_t, incr)
# derivatives
eqs = map(lambda e: e.subs(par.items()), crn.equations())
# integration
sol = integrate.odeint(lambda x, t: odes(x, t, map(sp.Symbol, crn.species), eqs),
[initial_cond[s] for s in crn.species],
times)
return dict(zip(crn.species, np.transpose(sol))), times
def newton_ang_func(L_val,c2,m):
L = L_val
slope = (m*(m+1)+c2-L)/(m+1)/2.
z0 = [1+step*slope,slope]
z = odeint(f,z0,x_ang,args=(c2,L,m))
temp=1-pow(x_ang,2.0)
temp=pow(temp,m/2.)
zz=temp*z[:,0]
first_zz = np.array([1])
zz=np.append(first_zz, zz)
sloper = -(m*(m+1)+c2-L)/(m+1)/2.
z0r = [1-step*sloper,sloper]
zr = odeint(f,z0r,x_angr,args=(c2,L,m))
zzr=temp*zr[:,0]
zzr=np.append(first_zz, zzr)
return z[:,1][-1]
def Main(self):
"""
Main demo for the Hodgkin Huxley neuron model
"""
X = odeint(self.dALLdt, [-65, 0.05, 0.6, 0.32], self.t, args=(self,))
V = X[:,0]
m = X[:,1]
h = X[:,2]
n = X[:,3]
ina = self.I_Na(V, m, h)
ik = self.I_K(V, n)
il = self.I_L(V)
plt.figure()
plt.subplot(3,1,1)
plt.title('Hodgkin-Huxley Neuron')
plt.plot(self.t, V, 'k')
plt.ylabel('V (mV)')
plt.xticks([])
# plt.subplot(4,1,2)
# plt.plot(self.t, ina, 'c', label='$I_{Na}$')
# plt.plot(self.t, ik, 'y', label='$I_{K}$')
# plt.plot(self.t, il, 'm', label='$I_{L}$')
# plt.ylabel('Current')
# plt.xticks([])
# plt.legend(loc='upper center', ncol=3, prop=fontP)
plt.subplot(3,1,2)
plt.plot(self.t, m, 'r', label='m')
plt.plot(self.t, h, 'g', label='h')
plt.plot(self.t, n, 'b', label='n')
plt.ylabel('Gating Value')
plt.xticks([])
plt.legend(loc='upper center', ncol=3, prop=fontP)
plt.subplot(3,1,3)
i_inj_values = [self.I_inj(t) for t in self.t]
plt.plot(self.t, i_inj_values, 'k')
plt.xlabel('t (ms)')
plt.ylabel('$I_{inj}$ ($\\mu{A}/cm^2$)')
plt.ylim(-2, 42)
plt.savefig('/tmp/hh_data_all.pdf')
plt.figure()
plt.plot(V, n, 'ok', alpha=0.2)
plt.xlabel('V')
plt.ylabel('n')
np.savetxt('hh_data.txt',
np.vstack((V, m, n, h, np.array(i_inj_values))).T,
fmt='%.5f')
plt.show()
plt.savefig('/tmp/hh_data_V_n.pdf')
def trajectory_from_magnetic_field_method_ODE(self, source,t=None):
gamma = source.Lorentz_factor()
B=source.magnetic_field
# trajectory =
# [t........]
# [ X/c......]
# [ Y/c ......]
# [ Z/c ......]
# [ Vx/c .....]
# [ Vy/c .....]
# [ Vz/c .....]
# [ Ax/c .....]
# [ Ay/c .....]
# [ Az/c .....]
if t is None :
time_calc = source.construct_times_vector(initial_contition=self.initial_condition,Nb_pts=self.Nb_pts)
else :
self.Nb_pts=len(t)
time_calc =t
trajectory = np.zeros((10, self.Nb_pts))
trajectory[0]=time_calc
cst = -codata.e / (codata.m_e * gamma)
initial_condition_for_ODE=self.initial_condition
#TODO rtol et a tol modifiable
#atol=np.array([1e-10,1e-10,1e-10,1e-10,1e-10,1e-10])
rtol=source.rtol_for_ODE_method()
atol=source.atol_for_ODE_method()
res = odeint(func=fct_ODE_magnetic_field,y0=initial_condition_for_ODE, t=trajectory[0],
args=(cst,B.Bx,B.By,B.Bz),rtol=rtol,atol=atol,mxstep=5000,full_output=True)
traj = res[0]
info = res[1]
print("1 : nonstiff problems, Adams . 2: stiff problem, BDF")
print(info.get('mused'))
traj = np.transpose(traj)
trajectory[4] = traj[0]
trajectory[5] = traj[1]
trajectory[6] = traj[2]
trajectory[1] = traj[3]
trajectory[2] = traj[4]
trajectory[3] = traj[5]
trajectory[7] = - cst * B.By(trajectory[3], trajectory[2],trajectory[1]) * trajectory[6]
trajectory[9] = cst * B.By(trajectory[3], trajectory[2],trajectory[1]) * trajectory[4]
T=self.create_from_array(trajectory)
T.multiply_by((1.0/codata.c))
return T
def make_banded_jacobian(self):
'''returns a banded jacobian list (in odeint's format), along with mu and ml parameters'''
assert not self.sparse
if self.dense_a_matrix is None:
self.dense_a_matrix = self.sparse_a_matrix.toarray()
matrix = self.dense_a_matrix
# first find the values of mu and ml
dims = matrix.shape[0]
assert dims == matrix.shape[1]
mu = 0
ml = 0
for row in xrange(dims):
for col in xrange(dims):
if matrix[row][col] != 0:
if col > row:
dif = col - row
mu = max(mu, dif)
else:
dif = row - col
ml = max(ml, dif)
banded = []
for yoffset in xrange(-mu, ml+1):
row = []
for diag in xrange(dims):
x_index = diag
y_index = diag + yoffset
if y_index < 0 or y_index >= dims:
row.append(0.0)
else:
row.append(matrix[y_index][x_index])
banded.append(row)
return (banded, mu, ml)
def Lorenz(x0, sigma, rho, beta, time):
"""
This small function runs a simulation of the Lorenz system.
Inputs
------
x0 : numpy array containing the initial condition.
sigma, rho, beta : parameters of the Lorenz system.
time : numpy array for the evaluation of the state of
the Lorenz system at some given time instants.
Outputs
-------
x : numpy two-dimensional array.
State vector of the vector for the time instants
specified in time.
xdot : corresponding derivatives evaluated using
central differences.
"""
def dynamical_system(y,t):
dy = np.zeros_like(y)
dy[0] = sigma*(y[1]-y[0])
dy[1] = y[0]*(rho - y[2]) - y[1]
dy[2] = y[0]*y[1] - beta*y[2]
return dy
x = odeint(dynamical_system,x0,time)
dt = time[1]-time[0]
from sparse_identification.utils import derivative
xdot = derivative(x,dt)
return x, xdot
def oderest(sizes,x,u,pi,t,constants,boundary,restrictions):
print("\nIn oderest.")
# get sizes
N = sizes['N']
n = sizes['n']
m = sizes['m']
p = sizes['p']
print("Calc phi...")
# calculate phi and psi
phi = calcPhi(sizes,x,u,pi,constants,restrictions)
# aux: phi - dx/dt
aux = phi - ddt(sizes,x)
# get gradients
print("Calc grads...")
Grads = calcGrads(sizes,x,u,pi,constants,restrictions)
phix = Grads['phix']
lam = 0*x
B = numpy.zeros((N,m))
C = numpy.zeros(p)
print("Integrating ODE for A...")
# integrate equation for A:
A = odeint(calcADotOdeRest,numpy.zeros(n),t,args= (t,phix,aux))
#optPlot = dict()
#optPlot['mode'] = 'var'
#plotSol(sizes,t,A,B,C,constants,restrictions,optPlot)
print("Calculating step...")
alfa = calcStepOdeRest(sizes,t,x,u,pi,A,B,C,constants,boundary,restrictions)
nx = x + alfa * A
print("Leaving oderest with alfa =",alfa)
return nx,u,pi,lam,mu
# ##################
# MAIN SEGMENT:
# ##################
def S_eta(L_vec, c2_val):
c2 = c2_val
traj = np.zeros([x_ang.shape[0]*2+1,L_vec.shape[0]])
n=0
for L in L_vec:
slope = (m*(m+1)+c2-L)/(m+1)/2.
z0 = [1+step*slope,slope]
z = odeint(f,z0,x_ang,args=(c2,L,m))
temp=1-pow(x_ang,2.0)
temp=pow(temp,m/2.)
zz=temp*z[:,0]
first_zz = np.array([1])
zz=np.append(first_zz, zz)
sloper = -(m*(m+1)+c2-L)/(m+1)/2.
z0r = [1-step*sloper,sloper]
zr = odeint(f,z0r,x_angr,args=(c2,L,m))
zzr=temp*zr[:,0]
zzr=np.append(first_zz, zzr)
traj[:,n] = np.append(zz,zzr[::-1][1:])
n += 1
x_tot = np.arange(-1,1.+step,step)
figure = plt.figure(figsize=(12, 11))
plt.plot(x_tot,traj, linewidth=2.0, label = '')
plt.ylabel('S($\\eta$)')#, fontdict=font)
plt.xlabel('$\\eta$')#, fontdict=font)
#plt.xlim(-1.0,1.0)
#plt.ylim(0.3,1.0)
plt.locator_params(axis='x', nbins=10)
plt.locator_params(axis='y', nbins=10)
plt.tick_params(axis='x', pad=15)
#plt.legend(loc=1)
plt.show()
plt.close()
def R_xi(E_vec,L0):
traj = np.zeros([x_rad.shape[0]+1,E_vec.shape[0]])
n=0
for E in E_vec:
c2 = D*D*E/4.
L = newton(newton_ang_func,L0,args=(c2,m),tol=1e-8,maxiter=200)
slope = -(-L+m*(m+1.)+2.*D+c2)/(2.*(m+1.))
z0 = [1+step*slope,slope]
z = odeint(g,z0,x_rad,args=(c2,L,D,m))
temp=pow(x_rad,2.0)-1.
temp=pow(temp,m/2.)
zz=temp*z[:,0]
first_zz = np.array([1])
zz=np.append(first_zz, zz)
traj[:,n] = zz
n += 1
xt = np.append(np.array([1]),x_rad)
figure = plt.figure(figsize=(12, 11))
plt.plot(xt,traj, linewidth=2.0, label = '')
plt.ylabel('R($\\xi$)')#, fontdict=font)
plt.xlabel('$\\xi$')#, fontdict=font)
#plt.xlim(1.0,10.0)
#plt.ylim(-1.0,1.0)
plt.locator_params(axis='x', nbins=10)
plt.locator_params(axis='y', nbins=10)
plt.tick_params(axis='x', pad=15)
#plt.legend(loc=1)
plt.show()
plt.close()
