我们从Python开源项目中，提取了以下40个代码示例，用于说明如何使用scipy.integrate.odeint()。
def likelihood(parameter_vector): parameter_vector = 10**np.array(parameter_vector) #Solve ODE system given parameter vector yout = odeint(odefunc, y0, tspan, args=(parameter_vector,)) cout = yout[:, 2] #Calculate log probability contribution given simulated experimental values. logp_ctotal = np.sum(like_ctot.logpdf(cout)) #If simulation failed due to integrator errors, return a log probability of -inf. if np.isnan(logp_ctotal): logp_ctotal = -np.inf return logp_ctotal # Add vector of rate parameters to be sampled as unobserved random variables in DREAM with uniform priors.
def velocity(stateVec, t): """ return the velocity field of Lorentz system. stateVec : the state vector in the full space. [x, y, z] t : time is used since odeint() requires it. """ x = stateVec[0] y = stateVec[1] z = stateVec[2] # complete the flowing 3 lines. vx = G_sigma*(y - x) vy = G_rho*x - y - x*z vz = x*y - G_b*z return np.array([vx, vy, vz])
def velocity(stateVec, t): """ velocity in the full state space. stateVec: state vector [x1, y1, x2, y2] t: just for convention of odeint, not used. return: velocity at stateVec. Dimension [1 x 4] """ x1 = stateVec[0] y1 = stateVec[1] x2 = stateVec[2] y2 = stateVec[3] r2 = x1**2 + y1**2 velo = np.array([(G_mu1-r2) * x1 + G_c1 * (x1*x2 + y1*y2), (G_mu1-r2) * y1 + G_c1 * (x1*y2 - x2*y1), x2 + y2 + x1**2 - y1**2 + G_a2 * x2 * r2, -x2 + y2 + 2.0 * x1 * y1 + G_a2 * y2 * r2]) return velo
def solve_coupled_constecc_solution(F0, e0, phase0, mc, t): """ Compute the solution to the coupled system of equations from from Peters (1964) and Barack & Cutler (2004) at a given time. :param F0: Initial orbital frequency [Hz] :param mc: Chirp mass of binary [Solar Mass] :param t: Time at which to evaluate solution [s] :returns: (F(t), phase(t)) """ y0 = np.array([F0, phase0]) y, infodict = odeint(get_coupled_constecc_eqns, y0, t, args=(mc,e0), full_output=True) if infodict['message'] == 'Integration successful.': ret = y else: ret = 0 return ret
def step(self, action): action = self._discrete_actions[action[0]] sa = np.append(self._state, action) new_state = odeint(self._dpds, sa, [0, self._dt]) self._state = new_state[-1, :-1] if self._state[0] < -self.max_pos or \ np.abs(self._state[1]) > self.max_velocity: reward = -1 absorbing = True elif self._state[0] > self.max_pos and \ np.abs(self._state[1]) <= self.max_velocity: reward = 1 absorbing = True else: reward = 0 absorbing = False return self._state, reward, absorbing, {}
def step(self, action): action = self._discrete_actions[action[0]] action += np.random.uniform(-10., 10.) sa = np.append(self._state, action) new_state = odeint(self._dpds, sa, [0, self._dt]) self._state = new_state[-1, :-1] self._state[0] = self._range_pi(self._state[0]) if np.abs(self._state[0]) > np.pi / 2.: reward = -1 absorbing = True else: reward = 0 absorbing = False return self._state, reward, absorbing, {}
def Wdesired(x): ''' x is the augmented variable represented in the network it obeys \tau_s x_\alpha = -x_\alpha + Wdesired_\alpha(x) x is related to the original augmented variable \tilde{x} by x_\alpha = varFactors_\alpha \tilde{x}_\alpha where varFactors_alpha = angleFactor | velocityFactor | torqueFactor now, original augmented variable obeys \dot{\tilde{x}}=f(\tilde{x}) so, we have Wdesired_\alpha(x) = \tau_s * varFactor_alpha * f_\alpha(\tilde{x}) + x ''' # \tilde{x} (two zeroes at x[N:N+N//2] are ignored xtilde = x/varFactors if XY: angles = armAngles(xtilde[:N]) else: angles = xtilde[:N//2] # f(\tilde{x}), \dot{u} part is not needed qdot,dqdot = evolveFns(angles,xtilde[Nobs-N//2:Nobs],xtilde[Nobs:],XY,dt) # returns deltaposn if XY else deltaangles # \tau_s * varFactors_alpha * f_\alpha(\tilde{x}) + x return np.append(np.append(qdot,dqdot),np.zeros(N//2))*varFactors*tau + x # integral on torque u also # VERY IMP to compensate for synaptic decay on torque #return np.append( np.append(qdot,dqdot)*tau*varFactors[:Nobs] + x[:Nobs], np.zeros(N//2) ) # normal synaptic decay on torque u ##### For the reference, choose EITHER robot simulation rateEvolveProbe above ##### OR evolve Wdesired inverted / evolveFns using odeint as below -- both should be exactly same
def reaction(y0, t, k0, k0r): """ Wrapper for the reaction. It receives an `np.array` `y0` with the initial concentrations and an `np.array` `t` with timestamps (it should also include `0`). This function solves the corresponding ODE system and returns an `np.array` `Y` in which each column represents a chemical species and each line a timestamp. """ def dydt(y, t): return np.array([ -1*v_0(y[0], y[1]), +1*v_0(y[0], y[1]), ]) # A <-> B def v_0(A, B): return k0 * A**1 - k0r * B**1 return odeint(dydt, y0, t) # Reaction rates: # A <-> B
def reaction(y0, t, k0, k0r): """ Wrapper for the reaction. It receives an `np.array` `y0` with the initial concentrations and an `np.array` `t` with timestamps (it should also include `0`). This function solves the corresponding ODE system and returns an `np.array` `Y` in which each column represents a chemical species and each line a timestamp. """ def dydt(y, t): return np.array([ -2*v_0(y[0], y[1]), +1*v_0(y[0], y[1]), ]) # 2*A <-> C def v_0(A, C): return k0 * A**2 - k0r * C**1 return odeint(dydt, y0, t) # Reaction rates: # 2*A <-> C
def reaction(y0, t, k0, k0r): """ Wrapper for the reaction. It receives an `np.array` `y0` with the initial concentrations and an `np.array` `t` with timestamps (it should also include `0`). This function solves the corresponding ODE system and returns an `np.array` `Y` in which each column represents a chemical species and each line a timestamp. """ def dydt(y, t): return np.array([ -1*v_0(y[1], y[0], y[2]), +2*v_0(y[1], y[0], y[2]), +2*v_0(y[1], y[0], y[2]), ]) # S <-> 2*A + 2*C def v_0(A, S, C): return k0 * S**1 - k0r * A**2 * C**2 return odeint(dydt, y0, t) # Reaction rates: # S <-> 2*A + 2*C
def raw_sim_one(start, steps, dy_data, settings, include_step_zero=False): ''' simulate from a single point at the given times return an nparray of states at those times, possibly excluding time zero ''' start.shape = (dy_data.num_dims,) times = np.linspace(0, settings.step * steps, num=steps+1) der_func = dy_data.make_der_func() jac_func, max_upper, max_lower = dy_data.make_jac_func() sim_tol = settings.sim_tol result = odeint(der_func, start, times, Dfun=jac_func, col_deriv=True, \ mxstep=int(1e8), mu=max_upper, ml=max_lower, \ atol=sim_tol, rtol=sim_tol) if not include_step_zero: result = result[1:] return result
def sim(start, der_func, time_amount, num_steps, quick): 'simulate for some fixed time, and return the resultant (states, times) tuple' tol = 1.49012e-8 if not quick: tol /= 1e5 # more accurate simulation times = np.linspace(0, time_amount, num=1 + num_steps) states = odeint(der_func, start, times, col_deriv=True, rtol=tol, atol=tol, mxstep=50000) states = states[1:] times = times[1:] assert len(states) == num_steps assert len(times) == num_steps return (states, times)
def run(self,t=1000,ts=100000,plot=False): """ Run it for time t with ts number of time steps Outputs the concentration profile for the run default values are 100000 and 1000000 """ t_index = np.linspace(0, t, ts) y = self.concentrations output = odeint(odes, y, t_index, args= (self.reactions,self.rates)) self.concentrations = output[-1,:] if plot: for i in xrange(output.shape[1]): label = self.species[i] plt.plot(t_index, output[:, i], label=label) plt.legend(loc='best') plt.xlabel('t') plt.title('Concentration Profile') plt.grid() plt.show() return output
def solveUnst(self,t): solution = [None]*len(self.mapping) # This has the unsteady part # Solver, just live and let live for ix, (i,k) in enumerate(self.mapping): solution[ix] = self.b[i][k] soln = odeint(self.rUnst, solution, t,hmax=(t[-1]-t[0])/len(t), \ rtol = 1e-4, atol = 1e-4) # final update self.updateUnst(t[-1]) self.update(soln[-1,:]) # Lets collect all the steps fullSolution = dict([(b.name + '_'+s,[]) for b in self.b for s in b.state]) for j in range(0,len(t)): for ix, (i,k) in enumerate(self.mapping): fullSolution[self.b[i].name+'_'+k].append(soln[j,ix]) fullSolution['t'] = t return fullSolution
def newton_rad_func(E_val,D,m,L0): E = E_val c2 = D*D*E/4. L = newton(newton_ang_func,L0,args=(c2,m),tol=1e-8,maxiter=500) 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) return z[:,1][-1]
def plottrajectories(fs, x0, t=np.linspace(1,400,10000), **kw): """ plots trajectory of the solution f -- must accept an array of X and t=0, and return a 2D array of \dot y and \dot x x0 -- vector Example ======= plottrajectories(lambda X,t=0:array([ X[0] - X[0]*X[1] , -X[1] + X[0]*X[1] ]), [ 5,5], color='red') """ x0 = np.array(x0) #f = lambda X,t=0: array(fs(X[0],X[1])) #fa = lambda X,t=0:array(fs(X[0],X[1])) X = integrate.odeint( fs, x0, t) plt.plot(X[:,0], X[:,1], **kw)
def solve(self): self.psi = odeint(self.RichardsEquation, self.psi0, self.t, args=( self.dz, self.n, self.p, self.qTop, self.qBot, self.psiTop, self.psiBot), mxstep=500) self.psi0 = self.psi[-1, :]
def calculate_solution(self): solution = odeint(func=self.set_up_equation, y0=self.initial_condition, t=self.steps, args=(self.coefficient_matrix, self.steps)) return solution
def update(self, dt, F, M): # limit thrust and Moment L = params.arm_length r = params.r prop_thrusts = params.invA.dot(np.r_[np.array([[F]]), M]) prop_thrusts_clamped = np.maximum(np.minimum(prop_thrusts, params.maxF/4), params.minF/4) F = np.sum(prop_thrusts_clamped) M = params.A[1:].dot(prop_thrusts_clamped) self.state = integrate.odeint(self.state_dot, self.state, [0,dt], args = (F, M))[1]
def TC_Simulate(Mode,initialCondition,time_bound): time_step = 0.05; time_bound = float(time_bound) initial = [float(tmp) for tmp in initialCondition] number_points = int(np.ceil(time_bound/time_step)) t = [i*time_step for i in range(0,number_points)] if t[-1] != time_step: t.append(time_bound) y_initial = initial[0] if Mode == 'On': rate = 0.1 elif Mode == 'Off': rate = -0.1 else: print('Wrong Mode name!') sol = odeint(thermo_dynamic,y_initial,t,args=(rate,),hmax = time_step) # Construct the final output trace = [] for j in range(len(t)): #print t[j], current_psi tmp = [] tmp.append(t[j]) tmp.append(sol[j,0]) trace.append(tmp) return trace # sol = TC_Simulate('Off',[60],10) # print(sol)
def integrator(init_x, dt, nstp): """ The integator of the Lorentz system. init_x: the intial condition dt : time step nstp: number of integration steps. return : a [ nstp x 3 ] vector """ state = odeint(velocity, init_x, np.arange(0, dt*nstp, dt)) return state
def Jacobian(ssp, dt, nstp): """ Jacobian function for the trajectory started on ssp, evolved for time t Inputs: ssp: Initial state space point. dx1 NumPy array: ssp = [x, y, z] t: Integration time Outputs: J: Jacobian of trajectory f^t(ssp). dxd NumPy array """ #CONSTRUCT THIS FUNCTION #Hint: See the Jacobian calculation in CycleStability.py #J = None Jacobian0 = np.identity(3) # COMPLETE THIS LINE. HINT: Use np.identity(DIMENSION) #Initial condition for Jacobian integral is a d+d^2 dimensional matrix #formed by concatenation of initial condition for state space and the #Jacobian: sspJacobian0 = np.zeros(3 + 3 ** 2) # Initiate sspJacobian0[0:3] = ssp # First 3 elemenets sspJacobian0[3:] = np.reshape(Jacobian0, 9) # Remaining 9 elements tInitial = 0 # Initial time tFinal = dt*nstp # Final time Nt = nstp # Number of time points to be used in the integration tArray = np.linspace(tInitial, tFinal, Nt) # Time array for solution sspJacobianSolution = odeint(JacobianVelocity, sspJacobian0, tArray) xt = sspJacobianSolution[:, 0] # Read x(t) yt = sspJacobianSolution[:, 1] # Read y(t) zt = sspJacobianSolution[:, 2] # Read z(t) #Read the Jacobian for the periodic orbit: J = sspJacobianSolution[-1, 3:].reshape((3, 3)) return J
def integrator(init_state, dt, nstp): """ integrate two modes system in the full state sapce. init_state: initial state [x1, y1, x2, y2] dt: time step nstp: number of time step """ states = odeint(velocity, init_state, np.arange(0, dt*nstp, dt)) return states
def solve_coupled_ecc_solution(F0, e0, gamma0, phase0, mc, q, t): """ Compute the solution to the coupled system of equations from from Peters (1964) and Barack & Cutler (2004) at a given time. :param F0: Initial orbital frequency [Hz] :param e0: Initial orbital eccentricity :param gamma0: Initial angle of precession of periastron [rad] :param mc: Chirp mass of binary [Solar Mass] :param q: Mass ratio of binary :param t: Time at which to evaluate solution [s] :returns: (F(t), e(t), gamma(t), phase(t)) """ y0 = np.array([F0, e0, gamma0, phase0]) y, infodict = odeint(get_coupled_ecc_eqns, y0, t, args=(mc,q), full_output=True) if infodict['message'] == 'Integration successful.': ret = y else: ret = 0 return ret
def matevolve2(y,t): ''' the reference y is only N-dim i.e. (q,dq), not 2N-dim, as inpfn is used directly as reference for torque u ''' ## invert the nengo function transformation with angleFactor, tau, +x, etc. in Wdesired() ## -- some BUG, inversion is not working correctly #xfull = np.append(y,inpfn(t))*varFactors #return ( (Wdesired(xfull)[:N]/tau/varFactors[:N] - xfull[:N]) \ # if (t%Tperiod)<(Tperiod-Tclamp) else -x/tau ) # instead of above, directly use evolveFns() ######### DOESN'T WORK: should only use armAngles() with valid posn, not all posn-s are valid for an arm!!! ########### if XY: angles = armAngles(y[:N]) else: angles = y[:N//2] # evolveFns returns deltaposn if XY else deltaangles if trialClamp: return ( evolveFns( angles, y[Nobs-N//2:Nobs], inpfn(t), XY, dt).flatten()\ if (t%Tperiod)<(Tperiod-Tclamp) else -y/dt ) else: return evolveFns( angles, y[Nobs-N//2:Nobs], inpfn(t), XY, dt).flatten() ##### uncomment below to override rateEvolveProbe by matevolve2-computed Wdesired-inversion / evolveFns-evolution, as reference signal #trange = np.arange(0.0,Tmax,dt) #y = odeint(matevolve2,0.001*np.ones(N),trange,hmax=dt) # set hmax=dt, to avoid adaptive step size # # some systems (van der pol) have origin as a fixed pt # # hence start just off-origin #rateEvolveProbe = y # only copies pointer, not full array (no need to use np.copy() here) ### ### Reference evolution used when copycat layer is not used for reference ### ### # scale the output of the robot simulation or odeint to cover the representation range of the network # here I scale by angle/velocity factors, below at nodeIn I scale by torque factors.
def reaction(y0, t, k0, k0r, k1, k1r): """ Wrapper for the reaction. It receives an `np.array` `y0` with the initial concentrations and an `np.array` `t` with timestamps (it should also include `0`). This function solves the corresponding ODE system and returns an `np.array` `Y` in which each column represents a chemical species and each line a timestamp. """ def dydt(y, t): return np.array([ -2*v_0(y[0], y[2], y[1])-1*v_1(y[0], y[2], y[3]), -1*v_0(y[0], y[2], y[1]), +1*v_0(y[0], y[2], y[1])-1*v_1(y[0], y[2], y[3]), +1*v_1(y[0], y[2], y[3]), ]) # 2*A + B <-> C def v_0(A, C, B): return k0 * A**2 * B**1 - k0r * C**1 # A + C <-> D def v_1(A, C, D): return k1 * A**1 * C**1 - k1r * D**1 return odeint(dydt, y0, t) # Reaction rates: # 2*A + B <-> C
def reaction(y0, t, k0, k0r, k1, k1r, k2, k2r): """ Wrapper for the reaction. It receives an `np.array` `y0` with the initial concentrations and an `np.array` `t` with timestamps (it should also include `0`). This function solves the corresponding ODE system and returns an `np.array` `Y` in which each column represents a chemical species and each line a timestamp. """ def dydt(y, t): return np.array([ -1*v_0(y[1], y[0], y[2])+1*v_1(y[3], y[0], y[2]), -1*v_0(y[1], y[0], y[2])-1*v_2(y[3], y[1]), +1*v_0(y[1], y[0], y[2])-1*v_1(y[3], y[0], y[2]), +1*v_1(y[3], y[0], y[2])+1*v_2(y[3], y[1]), ]) # E + S <-> ES def v_0(S, E, ES): return k0 * E**1 * S**1 - k0r * ES**1 # ES <-> E + P def v_1(P, E, ES): return k1 * ES**1 - k1r * E**1 * P**1 # S <-> P def v_2(P, S): return k2 * S**1 - k2r * P**1 return odeint(dydt, y0, t) # Reaction rates: # E + S <-> ES
def reaction(y0, t, k0, k0r, k1, k1r): """ Wrapper for the reaction. It receives an `np.array` `y0` with the initial concentrations and an `np.array` `t` with timestamps (it should also include `0`). This function solves the corresponding ODE system and returns an `np.array` `Y` in which each column represents a chemical species and each line a timestamp. """ def dydt(y, t): return np.array([ -1*v_0(y[1], y[0], y[2])+1*v_1(y[3], y[0], y[2]), -1*v_0(y[1], y[0], y[2]), +1*v_0(y[1], y[0], y[2])-1*v_1(y[3], y[0], y[2]), +1*v_1(y[3], y[0], y[2]), ]) # E + S <-> ES def v_0(S, E, ES): return k0 * E**1 * S**1 - k0r * ES**1 # ES <-> E + P def v_1(P, E, ES): return k1 * ES**1 - k1r * E**1 * P**1 return odeint(dydt, y0, t) # Reaction rates: # E + S <-> ES
def _transform_eigenfunction(self): eigenfunction = si.odeint(self._ff, self._init_state_vect, self._domain) return [eigenfunction[:, 0], eigenfunction[:, 1], eigenfunction[:, 2], eigenfunction[:, 3]]
def run_till(self,delta_time=0.001,eps = 10**-12,every=10, stop=100000): """ Runs till the gradients are less than eps returns time taken with array of concentrations default values of time_step and eps are 10**-3 and 10**-12 """ y = self.concentrations output = [y] t=0 t_index = np.linspace(0,every,int(every/(delta_time))) if every>10*delta_time: while True: gy = odes(y,t,self.reactions,self.rates) if (np.abs(gy) <eps).all(): break o = odeint(odes, y, t_index, args= (self.reactions,self.rates)) y = o[-1,:] t+=every output = output + o[1:] if t>stop: break else: while True: gy = odes(y,t,self.reactions,self.rates) if (np.abs(gy) <eps).all(): break y+= gy*delta_time t+=delta_time output.append(y) if t>stop: break output = np.array(output) self.concentrations =y return output,t
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))
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()
