solve_tour.ipynb

Kernel: SageMath 8.1

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)

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

eq.lhs()

eq.rhs()

solve(eq, x)

The solutions are returned as a list:

sol = solve(eq, x)
sol

sol[1]

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

A numerical approximation, n being a shortcut for numerical_approx:

n(sol[1].rhs())

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

A new equation involving the above solution:

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

Asking Sage whether this equation always holds:

bool(_)

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)

sol = solve((eq1, eq2), x, y)
sol

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]

sol[0][1]

sol[0][1].rhs()

n(sol[0][1].rhs())

Approximate solution: find_root

Let us consider the following transcendental equation:

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

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

solve(eq, x)

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

find_root(eq, 0, 10)

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)

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

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)

f1(x,y) = eq1.lhs()
f1

f2(x,y) = eq2.lhs()
f2

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

sol.tolist()

x1 = RR(sol.tolist()[0][0])
y1 = RR(sol.tolist()[1][0])
x1, y1

f1(x1,y1)

f2(x1,y1)

sol2 = minimize(eq1.lhs()^2 + eq2.lhs()^2, (-0.6,0.4))
sol2

x2, y2 = sol2[0], sol2[1]
x2, y2

f1(x2,y2), f2(x2,y2)

x1 - x2, y1 - y2