Find root of a transcendental equation with python

Question

I have to solve the following transcendental equation

cos(x)/x=c

for given constant c.

For example I did a short code in Mathematica, where I generated a list of random values for constant c

const = Table[RandomReal[{0, 5}], {i, 1, 10}]

(*{1.67826, 0.616656, 0.290878, 1.10592, 0.0645222, 0.333932, 3.59584, \
2.70337, 3.91535, 2.78268}*)

Than I defined the function

f[x_, i_] := Cos[x]/x - const[[i]]

and started looking for the roots:

Table[FindRoot[f[x, i] == 0, {x, 0.1}][[1, 2]], {i, 1, Length[const]}]
(*{0.517757, 0.947103, 1.21086, 0.694679, 1.47545, 1.16956, 0.26816, \
0.347764, 0.247615, 0.338922}*)

---

Now I would love to programme something similar in python (probably using numpy?) but I can't really find any good existing answer to a problem like that. Could somebody help?

Answer 1

One way that I have achieved this in the past is to use scipy.optimize.minimize to find the minima of the squared function.

from scipy.optimize import minimize
from numpy import cos

def opt_fun(x, c):
    return (cos(x)/x - c)**2

const = 1.2
res = minimize(lambda x: opt_fun(x, const), x0=0.001)

# Check if the optimization was successful
print(res.success)
# >> True

# Extract the root from the minimization result
print(res.x[0])
# >> 0.65889256782472172

This is by no means fool-proof, but it can be quick and accurate. If there are multiple roots, for instance, minimize will find the one in the 'downhill direction' from the initial point you select which is why I've chosen a small positive value above.

One other issue to keep an eye out for, which is always true with minimization problems, is numbers with dramatically different orders of magnitude. In your equation, as c gets very large, the first positive root gets very small. If you wind up trying to find roots in that circumstance, you may need to scale both x to be near to 1 in order to get accurate results (an example here).

Answer 2

Alternatively, you can use root:

import numpy as np
from scipy.optimize import root

def func_cos(x, c):
    return np.cos(x) / x - c

crange = range(1, 11)

res = [root(func_cos, 0.5, args=(ci, )).x[0] for ci in crange]

Then res looks as follows:

[0.73908513321516056,
 0.45018361129487355,
 0.31675082877122118,
 0.24267468064089021,
 0.19616428118784215,
 0.16441893826043114,
 0.14143076140757282,
 0.12403961812459068,
 0.11043425911223313,
 0.099505342687387879]

If you are interested in solving systems of equations using root you can check this answer.

Answer 3

For this type of simple, univariate functions, you can easily find all the roots within an interval of interest using a python implementation of Chebfun. I am aware of two, Chebpy and pychebfun, which are both excellent.

For example, using Chebpy, one would do the following to find the roots of cos(x)/x - 0.05 in the interval [0.5, 12]:

from chebpy import chebfun

x = chebfun('x', [0.5, 12])
c = 0.05
f = np.cos(x)/x - c

rts = f.roots()
print(rts)

[ 1.4959 4.9632 7.4711 11.6152]

Answer 4

You can do this with sympy:

>>> from sympy import cos, Symbol, nsolve
>>> x = Symbol('x')
>>> consts = [random.random() for _ in range(10)]
>>> [nsolve(cos(x)/x - c, x, 1) for c in consts]
[mpf('0.89659506789294669'),
 mpf('0.96201114853313738'),
 mpf('0.74186728791161379'),
 mpf('1.1720944924353926'),
 mpf('0.92953351945607071'),
 mpf('0.96626530553984035'),
 mpf('1.4270719610604761'),
 mpf('0.85968954499458035'),
 mpf('0.86682911058530746'),
 mpf('0.91591678333479274')]