1
votes

I can be considered pretty much new to python and coding in general so forgive me for my ignorance.

I'm trying to solve a system of trigonometric functions in python, and I'm doing so using the solve command from sympy. However, this method returns only a finite number of solutions, two in this particular case. I've read through the documentation and it seems that to get an expression for all the solutions solveset is to be used instead. However, I do not want all the solutions to be displayed, but rather only a finite amount which is contained within a certain range.

Here's the example:

from sympy import *


x, y = symbols('x, y')

eq1 = Eq(y - sin(x), 0)
eq2 = Eq(y - cos(x), 0)

sol = solve([eq1, eq2], [x, y])
print(sol)

which only returns the first two solutions in the positive x range.

How could I do to, for example, display all the solutions within the x range [-2pi, 2pi]?

I'd want them in explicit form rather than written in term of some multiplier since I then need to convert them into numerical form.

Thank you in advance.

3
are you looking for the zero-crossings of sin and cos?curlycharcoal
No, I'm looking for the intersections between sin and cos.Edoardo Serra

3 Answers

2
votes

By using solveset you can restrict the solutions with domain argument. For evaluating numerical results use .evalf() or another similar method.

from sympy import Interval, symbols, solveset, sin, cos, pi


x = symbols('x')
sol = solveset(cos(x) - sin(x), x, domain=Interval(-2*pi, 2*pi))
print(sol)
print(sol.evalf())

Output

FiniteSet(-7*pi/4, -3*pi/4, pi/4, 5*pi/4)
FiniteSet(-5.49778714378214, -2.35619449019234, 0.785398163397448, 3.92699081698724)

I hope this helps!

1
votes

SymPy can really take you down rabbit holes. I agree with kampmani's solution, only if you can easily solve for y on your own. However, in more general cases and in higher dimensions, his solution does not hold.

For example, the following will be slightly more tricky:

eq1 = Eq(z - x*y, 0)
eq2 = Eq(z - cos(x) - sin(y), 0)
eq3 = Eq(z + x*y, 0)

So here I am; killing a fly with a bazooka. The problem is that one was able to simplify the set of equations into a single equation with a single variable. But what if you can't do that (for example, if it was a larger system)?

In this case, one needs to use nonlinsolve to solve the system of equations. But this does not provide numeric solutions directly and does not have a domain argument.

So the following code unpacks the solutions. It goes through each tuple in the set of solutions and finds the numeric solutions for each component in the tuple. Then in order to get the full list, you need a Cartesian Product of each of those components. Repeat this for each tuple in the set of solutions.

The following should work for almost any system of equations in any dimension greater than 1. It produces numeric solutions in the cube whose boundaries are the domains variable.

from sympy import *
import itertools  # used for cartesian product

x, y, z = symbols('x y z', real=True)
domains = [Interval(-10, 10), Interval(-10, 10), Interval(-10, 10)]  # The domain for each variable

eq1 = z - x*y
eq2 = z - cos(x) - sin(y)
eq3 = z + x*y


solutions = nonlinsolve([eq1, eq2, eq3], [x, y, z])  # the recommended function for this situation
print("---------Solution set----------")
print(solutions)  # make sure the solution set is reasonable. If not, assertion error will occur

_n = Symbol("n", integer=True)  # the solution set often seems to contain these symbols
numeric_solutions = []
assert isinstance(solutions, Set)  # everything that I had tried resulted in a FiniteSet output

for solution in solutions.args:  # loop through the different kinds of solutions
    assert isinstance(solution, Tuple)  # each solution should be a Tuple if in 2D or higher

    list_of_numeric_values = []  # the list of lists of a single numerical value
    for i, element in enumerate(solution):

        if isinstance(element, Set):
            numeric_values = list(element.intersect(domains[i]))
        else:  # assume it is an Expr
            assert isinstance(element, Expr)
            if _n.name in [s.name for s in element.free_symbols]:  # if n is in the expression
                # change our own _n to the solutions _n since they have different hidden
                # properties and they cannot be substituted without having the same _n
                _n = [s for s in element.free_symbols if s.name == _n.name][0]
                numeric_values = [element.subs(_n, n)
                                  for n in range(-10, 10)  # just choose a bunch of sample values
                                  if element.subs(_n, n) in domains[i]]
            elif len(element.free_symbols) == 0:  # we just have a single, numeric number
                numeric_values = [element] if element in domains[i] else []
            else:  # otherwise we just have an Expr that depends on x or y
                # we assume this numerical value is in the domain
                numeric_values = [element]
        # note that we may have duplicates, so we remove them with `set()`
        list_of_numeric_values.append(set(numeric_values))

    # find the resulting cartesian product of all our numeric_values
    numeric_solutions += itertools.product(*list_of_numeric_values)

# remove duplicates again to be safe with `set()` but then retain ordering with `list()`
numeric_solutions = list(set(numeric_solutions))
print("--------`Expr` values----------")
for i in numeric_solutions:
    print(list(i))  # turn it into a `list` since the output below is also a `list`.

print("--------`float` values---------")
for i in numeric_solutions:
    print([N(j) for j in i])  # could have been converted into a `tuple` instead

In particular, it produces the following output for the new problem:

---------Solution set----------
FiniteSet((0, ImageSet(Lambda(_n, 2*_n*pi + 3*pi/2), Integers), 0))
--------`Expr` values----------
[0, -5*pi/2, 0]
[0, -pi/2, 0]
[0, 3*pi/2, 0]
--------`float` values---------
[0, -7.85398163397448, 0]
[0, -1.57079632679490, 0]
[0, 4.71238898038469, 0]

It was a lot of effort and probably not worth it but oh well.

0
votes

Thanks to the brilliant suggestion from @kampmani it is possible to achieve the desired result.

For start, the FiniteSet elements are not indexed and cannot be used, so the FiniteSet has to be converted into a list:

solx_array = []
#
#
#
solx = solveset(cos(x) - sin(x), x, domain=Interval(-2*pi, 2*pi))
solx_array=list(solx)

The next step is to find the y coordinate of the intersection point given its x coordinate. The final code should look somewhat similar to this:

from sympy import Interval, symbols, solveset, sin, cos, pi

sol_array = []
x = symbols('x')

solx = solveset(cos(x) - sin(x), x, domain=Interval(-2*pi, 2*pi))
solx_array=list(solx)

for i in range(len(solx_array)):
    soly = cos(solx_array[i])
    sol_array.append(str(solx_array[i] + soly))

for i in range(len(sol_array)):
    print(sol_array[i])

Still don't know how to convert the results into numerical form though, any idea is very appreciated.