%display latex
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 (recall that n
is a shortcut alias 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(_)
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())
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)
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