5

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?

Willem Van Onsem
  • 443,496
  • 30
  • 428
  • 555
skrat
  • 648
  • 2
  • 10
  • 27
  • Is it critical to solve this using random samples? If no, have a look at [Newton's method](https://en.wikipedia.org/wiki/Newtons_method) – Michail Mar 27 '17 at 13:21
  • @Michail, No, definitely not! I only used random generator to provide some synthetic data. In real case I will have some experimental data. – skrat Mar 27 '17 at 13:27
  • 1
    How do you want to deal with multiple roots? Is it important to find all of them be found or do you want to find the one closes to 0? – Chris Mueller Mar 27 '17 at 13:28
  • 1
    @ChrisMueller I am only interested in the first positive root. No more. – skrat Mar 27 '17 at 13:30
  • There is a Bessel function series solution. See [Kepler equation](https://en.wikipedia.org/wiki/Kepler%27s_equation) – Tyma Gaidash Jun 12 '23 at 18:50

4 Answers4

6

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.

Community
  • 1
  • 1
Cleb
  • 25,102
  • 20
  • 116
  • 151
5

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

Community
  • 1
  • 1
Chris Mueller
  • 6,490
  • 5
  • 29
  • 35
  • 2
    scipy.optimize also offers various [root finding functions](https://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding), although I doubt that, for this example, your proposed optimization-based approach performs any worse. – Stelios Mar 27 '17 at 14:13
  • @Stelios: I used such a root finding algorithm below. The value for the example with `c=1.2` is indeed almost identical. – Cleb Apr 03 '17 at 15:46
  • Finding a root of an equation by looking for minima of the function squared isn't really recommendable, since it introduces an ambiguity between roots and zero-value minima of the original function, and the accuracy of the minima is something like the square root fo the accuracy to which the roots can be determined. Directly finding the roots via methods mentioned by other answers is to be preferred. – Robert Dodier Mar 25 '21 at 18:23
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]

enter image description here

Stelios
  • 5,271
  • 1
  • 18
  • 32
  • This looks like a really elegant solution but I am having so much troubles installing the ChebPy on windows...you can't even imagine. – skrat Mar 27 '17 at 16:02
  • Do you have any experience with installing it on windows 10, python 3.6.0? – skrat Mar 27 '17 at 16:21
  • @skrat Judging from the installation instructions of cheby, windows may not be supported. Maybe try pychebfun, which I had successfully installed on a windows machine with python 3.5. – Stelios Mar 27 '17 at 17:01
2

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')]
AChampion
  • 29,683
  • 4
  • 59
  • 75
  • 1
    Do you have any idea how this handles multiple roots? If you set `c=0.1`, for instance, this equation appears to have ~10 roots. – Chris Mueller Mar 27 '17 at 13:43
  • `nsolve` only solves for one root, behind the scenes it uses `mpmath.findroot()` so you maybe able to use that lower level library directly. – AChampion Mar 27 '17 at 13:51