我们从Python开源项目中，提取了以下15个代码示例，用于说明如何使用scipy.integrate.ode()。
def __init__(self, settings): settings.update(output_dim=4) super().__init__(settings) self.a_mat, self.b_mat, self.c_mat = linearise_system( self._settings["steady state"], self._settings["steady tau"]) self.L = pm.place_siso(self.a_mat.T, self.c_mat.T, self._settings["poles"]).T self.step_width = self._settings["tick divider"] * self._settings["step width"] self.output = np.array(self._settings["initial state"], dtype=float) self.solver = ode(self.linear_state_function) self.solver.set_integrator("dopri5", method=self._settings["Method"], rtol=self._settings["rTol"], atol=self._settings["aTol"]) self.solver.set_initial_value(self._settings["initial state"]) self.nextObserver_output = self.output
def __init__(self, settings): Solver.__init__(self, settings) # setup solver if hasattr(self._model, "jacobian"): self._solver = ode(self._model.state_function, jac=self._model.jacobian) else: self._solver = ode(self._model.state_function) self._solver.set_integrator(self._settings["Mode"], method=self._settings["Method"], rtol=self._settings["rTol"], atol=self._settings["aTol"], max_step=self._settings["step size"] ) self._solver.set_initial_value(np.atleast_1d(self._model.initial_state), t=self._settings["start time"])
def odeint(fun, t, y0, integrator=0, method=0, rtol=0, atol=1e-12): s = ode(fun) integrators = ['vode', 'zvode', 'lsoda', 'dopri5', 'dop853'] methods = ['adams', 'bdf'] s.set_integrator(integrators[0], method=methods[0], order=10, rtol=rtol, atol=atol, with_jacobian=False) t0 = t[0] dt = t[1]-t0 y = [y0] s.set_initial_value(y0, t0) while s.successful() and s.t < t[-1]: s.integrate(s.t+dt) y.append(s.y) y = array(y) return y
def implicit_black_box(propensities,V,X,w,h,deter_vector,stoc_positions, positions, valid,deriv): from scipy.integrate import ode #pdb.set_trace() deter_ode = ode(f).set_integrator('lsoda',method='bdf', with_jacobian=False) deter_ode.set_initial_value(X[deter_vector,:].flatten(), 0).set_f_params([propensities,V,X,deter_vector,stoc_positions, positions, valid,w]) #pdb.set_trace() while deter_ode.successful() and deter_ode.t < h: deter_ode.integrate(h) #print("Black Box: \n"+ str(deter_ode.y)) #print("iterator : \n:"+str(next_X[deter_vector,:])) X[deter_vector,:] = deter_ode.y.reshape(( np.sum(deter_vector),X.shape[1])) return X
def perform_action(self, action, perturb_params=False, p_lambda1=0, p_lambda2=0, p_k1=0, \ p_k2=0, p_f=0, p_m1=0, p_m2=0, p_lambdaE=0, p_bE=0, p_Kb=0, p_d_E=0, p_Kd=0, **kw): """Perform the specifed action and upate the environment. Arguments: action -- action to be taken Keyword Arguments: perturb_params -- boolean indicating whether to perturb dynamics (default: False) p_lambda1 -- hidden parameter (default: 0) p_lambda2 -- hidden parameter (default: 0) p_k1 -- hidden parameter (default: 0) p_k2 -- hidden parameter (default: 0) p_f -- hidden parameter (default: 0) p_m1 -- hidden parameter (default: 0) p_m2 -- hidden parameter (default: 0) p_lambdaE -- hidden parameter (default: 0) p_bE -- hidden parameter (default: 0) p_Kb -- hidden parameter (default: 0) p_d_E -- hidden parameter (default: 0) p_Kd -- hidden parameter (default: 0) """ self.t += 1 self.action = action eps1, eps2 = self.eps_values_for_actions[action] r = ode(self.model_derivatives).set_integrator('vode',nsteps=10000,method='bdf') t0 = 0 deriv_args = (eps1, eps2, perturb_params, p_lambda1, p_lambda2, p_k1, p_k2, p_f, p_m1, p_m2, p_lambdaE, p_bE, p_Kb, p_d_E, p_Kd) r.set_initial_value(self.state, t0).set_f_params(deriv_args) self.state = r.integrate(self.dt) reward = self.calc_reward(action=action) return reward, self.observe()
def line(date_dt, pos0, dt=1e2): """ ??? https://www.mathworks.com/matlabcentral/fileexchange/34388-international-geomagnetic-reference-field--igrf--model/content/igrfline.m """ r_km0 = pos0.radius / 1e3 theta0 = math.radians(pos0.geocentricLatitude) phi0 = math.radians(pos0.longitude) y0, t0 = [r_km0, theta0, phi0], 0 tracer = ode(differential).set_integrator('dopri5') tracer.set_initial_value(y0, t0) pos_path = [PyPosition(math.degrees(theta0), math.degrees(phi0), r_km0 * 1e3, PyPosition.CoordinateSystem['geocentric'])] while True: assert tracer.successful() r_km_i, theta_i, phi_i = tracer.integrate(tracer.t + dt) pos = PyPosition(math.degrees(theta_i), math.degrees(phi_i), r_km_i * 1e3, PyPosition.CoordinateSystem['geocentric']) #print(pos.height) if pos.height < 0: break pos_path.append(pos) return pos_path
def create_solver(dydt): solver = ode(dydt).set_integrator('lsoda', method='bdf', rtol=1e-2) return solver
def __init__(self, model, jac=False, integrator='vode', **integrator_kwargs): def func_scipy_proxy(t, U, fields, pars, hook): fields.fill(U) fields, pars = hook(t, fields, pars) return model.F(fields, pars) def jacob_scipy_proxy(t, U, fields, pars, hook): fields.fill(U) fields, pars = hook(t, fields, pars) return model.J(fields, pars, sparse=False) self._solv = ode(func_scipy_proxy, jac=jacob_scipy_proxy if jac else None) self._solv.set_integrator(integrator, **integrator_kwargs)
def implicit_black_box(propensities,V,X,w,h,deter_vector,stoc_positions, positions, valid,deriv): # Adjustment for systems reaching steady state temp = derivative_G(propensities,V,X,w,deter_vector,stoc_positions,positions,valid) #pdb.set_trace() valid_adjust_pos = np.where(np.sum(np.abs(temp),axis=0) < 1e-10,True,False) valid_adjust = valid[:,:] valid_adjust[valid_adjust_pos,:] = False #print(" Reached Steady State %d"%(np.sum(valid_adjust_pos))) from scipy.integrate import ode #pdb.set_trace() deter_ode = ode(f).set_integrator('lsoda',method='adams', with_jacobian=False) deter_ode.set_initial_value(X[deter_vector,:].flatten(), 0).set_f_params([propensities,V,X,deter_vector,stoc_positions, positions, valid_adjust,w]) #pdb.set_trace() while deter_ode.successful() and deter_ode.t < h: deter_ode.integrate(h) #print("Black Box: \n"+ str(deter_ode.y)) #print("iterator : \n:"+str(next_X[deter_vector,:])) X[deter_vector,:] = deter_ode.y.reshape(( np.sum(deter_vector),X.shape[1])) ## Another adjust to compensate for non negative X = np.where(X < 0.0,0.0,X) return X
def implicit_black_box(propensities,V,X,w,h,deter_vector,stoc_positions, positions, valid,jac): # Adjustment for systems reaching steady state """ temp = derivative_G(propensities,V,X,w,deter_vector,stoc_positions,positions,valid,jac) #pdb.set_trace() valid_adjust_pos = np.where(np.sum(np.abs(temp),axis=0) < 1e-10,True,False) valid_adjust = valid[:,:] valid_adjust[valid_adjust_pos,:] = False print(" Reached Steady State %d"%(np.sum(valid_adjust_pos))) """ from scipy.integrate import ode #pdb.set_trace() deter_ode = ode(f).set_integrator('vode',method='bdf', with_jacobian=False) deter_ode.set_initial_value(X[deter_vector,:].flatten(), 0) #deter_ode.set_f_params([propensities,V,X,deter_vector,stoc_positions, positions,valid_adjust,w,jac]) deter_ode.set_f_params([propensities,V,X,deter_vector,stoc_positions, positions,valid,w,jac]) #pdb.set_trace() while deter_ode.successful() and deter_ode.t < h: deter_ode.integrate(h) #print("Black Box: \n"+ str(deter_ode.y)) #print("iterator : \n:"+str(next_X[deter_vector,:])) X[deter_vector,:] = deter_ode.y.reshape(( np.sum(deter_vector),X.shape[1])) # Another adjust to compensate for non negative X = np.where(X < 0.0,0.0,X) return X
def setup(self): self.integrator = ode(self.constitutive_relation) self.integrator.set_integrator("lsoda") # Main integrator function. Output vector: # 0: friction, 1+: state parameters
def _simulate_with_scipy(self, initcon, tsim, switchpts, varying_inputs, integrator, integrator_options): scipyint = \ scipy.ode(self._rhsfun).set_integrator(integrator, **integrator_options) scipyint.set_initial_value(initcon, tsim[0]) profile = np.array(initcon) i = 1 while scipyint.successful() and scipyint.t < tsim[-1]: # check if tsim[i-1] is a switching time and update value if tsim[i - 1] in switchpts: for v in self._siminputvars.keys(): if tsim[i - 1] in varying_inputs[v]: p = self._templatemap[self._siminputvars[v]] p.set_value(varying_inputs[v][tsim[i - 1]]) profilestep = scipyint.integrate(tsim[i]) profile = np.vstack([profile, profilestep]) i += 1 if not scipyint.successful(): raise DAE_Error("The Scipy integrator %s did not terminate " "successfully." % integrator) return [tsim, profile]
def _simulate_with_casadi_no_inputs(self, initcon, tsim, integrator, integrator_options): # Old way (10 times faster, but can't incorporate time # varying parameters/controls) xalltemp = [self._templatemap[i] for i in self._diffvars] xall = casadi.vertcat(*xalltemp) odealltemp = [value(self._rhsdict[i]) for i in self._derivlist] odeall = casadi.vertcat(*odealltemp) dae = {'x': xall, 'ode': odeall} if len(self._algvars) != 0: zalltemp = [self._templatemap[i] for i in self._simalgvars] zall = casadi.vertcat(*zalltemp) algalltemp = [value(i) for i in self._alglist] algall = casadi.vertcat(*algalltemp) dae['z'] = zall dae['alg'] = algall integrator_options['grid'] = tsim integrator_options['output_t0'] = True F = casadi.integrator('F', integrator, dae, integrator_options) sol = F(x0=initcon) profile = sol['xf'].full().T if len(self._algvars) != 0: algprofile = sol['zf'].full().T profile = np.concatenate((profile, algprofile), axis=1) return [tsim, profile]
def calculateTurbulent_cf(self, S_0): """ Return the respective streamline given a starting position. point : [x y z] defining the starting point. MUST be one of the starting points used to generate the streamlines in the first place. """ self.getProperties(S_0, 'turbulent') cf_0 = 0.455 / ( self.Q**2. * log(0.06/self.Q * self.localReynolds * self.mu/self.mu_wall * (self.TwOnT)**-0.5 )**2.) cf_0 = 0.00165 self.chi_0 = (2. / cf_0)**0.5 self.streamlineCfs = [] self.streamlineXs = [] r1 = ode(self.whiteChristoph1).set_integrator('dopri5', atol=1e-2, rtol=1e-2, ) r1.set_initial_value(self.chi_0, S_0) r2 = ode(self.whiteChristoph2).set_integrator('dopri5', atol=1e-2, rtol=1e-2, ) r2.set_initial_value(self.chi_0, S_0) dt = 0.01 # Integrate along the streamline while r1.t < self.parameterisedStreamline[-1] and r2.t < self.parameterisedStreamline[-1]: # print r1.t + dt # print r2.t + dt r1.integrate(r1.t+dt) if self.ReynoldsStar < 0: self.streamlineCfs.append(2. / r1.y**2.) self.streamlineXs.append(r1.t) else: # print "r2" r2.integrate(r2.t+dt) if r1.y/self.chi_max < 0.36: self.streamlineCfs.append(2. / r1.y**2.) self.streamlineXs.append(r1.t) elif r2.y/self.chi_max > 0.36: self.streamlineCfs.append(2. / r2.y**2.) self.streamlineXs.append(r2.t) else: print "Something stuffed up" return self.streamlineCfs, self.streamlineXs
def line(self, x0): """Returns the field line passing through x0. Refs: http://folk.uib.no/fcihh/seminar/lec1.pdf and lect2.pdf http://numbercrunch.de/blog/2013/05/visualizing-streamlines/ and especially: "Electric field lines don't work", http://scitation.aip.org/content/aapt/journal/ajp/64/6/10.1119/1.18237 """ if None in [XMIN, XMAX, YMIN, YMAX]: raise ValueError('Domain must be set using init().') # Set up integrator for the field line streamline = lambda t, y: list(self.direction(y)) solver = ode(streamline).set_integrator('vode') # Initialize the coordinate lists x = [x0] # Integrate in both the forward and backward directions dt = 0.008 # Solve in both the forward and reverse directions for sign in [1, -1]: # Set the starting coordinates and time solver.set_initial_value(x0, 0) # Integrate field line over successive time steps while solver.successful(): # Find the next step solver.integrate(solver.t + sign*dt) # Save the coordinates if sign > 0: x.append(solver.y) else: x.insert(0, solver.y) # Check if line connects to a charge flag = False for c in self.charges: if c.is_close(solver.y): flag = True break # Terminate line at charge or if it leaves the area of interest if flag or not (XMIN < solver.y[0] < XMAX) or \ not YMIN < solver.y[1] < YMAX: break return FieldLine(x)
