我们从Python开源项目中，提取了以下28个代码示例，用于说明如何使用sympy.solve()。
def _nthroot_solve(p, n, prec): """ helper function for ``nthroot`` It denests ``p**Rational(1, n)`` using its minimal polynomial """ from sympy.polys.numberfields import _minimal_polynomial_sq from sympy.solvers import solve while n % 2 == 0: p = sqrtdenest(sqrt(p)) n = n // 2 if n == 1: return p pn = p**Rational(1, n) x = Symbol('x') f = _minimal_polynomial_sq(p, n, x) if f is None: return None sols = solve(f, x) for sol in sols: if abs(sol - pn).n() < 1./10**prec: sol = sqrtdenest(sol) if _mexpand(sol**n) == p: return sol
def testAlgebraInverse(self): dataset_objects = algorithmic_math.math_dataset_init(26) counter = 0 for d in algorithmic_math.algebra_inverse(26, 0, 3, 10): counter += 1 decoded_input = dataset_objects.int_decoder(d["inputs"]) solve_var, expression = decoded_input.split(":") lhs, rhs = expression.split("=") # Solve for the solve-var. result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var) target_expression = dataset_objects.int_decoder(d["targets"]) # Check that the target and sympy's solutions are equivalent. self.assertEqual( 0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression)) self.assertEqual(counter, 10)
def execute_devito(ui, spacing=0.01, a=0.5, timesteps=500): """Execute diffusion stencil using the devito Operator API.""" nx, ny = ui.shape dx2, dy2 = spacing**2, spacing**2 dt = dx2 * dy2 / (2 * a * (dx2 + dy2)) # Allocate the grid and set initial condition # Note: This should be made simpler through the use of defaults grid = Grid(shape=(nx, ny)) u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2) u.data[0, :] = ui[:] # Derive the stencil according to devito conventions eqn = Eq(u.dt, a * (u.dx2 + u.dy2)) stencil = sympy.solve(eqn, u.forward)[0] op = Operator(Eq(u.forward, stencil)) # Execute the generated Devito stencil operator tstart = time.time() op.apply(u=u, t=timesteps, dt=dt) runtime = time.time() - tstart log("Devito: Diffusion with dx=%0.4f, dy=%0.4f, executed %d timesteps in %f seconds" % (spacing, spacing, timesteps, runtime)) return u.data[1, :], runtime
def run_simulation(save=False, dx=0.01, dy=0.01, a=0.5, timesteps=100): nx, ny = int(1 / dx), int(1 / dy) dx2, dy2 = dx**2, dy**2 dt = dx2 * dy2 / (2 * a * (dx2 + dy2)) grid = Grid(shape=(nx, ny)) u = TimeFunction( name='u', grid=grid, save=timesteps, initializer=initializer, time_order=1, space_order=2 ) eqn = Eq(u.dt, a * (u.dx2 + u.dy2)) stencil = solve(eqn, u.forward)[0] op = Operator(Eq(u.forward, stencil), time_axis=Forward) op.apply(time=timesteps, dt=dt) if save: return u.data[timesteps - 1, :] else: return u.data[(timesteps+1) % 2, :]
def _new_operator3(shape, time_order, **kwargs): grid = Grid(shape=shape) spacing = 0.1 a = 0.5 c = 0.5 dx2, dy2 = spacing**2, spacing**2 dt = dx2 * dy2 / (2 * a * (dx2 + dy2)) # Allocate the grid and set initial condition # Note: This should be made simpler through the use of defaults u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2) u.data[0, :] = np.arange(reduce(mul, shape), dtype=np.int32).reshape(shape) # Derive the stencil according to devito conventions eqn = Eq(u.dt, a * (u.dx2 + u.dy2) - c * (u.dxl + u.dyl)) stencil = solve(eqn, u.forward, rational=False)[0] op = Operator(Eq(u.forward, stencil), **kwargs) # Execute the generated Devito stencil operator op.apply(u=u, t=10, dt=dt) return u.data[1, :], op
def solveDAEs(self): for i, dae in enumerate(self): var, res = dae.solve() if len(res) > 0: t_var = self.__model.listOfVariables.getBySymbol(var) t_formula = MathFormula(self.__model) if isinstance(res[0], dict): t_formula.setInternalMathFormula(res[0].values()[0]) else: t_formula.setInternalMathFormula(res[0]) cfe = CFE(self.__model) cfe.new(t_var, t_formula) self.__model.listOfCFEs.append(cfe) self.__model.listOfVariables.changeVariableType(t_var, MathVariable.VAR_ASS) list.__delitem__(self, i) self.__model.listOfCFEs.developCFEs()
def eval_trigsubstitution(theta, func, rewritten, substep, integrand, symbol): func = func.subs(sympy.sec(theta), 1/sympy.cos(theta)) trig_function = list(func.find(TrigonometricFunction)) assert len(trig_function) == 1 trig_function = trig_function[0] relation = sympy.solve(symbol - func, trig_function) assert len(relation) == 1 numer, denom = sympy.fraction(relation[0]) if isinstance(trig_function, sympy.sin): opposite = numer hypotenuse = denom adjacent = sympy.sqrt(denom**2 - numer**2) inverse = sympy.asin(relation[0]) elif isinstance(trig_function, sympy.cos): adjacent = numer hypotenuse = denom opposite = sympy.sqrt(denom**2 - numer**2) inverse = sympy.acos(relation[0]) elif isinstance(trig_function, sympy.tan): opposite = numer adjacent = denom hypotenuse = sympy.sqrt(denom**2 + numer**2) inverse = sympy.atan(relation[0]) substitution = [ (sympy.sin(theta), opposite/hypotenuse), (sympy.cos(theta), adjacent/hypotenuse), (sympy.tan(theta), opposite/adjacent), (theta, inverse) ] return _manualintegrate(substep).subs(substitution).trigsimp()
def control_systems(request): ct_sys, ref = request.param Ac, Bc, Cc = ct_sys.data Dc = np.zeros((Cc.shape[0], 1)) Q = np.eye(Ac.shape[0]) R = np.eye(Bc.shape[1] if len(Bc.shape) > 1 else 1) Sc = linalg.solve_continuous_are(Ac, Bc.reshape(-1, 1), Q, R,) Kc = linalg.solve(R, Bc.T @ Sc).reshape(1, -1) ct_ctr = LTISystem(Kc) evals = np.sort(np.abs( linalg.eig(Ac, left=False, right=False, check_finite=False) )) dT = 1/(2*evals[-1]) Tsim = (8/np.min(evals[~np.isclose(evals, 0)]) if np.sum(np.isclose(evals[np.nonzero(evals)], 0)) > 0 else 8 ) dt_data = signal.cont2discrete((Ac, Bc.reshape(-1, 1), Cc, Dc), dT) Ad, Bd, Cd, Dd = dt_data[:-1] Sd = linalg.solve_discrete_are(Ad, Bd.reshape(-1, 1), Q, R,) Kd = linalg.solve(Bd.T @ Sd @ Bd + R, Bd.T @ Sd @ Ad) dt_sys = LTISystem(Ad, Bd, dt=dT) dt_sys.initial_condition = ct_sys.initial_condition dt_ctr = LTISystem(Kd, dt=dT) yield ct_sys, ct_ctr, dt_sys, dt_ctr, ref, Tsim
def test_events(): # use bouncing ball to test events work # simulate in block diagram int_opts = block_diagram.DEFAULT_INTEGRATOR_OPTIONS.copy() int_opts['rtol'] = 1E-12 int_opts['atol'] = 1E-15 int_opts['nsteps'] = 1000 int_opts['max_step'] = 2**-3 x = x1, x2 = Array(dynamicsymbols('x_1:3')) mu, g = sp.symbols('mu g') constants = {mu: 0.8, g: 9.81} ic = np.r_[10, 15] sys = SwitchedSystem( x1, Array([0]), state_equations=r_[x2, -g], state_update_equation=r_[sp.Abs(x1), -mu*x2], state=x, constants_values=constants, initial_condition=ic ) bd = BlockDiagram(sys) res = bd.simulate(5, integrator_options=int_opts) # compute actual impact time tvar = dynamicsymbols._t impact_eq = (x2*tvar - g*tvar**2/2 + x1).subs( {x1: ic[0], x2: ic[1], g: 9.81} ) t_impact = sp.solve(impact_eq, tvar)[-1] # make sure simulation actually changes velocity sign around impact abs_diff_impact = np.abs(res.t - t_impact) impact_idx = np.where(abs_diff_impact == np.min(abs_diff_impact))[0] assert np.sign(res.x[impact_idx-1, 1]) != np.sign(res.x[impact_idx+1, 1])
def equilibrium_points(self, input_=None): return sp.solve(self.state_equation, self.state, dict=True)
def inverse_eqn(eqn): """ Inverse a symbolic expression with variable x Keeps x as the variable in the inverted expression > inverse('4*x/2') 'y = x/2' """ e = sympify('-x + ' + eqn.replace('x', 'y')) return str(solve(e, 'y')[0]).replace('x', 'float(x)')
def solow_residual(self): """ Symbolic expression for the Solow residual which is used as a measure of technology. :getter: Return the symbolic expression. :type: sym.Basic """ return sym.solve(Y - self.output, A)[0]
def subs(self, substitution_dict): """The subs() method allows for unknowns found from other ControlVolumes to be substituted into the equations for this system, reducing the total number of unknowns and the total degrees of freedom (hopefully until the degrees of freedom for the system reach zero and it becomes solvable). It requires a dictionary of variable:solution format from which it will substitute values into the ControlVolume equations_dict.""" # This list is created to store equations that are no longer useful and remove them # This occurs when an equation (generally an info equation) is used in another ControlVolume which implies that # all variables in that equation have been solved and it cannot provide any new relationships to the system remove_equation_list = [] # For each solved variable in substitution_dict for substitution, solution in substitution_dict.items(): # The substitution must meet 3 characteristics: it must exist in the current ControlVolume as a variable, # the variable in the ControlVolume must be unknown (a sympy Symbol, not a value), and the substitution must # be solved (it itself has zero unknowns in the form of sympy Symbols) if substitution in self.dict_of_variables and type(self.dict_of_variables[substitution]) == sp.Symbol and \ len(solution.atoms(sp.Symbol)) < 1: # If this is true, then the ControlVolume can remove the variable from it's dict_of_variables as it has # already been solved and doesn't need to be solved again. The total unknowns decreases by one self.dict_of_variables.pop(substitution) self.unknowns -= 1 # Each equation needs to substitute the unknown for it's solved solution using the sympy subs() method for key, equation in self.equations_dict.items(): self.equations_dict[key] = equation.subs(substitution, solution) # This if statement checks if the equation has become irrelevant (nothing to solve, just 0) # If the equation lacks unknowns, it will be removed from the equations_dict for the ControlVolume if len(self.equations_dict[key].atoms(sp.Symbol)) == 0: remove_equation_list.append(key) # This loop removes every equation that is no longer useful for key in remove_equation_list: self.equations_dict.pop(key) # After a substitution is done, the degrees of freedom have likely changed, so the ControlVolume will update it self.degrees_of_freedom_update()
def degrees_of_freedom_update(self): """This method calculates and updates the degrees of freedom for a system. If the system has as many equations as it has unknowns, then it as zero degrees of freedom and will be solved. Otherwise, it must wait for a substitution to occur to lower the unknowns in the system.""" # The degrees of freedom equation self.degrees_of_freedom = self.unknowns - len(self.equations_dict) + 1 # If the system lacks any degrees of freedom, then it is solvable, and the solve() method will run if self.degrees_of_freedom == 0: print('Solving control volume {}'.format(self)) self.solve()
def solve(self): """The solve method will build a matrix that sympy can solve with the sympy.solve() function. It will return the values in a dict which will then be used to store all solved unknowns to the dict_of_variables of the system.""" # A pre-allocation for the matrix used to solve the system matrix = [] # Each unknown must be put into a list so sympy can solve it unknowns_list = list(self.dict_of_variables.keys()) # Each equation (except for the 'Total') will be appended to the matrix. This is done to allow for the user # or the code (when this feature is added) to easily double check the variables for accuracy for key, equation in self.equations_dict.items(): if key != 'Total': matrix.append(equation) # sympy does it's thing and returns a dict in the form of {symbol: solution} solutions = sp.solve(matrix, unknowns_list, dict=True) # This loop updates the dict_of_variables with the newly solved values for each for solutions_set in solutions: # This is done because the solutions are given in a list containing a dictionary: [{}], which is weird for count in range(len(solutions_set)): # The newly solved variables can be used to solve other ControlVolumes self.dict_of_variables[unknowns_list[count]] = solutions_set[unknowns_list[count]]
def eval_trigsubstitution(theta, func, rewritten, substep, restriction, integrand, symbol): func = func.subs(sympy.sec(theta), 1/sympy.cos(theta)) trig_function = list(func.find(TrigonometricFunction)) assert len(trig_function) == 1 trig_function = trig_function[0] relation = sympy.solve(symbol - func, trig_function) assert len(relation) == 1 numer, denom = sympy.fraction(relation[0]) if isinstance(trig_function, sympy.sin): opposite = numer hypotenuse = denom adjacent = sympy.sqrt(denom**2 - numer**2) inverse = sympy.asin(relation[0]) elif isinstance(trig_function, sympy.cos): adjacent = numer hypotenuse = denom opposite = sympy.sqrt(denom**2 - numer**2) inverse = sympy.acos(relation[0]) elif isinstance(trig_function, sympy.tan): opposite = numer adjacent = denom hypotenuse = sympy.sqrt(denom**2 + numer**2) inverse = sympy.atan(relation[0]) substitution = [ (sympy.sin(theta), opposite/hypotenuse), (sympy.cos(theta), adjacent/hypotenuse), (sympy.tan(theta), opposite/adjacent), (theta, inverse) ] return sympy.Piecewise( (_manualintegrate(substep).subs(substitution).trigsimp(), restriction) )
def solveEquation(expr): expr = expr.replace('^', '**') members = expr.split('=') if len(members) != 2: raise BadFormatException('Bad number of equals.') from sympy.abc import * eq = sympy.Eq(*map(eval, members)) return [{repr(j): repr(k) for j, k in i.items()} for i in _ensureList(sympy.solve(eq, dict=1))]
def _pos_dependent(matrix, vector): """Check if vector can be written as a positive linear combination of the columns of matrix. Return the coeffients if possible, None otherwise.""" if len(matrix) > 0: coeffs = Matrix(list(map(lambda i: Symbol("x" + str(i)), range(len(matrix))))) sols = solve(Matrix(matrix).T * coeffs - Matrix(vector), coeffs, dict = True, particular = True) for sol in sols: if all([sol[s] >= 0 for s in sol]): return [sol[s] if s in sol else 0 for s in coeffs] return None
def execute_lambdify(ui, spacing=0.01, a=0.5, timesteps=500): """Execute diffusion stencil using vectorised numpy array accesses.""" nx, ny = ui.shape dx2, dy2 = spacing**2, spacing**2 dt = dx2 * dy2 / (2 * a * (dx2 + dy2)) u = np.concatenate((ui, np.zeros_like(ui))).reshape((2, nx, ny)) def diffusion_stencil(): """Create stencil and substitutions for the diffusion equation""" p = sympy.Function('p') x, y, t, h, s = sympy.symbols('x y t h s') dx2 = p(x, y, t).diff(x, x).as_finite_difference([x - h, x, x + h]) dy2 = p(x, y, t).diff(y, y).as_finite_difference([y - h, y, y + h]) dt = p(x, y, t).diff(t).as_finite_difference([t, t + s]) eqn = sympy.Eq(dt, a * (dx2 + dy2)) stencil = sympy.solve(eqn, p(x, y, t + s))[0] return stencil, (p(x, y, t), p(x + h, y, t), p(x - h, y, t), p(x, y + h, t), p(x, y - h, t), s, h) stencil, subs = diffusion_stencil() kernel = sympy.lambdify(subs, stencil, 'numpy') # Execute timestepping loop with alternating buffers tstart = time.time() for ti in range(timesteps): t0 = ti % 2 t1 = (ti + 1) % 2 u[t1, 1:-1, 1:-1] = kernel(u[t0, 1:-1, 1:-1], u[t0, 2:, 1:-1], u[t0, :-2, 1:-1], u[t0, 1:-1, 2:], u[t0, 1:-1, :-2], dt, spacing) runtime = time.time() - tstart log("Lambdify: Diffusion with dx=%0.4f, dy=%0.4f, executed %d timesteps in %f seconds" % (spacing, spacing, timesteps, runtime)) return u[ti % 2, :, :], runtime
def iso_stencil(field, time_order, m, s, damp, **kwargs): """ Stencil for the acoustic isotropic wave-equation: u.dt2 - H + damp*u.dt = 0 :param field: Symbolic TimeFunction object, solution to be computed :param time_order: time order :param m: square slowness :param s: symbol for the time-step :param damp: ABC dampening field (Function) :param kwargs: forwad/backward wave equation (sign of u.dt will change accordingly as well as the updated time-step (u.forwad or u.backward) :return: Stencil for the wave-equation """ # Creat a temporary symbol for H to avoid expensive sympy solve H = Symbol('H') # Define time sep to be updated next = field.forward if kwargs.get('forward', True) else field.backward # Define PDE eq = m * field.dt2 - H - kwargs.get('q', 0) # Add dampening field according to the propagation direction eq += damp * field.dt if kwargs.get('forward', True) else -damp * field.dt # Solve the symbolic equation for the field to be updated eq_time = solve(eq, next, rational=False, simplify=False)[0] # Get the spacial FD lap = laplacian(field, time_order, m, s) # return the Stencil with H replaced by its symbolic expression return [Eq(next, eq_time.subs({H: lap}))]
def solve(self): to_solve = [] for var in self.__definition.getDeveloppedInternalMathFormula().atoms(SympySymbol): variable = self.__model.listOfVariables.getBySymbol(var) if variable is not None and variable.isAlgebraic(): to_solve.append(var) return (to_solve[0], solve(self.__definition.getDeveloppedInternalMathFormula(), to_solve))
def idiff(eq, y, x, n=1): """Return ``dy/dx`` assuming that ``eq == 0``. Parameters ========== y : the dependent variable or a list of dependent variables (with y first) x : the variable that the derivative is being taken with respect to n : the order of the derivative (default is 1) Examples ======== >>> from sympy.abc import x, y, a >>> from sympy.geometry.util import idiff >>> circ = x**2 + y**2 - 4 >>> idiff(circ, y, x) -x/y >>> idiff(circ, y, x, 2).simplify() -(x**2 + y**2)/y**3 Here, ``a`` is assumed to be independent of ``x``: >>> idiff(x + a + y, y, x) -1 Now the x-dependence of ``a`` is made explicit by listing ``a`` after ``y`` in a list. >>> idiff(x + a + y, [y, a], x) -Derivative(a, x) - 1 See Also ======== sympy.core.function.Derivative: represents unevaluated derivatives sympy.core.function.diff: explicitly differentiates wrt symbols """ if is_sequence(y): dep = set(y) y = y[0] elif isinstance(y, Symbol): dep = set([y]) else: raise ValueError("expecting x-dependent symbol(s) but got: %s" % y) f = dict([(s, Function( s.name)(x)) for s in eq.free_symbols if s != x and s in dep]) dydx = Function(y.name)(x).diff(x) eq = eq.subs(f) derivs = {} for i in range(n): yp = solve(eq.diff(x), dydx)[0].subs(derivs) if i == n - 1: return yp.subs([(v, k) for k, v in f.items()]) derivs[dydx] = yp eq = dydx - yp dydx = dydx.diff(x)
def ratint_ratpart(f, g, x): """ Horowitz-Ostrogradsky algorithm. Given a field K and polynomials f and g in K[x], such that f and g are coprime and deg(f) < deg(g), returns fractions A and B in K(x), such that f/g = A' + B and B has square-free denominator. Examples ======== >>> from sympy.integrals.rationaltools import ratint_ratpart >>> from sympy.abc import x, y >>> from sympy import Poly >>> ratint_ratpart(Poly(1, x, domain='ZZ'), ... Poly(x + 1, x, domain='ZZ'), x) (0, 1/(x + 1)) >>> ratint_ratpart(Poly(1, x, domain='EX'), ... Poly(x**2 + y**2, x, domain='EX'), x) (0, 1/(x**2 + y**2)) >>> ratint_ratpart(Poly(36, x, domain='ZZ'), ... Poly(x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2, x, domain='ZZ'), x) ((12*x + 6)/(x**2 - 1), 12/(x**2 - x - 2)) See Also ======== ratint, ratint_logpart """ from sympy import solve f = Poly(f, x) g = Poly(g, x) u, v, _ = g.cofactors(g.diff()) n = u.degree() m = v.degree() A_coeffs = [ Dummy('a' + str(n - i)) for i in range(0, n) ] B_coeffs = [ Dummy('b' + str(m - i)) for i in range(0, m) ] C_coeffs = A_coeffs + B_coeffs A = Poly(A_coeffs, x, domain=ZZ[C_coeffs]) B = Poly(B_coeffs, x, domain=ZZ[C_coeffs]) H = f - A.diff()*v + A*(u.diff()*v).quo(u) - B*u result = solve(H.coeffs(), C_coeffs) A = A.as_expr().subs(result) B = B.as_expr().subs(result) rat_part = cancel(A/u.as_expr(), x) log_part = cancel(B/v.as_expr(), x) return rat_part, log_part
def intersection(self, o): """The intersection of the parabola and another geometrical entity `o`. Parameters ========== o : GeometryEntity, LinearEntity Returns ======= intersection : list of GeometryEntity objects Examples ======== >>> from sympy import Parabola, Point, Ellipse, Line, Segment >>> p1 = Point(0,0) >>> l1 = Line(Point(1, -2), Point(-1,-2)) >>> parabola1 = Parabola(p1, l1) >>> parabola1.intersection(Ellipse(Point(0, 0), 2, 5)) [Point2D(-2, 0), Point2D(2, 0)] >>> parabola1.intersection(Line(Point(-7, 3), Point(12, 3))) [Point2D(-4, 3), Point2D(4, 3)] >>> parabola1.intersection(Segment((-12, -65), (14, -68))) [] """ x, y = symbols('x y', real=True) parabola_eq = self.equation() if isinstance(o, Parabola): if o in self: return [o] else: return list(ordered([Point(i) for i in solve([parabola_eq, o.equation()], [x, y])])) elif isinstance(o, Point2D): if simplify(parabola_eq.subs(([(x, o._args[0]), (y, o._args[1])]))) == 0: return [o] else: return [] elif isinstance(o, (Segment2D, Ray2D)): result = solve([parabola_eq, Line2D(o.points[0], o.points[1]).equation()], [x, y]) return list(ordered([Point2D(i) for i in result if i in o])) elif isinstance(o, (Line2D, Ellipse)): return list(ordered([Point2D(i) for i in solve([parabola_eq, o.equation()], [x, y])])) elif isinstance(o, LinearEntity3D): raise TypeError('Entity must be two dimensional, not three dimensional') else: raise TypeError('Wrong type of argument were put')
def idiff(eq, y, x, n=1): """Return ``dy/dx`` assuming that ``eq == 0``. Parameters ========== y : the dependent variable or a list of dependent variables (with y first) x : the variable that the derivative is being taken with respect to n : the order of the derivative (default is 1) Examples ======== >>> from sympy.abc import x, y, a >>> from sympy.geometry.util import idiff >>> circ = x**2 + y**2 - 4 >>> idiff(circ, y, x) -x/y >>> idiff(circ, y, x, 2).simplify() -(x**2 + y**2)/y**3 Here, ``a`` is assumed to be independent of ``x``: >>> idiff(x + a + y, y, x) -1 Now the x-dependence of ``a`` is made explicit by listing ``a`` after ``y`` in a list. >>> idiff(x + a + y, [y, a], x) -Derivative(a, x) - 1 See Also ======== sympy.core.function.Derivative: represents unevaluated derivatives sympy.core.function.diff: explicitly differentiates wrt symbols """ if is_sequence(y): dep = set(y) y = y[0] elif isinstance(y, Symbol): dep = {y} else: raise ValueError("expecting x-dependent symbol(s) but got: %s" % y) f = dict([(s, Function( s.name)(x)) for s in eq.free_symbols if s != x and s in dep]) dydx = Function(y.name)(x).diff(x) eq = eq.subs(f) derivs = {} for i in range(n): yp = solve(eq.diff(x), dydx)[0].subs(derivs) if i == n - 1: return yp.subs([(v, k) for k, v in f.items()]) derivs[dydx] = yp eq = dydx - yp dydx = dydx.diff(x)
def _algebraic_equations(self, init_time, text, debug, queue): result_data = calc_result_data(text, True) result_data.calc_type = calc_type.ALGEBRAIC_EQUATIONS text = text_calculator.formula_to_py(result_data.formula_str) try: start_time = init_time text_line = text.split(text_calculator.EQUATION_VAR_FORMULA_SEPARATOR) if len(text_line) < 2: result_data.success = False result_data.calc_result = error.string_calculator.wrong_format_to_calc_equations() else: var_org = text_line[0] var_init_field = var_org.replace(u' ', u',') var_init_symbol = var_org formula_list = text_line[1:] if any((not formula.endswith(text_calculator.EQUATION_KEYWORD)) for formula in formula_list): result_data.success = False result_data.calc_result = error.string_calculator.wrong_format_to_calc_equations() else: formula_list_replaced = [eq.replace(text_calculator.EQUATION_KEYWORD, u'') for eq in formula_list] exec_py = '{} = sympy.symbols(\'{}\', real=True)'.format(var_init_field, var_init_symbol) exec_py += '\nresult = sympy.solve([{}], {})'.format(','.join(formula_list_replaced), var_init_field) start_time = init_time exec(exec_py) in globals(), locals() result_data.auto_record_time(start_time) result_data.success = True start_time = time.time() str_calc_result = str(result) result_data.latex = sympy.latex(result) result_data.auto_record_time(start_time) result_data.formula_str = '\n'.join(formula_list) result_data.calc_result = str_calc_result except Exception as ex: result_data.success = False result_data.calc_result = '{} - {}'.format(type(ex), ex.message) result_data.auto_record_time(start_time) queue.put(result_data)
def remove_by_cons(self, species, cons_law, debug = False): """Remove a species using a conservation law. First replace removed_species in the conservation law with their expression. Then use the conservation expression to write the species concentration as function of the remaining species. :Example: >>> from crnpy.crn import CRN, from_react_strings >>> net = from_react_strings(["E + S (k_1)<->(k1) C", "C ->(k2) E + P"]) >>> net.qss("C") >>> net.reactions (r0_r1: E + S ->(k1*k2/(k2 + k_1)) E + P,) >>> net.removed_species (('C', E*S*k1/(k2 + k_1)),) >>> cl = ConsLaw("E + C", "etot") >>> net.remove_by_cons("E", cl) >>> net.reactions (r0_r1: S ->(etot*k1*k2/(S*k1 + k2 + k_1)) P,) References: Tonello et al. (2016), On the elimination of intermediate species in chemical reaction networks. """ conservation = cons_law.expression if debug: print("Removed species: {}".format(self.removed_species)) print("Conservation: {}".format(conservation)) for variable, expr in self.removed_species: if debug: print("Replacing {} with {}".format(variable, expr)) conservation = conservation.subs(variable, expr) if debug: print("Found {}".format(conservation)) print # The next is quicker, but not always applicable #conservation = (conservation / sp.Symbol(species)).cancel() #exp = cons_law.constant / conservation exp = sp.solve(conservation - cons_law.constant, sp.Symbol(species))[0] # remove species self.remove_constant(species, expr = exp) if debug: print("Remove by Conservation: added to removed_species {}".format(self.removed_species))
