Python integrate.ode函数代码示例

    def test_concurrent_ok(self):
        f = lambda t, y: 1.0

        for k in xrange(3):
            for sol in ('vode', 'zvode', 'dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.1)
                assert_allclose(r2.y, 0.2)

            for sol in ('dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.3)
                assert_allclose(r2.y, 0.2)

    def __init__(self, lat, lon, h, psi=0, x_earth=0, y_earth=0,
                 integrator='dopri5', use_jac=False, **integrator_params):
        """
        Initialize the equations of the chosen model and selects the
        integrator. Check `scipy.integrate.ode` to see available integrators.

        If use_jac = True the jacobian of the equations is used for the
        integration.
        """
        super().__init__(lat, lon, h, psi, x_earth, y_earth)
        # State vector must be initialized with set_initial_state_vector() method
        self.state_vector = None

        from pyfme.models.euler_flat_earth import lamceq, kaeq, kleq
        self.lamceq = lamceq

        if use_jac:
            from pyfme.models.euler_flat_earth import lamceq_jac, kaeq_jac
            jac_LM_and_AM = lamceq_jac
            jac_att = kaeq_jac
            jac_nav = None  # not implemented
        else:
            jac_LM_and_AM = None
            jac_att = None
            jac_nav = None

        self._ode_lamceq = ode(lamceq, jac=jac_LM_and_AM)
        self._ode_kaqeq = ode(kaeq, jac=jac_att)
        self._ode_kleq = ode(kleq, jac=jac_nav)

        self._ode_lamceq.set_integrator(integrator, **integrator_params)
        self._ode_kaqeq.set_integrator(integrator, **integrator_params)
        self._ode_kleq.set_integrator(integrator, **integrator_params)

def solve(si, y, infile):
    """Conducts the solution step, based on the dopri5 integrator in scipy

    :param si: the simulation info object
    :type si: SimInfo
    :param y: the solution vector
    :type y: np.ndarray
    :param infile: the imported infile module
    :type infile: imported module
    """
    n = ode(f_n).set_integrator('dopri5')
    n.set_initial_value(y0_n(si), si.timer.
                        t0.magnitude)
    n.set_f_params(si)
    th = ode(f_th).set_integrator('dopri5', nsteps=infile.nsteps)
    th.set_initial_value(y0_th(si), si.timer.t0.magnitude)
    th.set_f_params(si)
    while (n.successful()
           and n.t < si.timer.tf.magnitude
           and th.t < si.timer.tf.magnitude):
        si.timer.advance_one_timestep()
        si.db.record_all()
        n.integrate(si.timer.current_time().magnitude)
        update_n(n.t, n.y, si)
        th.integrate(si.timer.current_time().magnitude)
        update_th(th.t, n.y, th.y, si)
    return si.y

    def shooting_function_full(self, r0):
        r"""
        The function used in the shooting algorithm.

        This is the full algorithm from integrating over
        :math:`0 \le \theta \le \pi`. The difference between the
        solution and its derivative at the matching point is the
        error to be minimized.

        Parameters
        ----------

        r0 : list of float
            Initial guess for the horizon radius, as outlined above.

        Returns
        -------

        list of float
            The error at the matching point.
        """
        # First half of the horizon
        H0 = np.array([r0[0], 0.0])
        solver1 = ode(self.expansion)
        solver1.set_integrator("dopri5", atol=1.e-8, rtol=1.e-6)
        solver1.set_initial_value(H0, 0.0)
        solver1.integrate(np.pi / 2.0)
        # Second half of the horizon
        H0 = np.array([r0[1], 0.0])
        solver2 = ode(self.expansion)
        solver2.set_integrator("dopri5", atol=1.e-8, rtol=1.e-6)
        solver2.set_initial_value(H0, np.pi)
        solver2.integrate(np.pi / 2.0)

        return solver1.y - solver2.y

    def __init__(self, integrator="dopri5", model="euler_flat_earth", jac=False, **integrator_params):
        """
        Initialize the equations of the chosen model and selects the
        integrator. Check `scipy.integrate.ode` to see available integrators.

        If jac = True the jacobian of the equations is used for the
        integration.
        """

        if model not in models:
            raise ValueError(
                "The specified model is not available, please \
                             check available models in ..."
            )

        if jac:
            jac_LM_and_AM = euler_flat_earth.jac_linear_and_angular_momentum_eqs
            jac_att = euler_flat_earth.jac_linear_and_angular_momentum_eqs
            jac_nav = None  # not implemented
        else:
            jac_LM_and_AM = None
            jac_att = None
            jac_nav = None

        self._LM_and_AM_eqs = ode(euler_flat_earth.linear_and_angular_momentum_eqs, jac=jac_LM_and_AM)
        self._attitude_eqs = ode(euler_flat_earth.kinematic_angular_eqs, jac=jac_att)
        self._navigation_eqs = ode(euler_flat_earth.navigation_eqs, jac=jac_nav)

        self._LM_and_AM_eqs.set_integrator(integrator, **integrator_params)
        self._attitude_eqs.set_integrator(integrator, **integrator_params)
        self._navigation_eqs.set_integrator(integrator, **integrator_params)

        # State vector must be initialized with set_initial_values() method
        self.state_vector = None

    def s_ode(self,params):
        x0 = np.ones(3)
        s0 = np.zeros(18)
        solver_x = ode(self.dx).set_integrator('dopri5')
        solver_x.set_initial_value(x0,0).set_f_params(params,)

        solver_s = ode(self.ds).set_integrator('dopri5')
        solver_s.set_initial_value(s0,0).set_f_params(params,x0)

        dt = 1
        t1 = 250
        sensitivities = []
        for column in self.cols:
            for param in self.params:
                sensitivities.append('{0},{1}'.format(column,param))
        sol = pandas.DataFrame(np.empty((t1/dt,22)),columns=self.cols+sensitivities)
        i = 0
        #return solver_x,solver_s
        while solver_x.successful() and solver_s.successful() and solver_x.t < t1:
            solver_x.integrate(solver_x.t+dt)
            sol.iloc[i][self.cols] = solver_x.y
            solver_s.set_f_params(params,solver_x.y)
            solver_s.integrate(solver_s.t+dt)
            sol.iloc[i][sensitivities] = solver_s.y
            i += 1
        return sol

def run_trials(z, integrators, times, M):
    # Set up ODE solver
    for integrator in integrators:
        if integrator == 'dop853' or integrator == 'dopri5':
            solver = ode(system)
            solver.set_integrator(integrator, atol=1E-1, rtol=1E-3)
        elif integrator == 'bdf':
            solver = ode(system)
            solver.set_integrator('vode', method=integrator, atol=1E-1, rtol=1E-3, nsteps=2000)

        name = integrator + ' ' + str(times[-1]) + ' ' + str(M)
        try:
            z, c1, c2, t, c1_40km, c2_40km = load_variables(name)
            print "Loaded data for: " + name
        except:
            print "Starting solver: ", integrator, "with times", times
            c1, c2, c1_40km, c2_40km, t = solve(solver, c, times, integrator)
            save_variables(name, z, c1, c2, t, c1_40km, c2_40km)

        # And plot some things
        if times[-1] == 86400.0:
            labels = [str(int(time / 3600.)) + " hours" for time in times[1:]]
            plot_c1(z, c, c1, labels, name)
            plot_c2(z, c, c2, labels, name)
        elif times[-1] == 864000.0:
            plot_40km(t, c1_40km, c2_40km, name)

def _no_collapse_expect_out(num_times,expect_out):
    ##Calculates xpect.values at times tlist if no collapse ops. given
    #  
    #------------------------------------
    opt=odeconfig.options
    if odeconfig.tflag in array([1,10,11]):
        ODE=ode(odeconfig.tdfunc)
        code = compile('ODE.set_f_params('+odeconfig.string+')', '<string>', 'exec')
        exec(code)
    elif odeconfig.tflag==2:
        ODE=ode(_cRHStd)
    elif odeconfig.tflag in array([20,22]):
        ODE=ode(_tdRHStd)
    elif odeconfig.tflag==3:
        ODE=ode(_pyRHSc)
    else:
        ODE = ode(cyq_ode_rhs)
        ODE.set_f_params(odeconfig.h_data, odeconfig.h_ind, odeconfig.h_ptr)
    
    ODE.set_integrator('zvode',method=opt.method,order=opt.order,atol=opt.atol,rtol=opt.rtol,nsteps=opt.nsteps,first_step=opt.first_step,min_step=opt.min_step,max_step=opt.max_step) #initialize ODE solver for RHS
    ODE.set_initial_value(odeconfig.psi0,odeconfig.tlist[0]) #set initial conditions
    for jj in range(odeconfig.e_num):
        expect_out[jj][0]=mc_expect(odeconfig.e_ops_data[jj],odeconfig.e_ops_ind[jj],odeconfig.e_ops_ptr[jj],odeconfig.e_ops_isherm[jj],odeconfig.psi0)
    for k in range(1,num_times):
        ODE.integrate(odeconfig.tlist[k],step=0) #integrate up to tlist[k]
        if ODE.successful():
            state=ODE.y/norm(ODE.y)
            for jj in range(odeconfig.e_num):
                expect_out[jj][k]=mc_expect(odeconfig.e_ops_data[jj],odeconfig.e_ops_ind[jj],odeconfig.e_ops_ptr[jj],odeconfig.e_ops_isherm[jj],state)
        else:
            raise ValueError('Error in ODE solver')
    return expect_out #return times and expectiation values

def _no_collapse_psi_out(num_times,psi_out):
    ##Calculates state vectors at times tlist if no collapse AND no expectation values are given.
    #
    opt=odeconfig.options
    if odeconfig.tflag in array([1,10,11]):
        ODE=ode(odeconfig.tdfunc)
        code = compile('ODE.set_f_params('+odeconfig.string+')', '<string>', 'exec')
        exec(code)
    elif odeconfig.tflag==2:
        ODE=ode(_cRHStd)
    elif odeconfig.tflag in array([20,22]):
        ODE=ode(_tdRHStd)
    elif odeconfig.tflag==3:
        ODE=ode(_pyRHSc)
    else:
        ODE = ode(cyq_ode_rhs)
        ODE.set_f_params(odeconfig.h_data, odeconfig.h_ind, odeconfig.h_ptr)
        
    ODE.set_integrator('zvode',method=opt.method,order=opt.order,atol=opt.atol,rtol=opt.rtol,nsteps=opt.nsteps,first_step=opt.first_step,min_step=opt.min_step,max_step=opt.max_step) #initialize ODE solver for RHS
    ODE.set_initial_value(odeconfig.psi0,odeconfig.tlist[0]) #set initial conditions
    psi_out[0]=Qobj(odeconfig.psi0,odeconfig.psi0_dims,odeconfig.psi0_shape,'ket')
    for k in range(1,num_times):
        ODE.integrate(odeconfig.tlist[k],step=0) #integrate up to tlist[k]
        if ODE.successful():
            psi_out[k]=Qobj(ODE.y/norm(ODE.y,2),odeconfig.psi0_dims,odeconfig.psi0_shape,'ket')
        else:
            raise ValueError('Error in ODE solver')
    return psi_out

def Weightloss_Shooting():
	age, sex =  38. , 'female'
	H, BW    =  1.73, 72.7
	T		 =  12*7.   # Time frame = 3 months    
	
	PALf, EIf = 1.7, 2025
	def EI(t): return EIf
	
	def PAL(t): return PALf
	
	from solution import fat_mass, compute_weight_curve, weight_odes
	F = fat_mass(BW,age,H,sex)
	L = BW-F
	# Fix activity level and determine Energy Intake to achieve Target weight of 145 lbs = 65.90909 kg
	
	init_FL = np.array([F,L])
	EI0, EI1 = 2000, 1950					# Variable Energy Intake
	beta = 65.90909*2.2
	reltol, abstol = 1e-9,1e-8
	dim, iterate = 2, 15
	print '\nj = 1', '\nt = ', EI0
	for j in range(2,iterate):
		print '\nj = ', j, '\nt = ', EI1
		ode_f = lambda t,y: weight_odes(t,y,EI1,PAL(t))
		example1 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example1.set_initial_value(init_FL, 0.) 
		y1 = 2.2*sum(example1.integrate(T) )
		
		if abs(y1-beta)<1e-8: 
			print '\n--Solution computed successfully--'
			break
		# Update
		ode_f = lambda t,y: weight_odes(t,y,EI0,PAL(t))
		example0 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example0.set_initial_value(init_FL, 0.) 
		y0 = 2.2*sum(example0.integrate(T))
		
		EI2 = EI1 - (y1 - beta)*(EI1-EI0)/(y1- y0)
		EI0 = EI1
		EI1 = EI2
		
	# Here we plot the solution 
	ode_f = lambda t,y: weight_odes(t,y,EI1,PAL(t))
	example = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol)
	example.set_initial_value(init_FL, 0.) 
	
	X = np.linspace(0.,T,801)
	Y = np.zeros((len(X),dim))
	Y[0,:] = init_FL
	
	for j in range(1,len(X)): Y[j,:] = example.integrate(X[j]) 
	
	print '\n|y(T) - Target Weight| = ', np.abs(beta-2.2*sum(Y[-1,:]) ),'\n'
	plt.plot(X,2.2*(Y[:,0]+Y[:,1]),'-k',linewidth=2.0)
	plt.axis([-1,T+1,0,170])
	plt.show(); plt.clf()
	return 

 def addMethods(self):
     # override to add _solver function
     ODEsystem.addMethods(self, usePsyco=False)
     if self.haveJacobian():
         self._solver = ode(getattr(self,self.funcspec.spec[1]),
                         getattr(self,self.funcspec.auxfns["Jacobian"][1]))
         self._funcreg['_solver'] = ('self',
              'ode(getattr(self,self.funcspec.spec[1]),' \
                + 'getattr(self,self.funcspec.auxfns["Jacobian"][1]))')
     else:
         self._solver = ode(getattr(self,self.funcspec.spec[1]))
         self._funcreg['_solver'] = ('self', 'ode(getattr(self,' \
                                               + 'self.funcspec.spec[1]))')

def DimensionalForms(BoundConditions, NDParameters, X0=param.X0, d=param.d, rICp=param.rICp, \
                     rhoH=param.rhoH, alpha=param.alpha, rhoD=param.rhoD, g=param.g, hminus=param.hmin, \
                     hmaxi=param.hmax, dh=param.dt, geom="cart", AbsTole=param.AbsTol, RelTole=param.AbsTol):
    """ Integration RK (de type Dormand Prince 4-5) of the equations system for an horizontal F-layer.
    resolutionSystem(y0, t0, tmax, dt, *arg)
    Default parameters for the others arguments are :
    PeT=1.e6,
    MX=0.17,
    MP=0.016,
    Veff=1.e7,
    gamma=7.2,
    Atol=1.e-10,
    Rtol=1.e-12
    The result is given in fully dimensional form. 
    """

    if geom == "cart":
        sol = ode(systemdiff.cart).set_integrator(
            'dopri5', atol=param.AbsTol, rtol=param.AbsTol, nsteps=1e7)
    elif geom == "spher":
        sol = ode(systemdiff.spher).set_integrator(
            'dopri5', atol=param.AbsTol, rtol=param.AbsTol, nsteps=1e7)

    sol.set_initial_value(BoundConditions, param.hmin)
    sol.set_f_params(NDParameters)  # (PeT, Mx, Mp, Veff, gamma)
    X, dXdh, phiVs,  h = BoundConditions[0] * np.ones(1), BoundConditions[1] * np.ones(1), \
      BoundConditions[2] * np.ones(1), param.hmin * np.ones(1)

    if param.dt > 0:
        while sol.successful() and sol.t < param.hmax:
            sol.integrate(sol.t + param.dt)
        # print("%g %g" % (r.t, r.y))
            X = np.append(X, sol.y[0])
            dXdh = np.append(dXdh, sol.y[1])
            phiVs = np.append(phiVs, sol.y[2])
            h = np.append(h, sol.t)
    else:
        while sol.successful() and sol.t > param.hmax:
            sol.integrate(sol.t + param.dt)
        # print("%g %g" % (r.t, r.y))
            X = np.append(X, sol.y[0])
            dXdh = np.append(dXdh, sol.y[1])
            phiVs = np.append(phiVs, sol.y[2])
            h = np.append(h, sol.t)

    # back to fully dimensional:
    x = X0 * X
    z = h * d
    T = param.Tliq0 * (1 - NDParameters['Mp'] * h - X * NDParameters['Mx'])

    return z, x, phiVs, T

def run():
    data_x = []
    data_y = []

    dt = 0.1
    tmax = 10
    state = None

    m = model.model_push()
    ir = sciint.ode(m.step).set_integrator('vode', atol=0.01)
    ir.set_initial_value(m.init, 0)

    while state is None or state[2] < m.x_h and ir.t < tmax:
        state = ir.integrate(ir.t + dt/1000)
        data_x.append(ir.t)
        data_y.append(m.M)

    if ir.t >= tmax:
        print('front failed')
        plt.ylim(0, m.Mmax + 10)
        plt.plot(data_x, data_y)
        plt.xlabel('Время t [с]')
        plt.ylabel('Момент М(t) [Н м]')
        plt.show()
        return None

    t1 = ir.t
    m = model.model_torque()
    ir = sciint.ode(m.step).set_integrator('vode', atol=0.01)
    m.init = [
        state[0], 0,
        state[2] - m.L * cos(asin(((state[3] - m.y_f)/m.L))), m.r_k + m.y_f
    ]
    state2 = state[:]
    ir.set_initial_value(m.init, 0)
    state = None
    while state is None or state[2] < m.x_h and ir.t + t1 < tmax:
        state = ir.integrate(ir.t + dt/1000)
        data_x.append(ir.t + t1)
        data_y.append(m.M)

    if ir.t + t1 >= tmax:
        print('rear failed')
        return None
    print('all pass')
    plt.ylim(0, m.Mmax + 10)
    plt.plot(data_x, data_y)
    plt.xlabel('Время t [с]')
    plt.ylabel('Момент М(t) [Н м]')
    plt.show()
    return (m.init, state2)

    def test_concurrent_fail(self):
        for sol in ('vode', 'zvode'):
            f = lambda t, y: 1.0

            r = ode(f).set_integrator(sol)
            r.set_initial_value(0, 0)

            r2 = ode(f).set_integrator(sol)
            r2.set_initial_value(0, 0)

            r.integrate(r.t + 0.1)
            r2.integrate(r2.t + 0.1)

            assert_raises(RuntimeError, r.integrate, r.t + 0.1)

def Cannon_Shooting(t0,t1,za=0.,va=45,nu=.0003): 
	a, b = 0., 195
	beta =  0.
	
	# t0,t1 = np.pi/6., np.pi/7
	dim, iterate = 3,40
	reltol, abstol = 1e-9,1e-8
	
	# Initial_Conditions = np.array([z=za,v=va,phi=t])
	g = 9.8067
	def ode_f(x,y): 
		# y = [z,v,phi]
		return np.array([np.tan(y[2]), -(g*np.sin(y[2]) + nu*y[1]**2.)/(y[1]*np.cos(y[2])), 
		-g/y[1]**2.])
	
	
	print '\nj = 1'
	print 't = ', t0
	for j in range(2,iterate):
		print '\nj = ', j
		print 't = ', t1
		example1 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example1.set_initial_value(np.array([za,va,t1]),a) 
		y1 = example1.integrate(b)[0]
		
		if abs(y1-beta)<1e-8: 
			print '\n--Solution y computed successfully--'
			break
		# Update
		example0 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example0.set_initial_value(np.array([za,va,t0]),a) 
		y0 = example0.integrate(b)[0]
		
		t2 = t1 - (y1 - beta)*(t1-t0)/(y1- y0)
		t0 = t1
		t1 = t2
		
	# Here we plot the solution 
	example = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol)
	example.set_initial_value(np.array([za,va,t1]),a) 
	X = np.linspace(a,b,801)
	Y = np.zeros((len(X),dim))
	Y[0,:] = np.array([0.,0.5,t1])
	
	for j in range(1,len(X)): Y[j,:] = example.integrate(X[j]) 
	
	print '\n|y(b) - beta| = ', np.abs(beta-Y[-1,0]),'\n'
	# plt.plot(X,Y[:,0],'-k')
	return X,Y