# preface

In the first two articles, we learned the definition and use of SymPy's input and output, basic symbols, expressions, and functions, as well as the simplification, expansion, and merging of expressions.

This article introduces SymPy equation solving, including: linear/nonlinear equation solving, linear equation solving and nonlinear equation solving, and the solution results are divided into symbolic solutions and numerical solutions. 1 . Solving equations consists of two main functions: solve() and solveset().

In addition, by the way, learn how to solve summation formula, multiplication formula, function limit and sequence limit.

# solve equation(s)

There are two main functions that can find the root of the equation: solve() and solveset(), which can solve general equations, such as

from sympy.abc import x, y
from sympy import solve, solveset

solve(x**2 - y, x, dict=True)
# [{x: -sqrt(y)}, {x: sqrt(y)}]
solveset(x**2 - y, x)
# {-sqrt(y), sqrt(y)}


As for the difference between the two, the official document recommends that we choose according to the following requirements

1. solve()

1. The return value of the solve() function can be relatively easily replaced by subs() to participate in other operations
2. wish to explicitly obtain a symbolic expression that satisfies the equation
2. solveset()

1. A more accurate mathematical solution can be obtained (by setting mathematical sets)
2. Can show all solutions, even infinitely many solutions
3. The interval of the solution can be freely restricted

## solve function

solve(f, *symbols, **flags) can solve polynomials, transcendental functions, piecewise functions, linear equations, etc., for example:

x     = sympy.symbols('x')
expr1 = sympy.solve(x - 1, x)                   # find x that satisfies the equation x - 1 = 0
print(expr1)
# [1]

x, y = sympy.symbols('x y')
expr2 = sympy.solve([2*x-y-3,3*x+y-7],[x,y])    # system of linear equations
print(expr2)
# {x: 2, y: 1}


Symbolic solutions can also be obtained, such as y ( x ) = x + 1 x − 1 y(x) = \dfrac{x + 1}{x - 1} y(x)=x−1x+1​, find x ( y ) x(y) x(y)

x, y = symbols('x y')
expr = Eq((x + 1) / (x - 1), y)
print(solve(expr, x))                           # Feed Eq directly to the solve function
# [(y + 1)/(y - 1)]


Of course *symbols parameters can also use expressions or functions

solve(x + 2 + sqrt(3), x + 2)
# [-sqrt(3)]
solve(f(x) - x, f(x))
# [x]
solve(f(x).diff(x) - f(x) - x, f(x).diff(x))
#  [x + f(x)]

# Use indexable notation A[i]
from sympy import Indexed, IndexedBase, Tuple, sqrt
A = IndexedBase('A')
eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1)
solve(eqs, eqs.atoms(Indexed))
# {A[1]: 1, A[2]: 2}


The *symbols parameter can also be omitted, at which point SymPy automatically judges which symbols are to be solved

solve(x - 3)
# [3]
solve(x**2 - y**2)
# [{x: -y}, {x: y}]
solve(z**2*x**2 - z**2*y**2)
# [{x: -y}, {x: y}, {z: 0}]
solve(z**2*x - z**2*y**2)
# [{x: y**2}, {z: 0}]


If you don't want to get a complete and explicit solution, you can use the parameter implicit=True, such as solving x + e x = 0 x+e^x=0 x+ex=0, the solution can be obtained x = − e x x = -e^x x=−ex

solve(x + exp(x), x)
# [-LambertW(1)]
solve(x + exp(x), x, implicit=True)
# [-exp(x)]


Find the coefficient a , b a, b a,b such that the equation ( a + b ) x − b + 2 = 0 (a+b)x - b + 2 = 0 (a+b)x−b+2=0 is always established

solve((a + b)*x - b + 2, a, b)
# {a: -2, b: 2}


By default, SymPy will solve in the entire field of complex numbers. If you limit the solution to the field of real numbers, you can refer to the following method

from sympy import Symbol, solve
x = Symbol('x', real=True)
solve(x**4 - 256, x, dict=True)
[{x: -4}, {x: 4}]						# If real=True is not restricted, output {-4, 4, -4*I, 4*I}


Equations with Inequality Conditions

solve([x < 3, x - 2])
# Eq(x, 2)
solve([x > 3, x - 2])
# False


Some variable parameters **flags and default values:

• dict=False: returns a list of mappings formatted as solutions of equations
• set=False: the return format is a tuple composed of a symbol list and a solution
• check=True: Check the solution and exclude the solution that makes the denominator of the fraction zero
• simplify=True: Whether to simplify the solution, setting it to False can shorten the running time
• force=False: Whether to force simplification regardless of the sign's positive or negative

## solveset function

The function solveset(f, symbol=None, domain=Complexes) can solve equations or inequalities in the specified domain. If it is solved in the domain of real numbers, solveset_real can be used, and solveset_comples can be used in the domain of complex numbers.

x = Symbol('x')     # When x = Symbol('x', real=True), the solution results are the same (the solution domain must be specified in the solveset function)
pprint(solveset(exp(x) - 1, x), use_unicode=False)
# {2*n*I*pi | n in Integers}

R = S.Reals
x = Symbol('x')
solveset(exp(x) - 1, x, R)
# {0}

solveset_real(exp(x) - 1, x)
# {0}

# Solving Inequalities
solveset(exp(x) > 1, x, R)
# Interval.open(0, oo)


When solving a single equation, the solveset function returns a solution set (solveset), finite set (FiniteSet), interval (Interval) or imaginary number set (ImageSet). When the solution does not exist, it returns an EmpySet empty set. If there is no limit in the specified symbol Find understanding outside, return ConditionSet.

