CODEWITHSUNDEEP
pythonnumpywolfram-mathematicaroot
skrat
skrat

Reputation: 748

Solution to transcendental equation from Mathematica and Python do not match

I am solving a system of transcendental equations:

  1. cos(x) / x = 0.48283 + a*3.46891
  2. cos(y) / y = 0.47814 + b*28.6418
  3. a + b = 1
  4. 1.02 * sinc(x) = 1.03 * sinc(y)

And it just so happens I tried to solve the above system in two separated programming languages (Mathematica and Python)

Mathematica

Running the code

FindRoot[{Cos[x]/x == 0.482828 + a*3.46891, 
  Cos[y]/y == 0.47814 + b*28.6418, a + b == 1, 
  1.02*Sinc[x] == 1.03*Sinc[y]}, {{x, .2}, {y, .2}, {a, 0.3}, {b, 
   0.3}}, PrecisionGoal -> 6]

returns

{x -> 0.261727, y -> 0.355888, a -> 0.924737, b -> 0.0752628}

Python

Running the code:

import numpy as np
from scipy.optimize import root

def fuuu(X, en,dv,tri,sti):
        x, y, a, b = X

        F = [np.cos(x) / x - en-a*dv,
             np.cos(y) / y - tri-b*sti,
             a + b - 1,
             1.02 * np.sinc(x) - 1.03 * np.sinc(y)]

        return F
    root(fuuu, [0.2, 0.2, 0.3, 0.3], args=(0.482828,3.46891,0.47814,28.6418)).x

returns

array([ 0.26843418,  0.27872813,  0.89626625,  0.10373375])

Comparison

Let's say that the 'x' value is the same. Let's just ignore the small difference. But the y values differ by miles! The physical meaning completely changes. For some reason I believe the values from Mathematica more than I believe values from Python.

Questions:

  1. Why do the calculations differ?
  2. Which one is now correct? What do I have to change in python (assuming python is the problematic one)?

Upvotes: 1

Views: 1018

Answers (1)

kennytm
kennytm

Reputation: 523274

The calculation differ because of the sinc function.

(* Mathematica *)
In[1] := Sinc[0.26843418]
Out[1] = 0.988034

# Python
>>> np.sinc(0.26843418)
0.88561519683835599
>>> np.sin(0.26843418) / 0.26843418
0.98803370932709034

Huh? Well let's RTFM

numpy.sinc(x)

Return the sinc function.

The sinc function is sin(πx)/(πx).

Oops. NumPy's sinc has a different definition than Mathematica's Sinc.

  • Mathematica's Sinc uses the unnormalized definition sin(x)/x. This definition is usually used in mathematics and physics.
  • NumPy's sinc uses the normalized version sin(πx)/(πx). This definition is usually used in digital signal processing and information theory. It is called normalized because
    ∫-∞∞ sin(πx)/(πx) dx = 1.

Therefore, if you want NumPy to produce the same result as Mathematica, you need to divide x and y by np.pi.

def fuuu(X, en,dv,tri,sti):
    x, y, a, b = X
    F = [np.cos(x) / x - en-a*dv,
         np.cos(y) / y - tri-b*sti,
         a + b - 1,
         1.02 * np.sinc(x/np.pi) - 1.03 * np.sinc(y/np.pi)]    # <---
    return F

>>> root(fuuu, [0.2, 0.2, 0.3, 0.3], args=(0.482828,3.46891,0.47814,28.6418)).x
array([ 0.26172691,  0.3558877 ,  0.92473722,  0.07526278])

Upvotes: 7

Related Questions