Solving a system of transcendental equations with python

Question

Assuming I have the following four equations:

1. cos(x)/x = a
2. cos(y)/y = b
3. a + b = 1
4. c sinc(x) = d sinc(y)

for unknown variables x, y, a and b. Note that cos(x)/x=a has multiple solutions. Similar goes for variable y. I am only interested in x and y values, which are first positive roots (if that matters).

You can safely assume a, b, c and d are known real constants, all positive.

In Mathematica the code to solve this would look something like:

FindRoot[{Cos[x]/x == 0.2 a + 0.1, 
          Cos[y]/y == 0.2 b + 0.1, 
          a + b == 1.0, 
           1.03*Sinc[x] == Sinc[y]*1.02}, 
          {{x, .1}, {y, .1}, {a, .3}, {b, .1}}]

which as a result returns

{x -> 1.31636, y -> 1.29664, a -> 0.456034, b -> 0.543966}

While this was quite easy, I have no idea how to do anything like that in python. So if somebody could kinda guide me (or simply show me how) to solve this, I would highly appreciate it.

Answer 1

You can use root:

import numpy as np
from scipy.optimize import root

def your_funcs(X):

    x, y, a, b = X

    f = [np.cos(x) / x - 0.2 * a - 0.1,
         np.cos(y) / y - 0.2 * b - 0.1,
         a + b - 1,
         1.03 * np.sinc(x) - 1.02 * np.sinc(y)]

    return f

sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1])
print(sol2.x)

that will print

[ 1.30301572  1.30987969  0.51530547  0.48469453]

Your functions have to be defined in a way that they evaluate to 0, e.g. a + b - 1 instead of a + b = 1.

A quick check:

print(your_funcs(sol2.x))

gives

[-1.9356960478944529e-11, 1.8931356482454476e-11, 0.0, -4.1039033282785908e-11]

So, the solution should be ok (please note that e-11 is basically 0).

Alternatively, you can also use fsolve:

from scipy.optimize import fsolve

sol3 = fsolve(your_funcs, [0.1, 0.1, 0.3, 0.1])

which gives you the same result:

[ 1.30301572  1.30987969  0.51530547  0.48469453]

You can pass additional arguments using the args argument:

def your_funcs(X, fac_a, fac_b):

    x, y, a, b = X

    f = [np.cos(x) / x - fac_a * a - 0.1,
         np.cos(y) / y - fac_b * b - 0.1,
         a + b - 1,
         1.03 * np.sinc(x) - 1.02 * np.sinc(y)]

    return f

sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1], args=(0.2, 0.2))
print(sol2.x)

which gives you the "old" output:

[ 1.30301572  1.30987969  0.51530547  0.48469453]

If you run

sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1], args=(0.4, 0.2))
print(sol2.x)

then you receive:

[ 1.26670224  1.27158794  0.34096159  0.65903841]