i = Symbol('i', imaginary=True)
solveset(Eq(i**2 + i*sin(i), 1), i, domain=S.Reals)
# ConditionSet(_R, Eq(_R**2 + _R*sin(_R) - 1, 0), Reals)


The solveset function does not specifically give multiple roots, you need to use the roots function to get the multiple roots of the polynomial

solveset(x**3 - 6*x**2 + 9*x, x)
# {0, 3}
roots(x**3 - 6*x**2 + 9*x, x)
# {0: 1, 3: 2}


Inequalities are only solved in the domain of real numbers. If an inequality is given in the domain of complex numbers, the error NotImplementedError will be returned

solveset(exp(x) > 1, x, R)
# Interval.open(0, oo)


Solve the inverse function invert_real(f_x, y, x) (the complex domain can use the function invert_complex(f_x, y, x, domain=Complexes))

invert_real(exp(x), y, x)
# (x, Intersection({log(y)}, Reals))

# When does exp(x) == 1?
invert_real(exp(x), 1, x)
# (x, {0})


Solve equations in a specified domain

from sympy import Interval, pi, sin, solveset
from sympy.abc import x
solveset(sin(x), x, Interval(-pi, pi))
# {0, -pi, pi}


Linear equations can be solved exclusively with the function linsolve(system, *symbols)

linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))
# {(-y - 1, y, 2)}

# augmented matrix form
linsolve(Matrix(([1, 1, 1, 1], [1, 1, 2, 3])), (x, y, z))
# {(-y - 1, y, 2)}

# Ax=b form
M = Matrix(((1, 1, 1, 1), (1, 1, 2, 3)))
system = A, b = M[:, :-1], M[:, -1]
linsolve(system, x, y, z)
# {(-y - 1, y, 2)}


Nonlinear equations can be solved with a special function nonlinsove (for trigonometric functions, it is best to use the solve function)

from sympy import sqrt
system = [x**2 - 2*y**2 -2, x*y - 2]
vars = [x, y]
nonlinsolve(system, vars)
# {(-2, -1), (2, 1), (-√2⋅ⅈ, √2⋅ⅈ), (√2⋅ⅈ, -√2⋅ⅈ)}


# to sum ∑ \sum ∑

sympy.summation(f, *symbols, **kwargs) , usage:

                            b
____
\
summation(f, (i, a, b)) =  )    f
/___,
i = a


Go straight to the chestnuts: please ∑ n = 1 100 3 n \sum_{n=1}^{100} 3n ∑n=1100​3n

n     = sympy.Symbol('n', integer=True)
expr1 = sympy.summation(3 * n,(n,1,100))
print(expr1)
# 15150


begging ∑ i = 1 n 2 i − 1 \sum_{i=1}^n 2i-1 ∑i=1n​2i−1

from sympy import summation, oo, symbols, log
i, n, m = symbols('i n m', integer=True)
summation(2*i - 1, (i, 1, n))
# n**2


begging ∑ i = 0 ∞ ( 1 2 ) i \sum_{i=0}^\infty \left(\frac{1}{2}\right)^i ∑i=0∞​(21​)i

summation(1/2**i, (i, 0, oo))
# 2


begging ∑ n = 0 m ∑ i = 0 n i , i = 0 , . . . , n , n = 0 , . . . , m \sum_{n=0}^m \sum_{i=0}^n i, \quad i=0, ..., n,\quad n=0, ..., m ∑n=0m​∑i=0n​i,i=0,...,n,n=0,...,m

summation(i, (i, 0, n), (n, 0, m))
# m**3/6 + m**2/2 + m/3


# Multiply ∏ \prod ∏

Multiplication function product(*args, **kwargs), usage:

                             b
_____
product(f(n), (i, a, b)) = |   | f(n)
|   |
i = a


directly on the chestnuts, ∏ i = 1 k i \prod_{i=1}^k i ∏i=1k​i:

from sympy import product, symbols
i, n, m, k = symbols('i n m k', integer=True)

product(i, (i, 1, k))
# factorial(k)


begging ∏ i = 1 k m \prod_{i=1}^k m ∏i=1k​m

product(m, (i, 1, k))
# m**k


begging ∏ k = 1 n ∏ i = 1 k i \prod_{k=1}^n \prod_{i=1}^k i ∏k=1n​∏i=1k​i

product(i, (i, 1, k), (k, 1, n))
# Product(factorial(k), (k, 1, n))


# find function limit

The function is limit(e, z, z0, dir='+'), which returns the limit of the e(z) function at z0. Approach direction dir has +right approach, -left approach, +-left and right approach. chestnut

from sympy import limit, sin, oo
from sympy.abc import x
limit(sin(x)/x, x, 0)
# 1
limit(1/x, x, 0) # default dir='+'
# oo
limit(1/x, x, 0, dir="-")
# -oo
limit(1/x, x, 0, dir='+-')
# zoo
limit(1/x, x, oo)
# 0


# find sequence limit

The function is limit_seq(expr, n=None, trials=5), which returns the limit of the expression expr when n tends to infinity. trials is the number of recursions. If None is returned, you can try to increase the number of recursions to get the result. chestnut:

from sympy import limit_seq, Sum, binomial
from sympy.abc import n, k, m
limit_seq((5*n**3 + 3*n**2 + 4) / (3*n**3 + 4*n - 5), n)
# 5/3
limit_seq(binomial(2*n, n) / Sum(binomial(2*k, k), (k, 1, n)), n)        # binomial binomial
# 3/4
limit_seq(Sum(k**2 * Sum(2**m/m, (m, 1, k)), (k, 1, n)) / (2**n*n), n)   # Sum (finite) summation
# 4
`
1. Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 https://doi.org/10.7717/peerj-cs.103 ↩︎

Posted by tauchai83 on Wed, 18 Jan 2023 01:52:12 +0530