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
python类solve()的实例源码
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 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))
