| Hosted by CoCalc | Download
Kernel: SageMath 8.1

Solving equations with SageMath

%display latex

Exact solutions : solve

Let us start with a simple example: solving x2x1=0x^2-x-1=0 for xx:

solve(x^2-x-1 == 0, x)
[x=125+12,x=125+12]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = -\frac{1}{2} \, \sqrt{5} + \frac{1}{2}, x = \frac{1}{2} \, \sqrt{5} + \frac{1}{2}\right]

Note that the equation is formed with the operator ==.

The equation itself can be stored in a Python variable:

eq = x^2-x-1 == 0
eq
x2x1=0\renewcommand{\Bold}[1]{\mathbf{#1}}x^{2} - x - 1 = 0
eq.lhs()
x2x1\renewcommand{\Bold}[1]{\mathbf{#1}}x^{2} - x - 1
eq.rhs()
0\renewcommand{\Bold}[1]{\mathbf{#1}}0
solve(eq, x)
[x=125+12,x=125+12]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = -\frac{1}{2} \, \sqrt{5} + \frac{1}{2}, x = \frac{1}{2} \, \sqrt{5} + \frac{1}{2}\right]

The solutions are returned as a list:

sol = solve(eq, x) sol
[x=125+12,x=125+12]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = -\frac{1}{2} \, \sqrt{5} + \frac{1}{2}, x = \frac{1}{2} \, \sqrt{5} + \frac{1}{2}\right]
sol[1]
x=125+12\renewcommand{\Bold}[1]{\mathbf{#1}}x = \frac{1}{2} \, \sqrt{5} + \frac{1}{2}

Each element of the solution list is itself an equation, albeit a trivial one. This can be seen by using the print function instead of the default LaTeX display:

print(sol[1])
x == 1/2*sqrt(5) + 1/2

The access to the value of the solution is via the rhs() method:

sol[1].rhs()
125+12\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, \sqrt{5} + \frac{1}{2}

A numerical approximation, n being a shortcut for numerical_approx:

n(sol[1].rhs())
1.61803398874989\renewcommand{\Bold}[1]{\mathbf{#1}}1.61803398874989
n(sol[1].rhs(), digits=50)
1.6180339887498948482045868343656381177203091798058\renewcommand{\Bold}[1]{\mathbf{#1}}1.6180339887498948482045868343656381177203091798058

A new equation involving the above solution:

sol[1].rhs() == golden_ratio()
125+12=ϕ\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, \sqrt{5} + \frac{1}{2} = \phi

Asking Sage whether this equation always holds:

bool(_)
True\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}

A system with various unknowns

Since x is the only predefined symbolic variable in Sage, we need to declare the other symbolic variables denoting the unknowns, here y:

y = var('y')

Then we may form the system:

eq1 = x^2 + x*y + 2 == 0 eq2 = y^2 == x*(x+y) (eq1, eq2)
(x2+xy+2=0,y2=(x+y)x)\renewcommand{\Bold}[1]{\mathbf{#1}}\left(x^{2} + x y + 2 = 0, y^{2} = {\left(x + y\right)} x\right)
sol = solve((eq1, eq2), x, y) sol
[[x=12i10+12i2,y=i2],[x=12i10+12i2,y=i2],[x=12i1012i2,y=i2],[x=12i1012i2,y=i2]]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[\left[x = -\frac{1}{2} i \, \sqrt{10} + \frac{1}{2} i \, \sqrt{2}, y = -i \, \sqrt{2}\right], \left[x = \frac{1}{2} i \, \sqrt{10} + \frac{1}{2} i \, \sqrt{2}, y = -i \, \sqrt{2}\right], \left[x = -\frac{1}{2} i \, \sqrt{10} - \frac{1}{2} i \, \sqrt{2}, y = i \, \sqrt{2}\right], \left[x = \frac{1}{2} i \, \sqrt{10} - \frac{1}{2} i \, \sqrt{2}, y = i \, \sqrt{2}\right]\right]

Here again the solutions are returned as a list; but now, each item of the list is itself a list, containing the values for x and y:

sol[0]
[x=12i10+12i2,y=i2]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = -\frac{1}{2} i \, \sqrt{10} + \frac{1}{2} i \, \sqrt{2}, y = -i \, \sqrt{2}\right]
sol[0][1]
y=i2\renewcommand{\Bold}[1]{\mathbf{#1}}y = -i \, \sqrt{2}
sol[0][1].rhs()
i2\renewcommand{\Bold}[1]{\mathbf{#1}}-i \, \sqrt{2}
n(sol[0][1].rhs())
1.41421356237310i\renewcommand{\Bold}[1]{\mathbf{#1}}-1.41421356237310i

Approximate solution: find_root

Let us consider the following transcendental equation:

eq = x*sin(x) == 4 eq
xsin(x)=4\renewcommand{\Bold}[1]{\mathbf{#1}}x \sin\left(x\right) = 4

solve returns a non-trivial equation equivalent to its input, meaning it could not find a solution:

solve(eq, x)
[x=4sin(x)]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = \frac{4}{\sin\left(x\right)}\right]

We use then find_root to get an approximate solution in a given interval, here [0,10][0, 10]:

find_root(eq, 0, 10)
8.96212620082\renewcommand{\Bold}[1]{\mathbf{#1}}8.96212620082

Actually find_root always return a single solution, even if there are more than one in the prescribed interval. In the current case, there is a second solution, as we can see graphically:

plot(x*sin(x)-4, (x, 0, 10))
Image in a Jupyter notebook

We can get it by forcing the search in the interval [6,8][6,8]:

find_root(eq, 6, 8)
6.90141260902\renewcommand{\Bold}[1]{\mathbf{#1}}6.90141260902
find_root(eq, 6, 8, full_output=True)
(6.90141260902,xxxxxxconverged:xTruexxxxxxxxxxxflag:x'converged'xfunction_calls:x10xxxxxiterations:x9xxxxxxxxxxxroot:x6.901412609021304)\renewcommand{\Bold}[1]{\mathbf{#1}}\left(6.90141260902, \begin{array}{l} \phantom{\verb!xxxxxx!}\verb|converged:|\phantom{\verb!x!}\verb|True|\\ \phantom{\verb!xxxxxxxxxxx!}\verb|flag:|\phantom{\verb!x!}\verb|'converged'|\\ \phantom{\verb!x!}\verb|function_calls:|\phantom{\verb!x!}\verb|10|\\ \phantom{\verb!xxxxx!}\verb|iterations:|\phantom{\verb!x!}\verb|9|\\ \phantom{\verb!xxxxxxxxxxx!}\verb|root:|\phantom{\verb!x!}\verb|6.901412609021304| \end{array}\right)

Approximate solution to a system: mpmath.findroot

eq1 = x*sin(x) - y == 0 eq2 = y*cos(y) - x - 1 == 0 for eq in [eq1, eq2]: show(eq)
xsin(x)y=0\renewcommand{\Bold}[1]{\mathbf{#1}}x \sin\left(x\right) - y = 0
ycos(y)x1=0\renewcommand{\Bold}[1]{\mathbf{#1}}y \cos\left(y\right) - x - 1 = 0
f1(x,y) = eq1.lhs() f1
(x,y)  xsin(x)y\renewcommand{\Bold}[1]{\mathbf{#1}}\left( x, y \right) \ {\mapsto} \ x \sin\left(x\right) - y
f2(x,y) = eq2.lhs() f2
(x,y)  ycos(y)x1\renewcommand{\Bold}[1]{\mathbf{#1}}\left( x, y \right) \ {\mapsto} \ y \cos\left(y\right) - x - 1

For solving a system, we use the mpmath toolbox:

import mpmath
f = [lambda a,b: f1(RR(a), RR(b)), lambda a,b: f2(RR(a), RR(b))] sol = mpmath.findroot(f, (0, 0)) sol
[-0.642871006593746][x0.385398467857987]\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{l} \verb|[-0.642871006593746]|\\ \verb|[|\phantom{\verb!x!}\verb|0.385398467857987]| \end{array}
sol.tolist()
[[0.642871006593746],[0.385398467857987]]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[\left[-0.642871006593746\right], \left[0.385398467857987\right]\right]
x1 = RR(sol.tolist()[0][0]) y1 = RR(sol.tolist()[1][0]) x1, y1
(0.642871006593746,0.385398467857987)\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-0.642871006593746, 0.385398467857987\right)
f1(x1,y1)
5.55111512312578×1016\renewcommand{\Bold}[1]{\mathbf{#1}}-5.55111512312578 \times 10^{-16}
f2(x1,y1)
1.11022302462516×1016\renewcommand{\Bold}[1]{\mathbf{#1}}1.11022302462516 \times 10^{-16}
sol2 = minimize(eq1.lhs()^2 + eq2.lhs()^2, (-0.6,0.4)) sol2
(0.642870824741,0.385398567605)\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-0.642870824741,\,0.385398567605\right)
x2, y2 = sol2[0], sol2[1] x2, y2
(0.642870824741,0.385398567605)\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-0.642870824741, 0.385398567605\right)
f1(x2,y2), f2(x2,y2)
(3.02336881641×1007,1.0387405891×1007)\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-3.02336881641 \times 10^{-07}, -1.0387405891 \times 10^{-07}\right)
x1 - x2, y1 - y2
(1.81852617853195×107,9.97466205743258×108)\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-1.81852617853195 \times 10^{-7}, -9.97466205743258 \times 10^{-8}\right)