Sharedsolve_tour.ipynbOpen in CoCalc

# Solving equations with SageMath¶

%display latex


## Exact solutions : solve¶

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

solve(x^2-x-1 == 0, x)

$\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

$x^{2} - x - 1 = 0$
eq.lhs()

$x^{2} - x - 1$
eq.rhs()

$0$
solve(eq, x)

$\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

$\left[x = -\frac{1}{2} \, \sqrt{5} + \frac{1}{2}, x = \frac{1}{2} \, \sqrt{5} + \frac{1}{2}\right]$
sol[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()

$\frac{1}{2} \, \sqrt{5} + \frac{1}{2}$

A numerical approximation, n being a shortcut for numerical_approx:

n(sol[1].rhs())

$1.61803398874989$
n(sol[1].rhs(), digits=50)

$1.6180339887498948482045868343656381177203091798058$

A new equation involving the above solution:

sol[1].rhs() == golden_ratio()

$\frac{1}{2} \, \sqrt{5} + \frac{1}{2} = \phi$

Asking Sage whether this equation always holds:

bool(_)

$\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)

$\left(x^{2} + x y + 2 = 0, y^{2} = {\left(x + y\right)} x\right)$
sol = solve((eq1, eq2), x, y)
sol

$\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]

$\left[x = -\frac{1}{2} i \, \sqrt{10} + \frac{1}{2} i \, \sqrt{2}, y = -i \, \sqrt{2}\right]$
sol[0][1]

$y = -i \, \sqrt{2}$
sol[0][1].rhs()

$-i \, \sqrt{2}$
n(sol[0][1].rhs())

$-1.41421356237310i$

## Approximate solution: find_root¶

Let us consider the following transcendental equation:

eq = x*sin(x) == 4
eq

$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)

$\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]$:

find_root(eq, 0, 10)

$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))


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

find_root(eq, 6, 8)

$6.90141260902$
find_root(eq, 6, 8, full_output=True)

$\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)

$x \sin\left(x\right) - y = 0$
$y \cos\left(y\right) - x - 1 = 0$
f1(x,y) = eq1.lhs()
f1

$\left( x, y \right) \ {\mapsto} \ x \sin\left(x\right) - y$
f2(x,y) = eq2.lhs()
f2

$\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

$\begin{array}{l} \verb|[-0.642871006593746]|\\ \verb|[|\phantom{\verb!x!}\verb|0.385398467857987]| \end{array}$
sol.tolist()

$\left[\left[-0.642871006593746\right], \left[0.385398467857987\right]\right]$
x1 = RR(sol.tolist()[0][0])
y1 = RR(sol.tolist()[1][0])
x1, y1

$\left(-0.642871006593746, 0.385398467857987\right)$
f1(x1,y1)

$-5.55111512312578 \times 10^{-16}$
f2(x1,y1)

$1.11022302462516 \times 10^{-16}$
sol2 = minimize(eq1.lhs()^2 + eq2.lhs()^2, (-0.6,0.4))
sol2

$\left(-0.642870824741,\,0.385398567605\right)$
x2, y2 = sol2[0], sol2[1]
x2, y2

$\left(-0.642870824741, 0.385398567605\right)$
f1(x2,y2), f2(x2,y2)

$\left(-3.02336881641 \times 10^{-07}, -1.0387405891 \times 10^{-07}\right)$
x1 - x2, y1 - y2

$\left(-1.81852617853195 \times 10^{-7}, -9.97466205743258 \times 10^{-8}\right)$