résolvant un système d'équations transcendantes avec python

Question

En supposant que je les quatre équations suivantes:

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

pour les variables inconnues x, y, a et b. Notez que cos(x)/x=a a plusieurs solutions. Semblable va pour la variable y. Je ne m'intéresse qu'aux valeurs x et y, qui sont les premières racines positives (si cela compte).

Vous pouvez supposer que a, b, c et d sont des constantes réelles connues, toutes positives.

En Mathematica le code pour résoudre ce ressemblerait à quelque chose comme:

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}}] 

qui à la suite retourne

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

Bien que ce fut assez facile, je ne sais rien faire comme ça en python. Donc, si quelqu'un pouvait me guider (ou simplement me montrer comment) pour résoudre cela, je l'apprécierais grandement.

Answer 1

Vous pouvez utiliser 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) 

qui imprimeront

[ 1.30301572 1.30987969 0.51530547 0.48469453] 

Vos fonctions doivent être définies d'une manière qu'ils évaluent à 0, par exemple a + b - 1 au lieu de a + b = 1.

Une vérification rapide:

print(your_funcs(sol2.x)) 

donne

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

Alors, devrait être ok la solution (s'il vous plaît noter que e-11 est essentiellement 0).

Vous pouvez également utiliser fsolve:

from scipy.optimize import fsolve 

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

qui vous donne le même résultat:

[ 1.30301572 1.30987969 0.51530547 0.48469453] 

Vous pouvez passer des arguments supplémentaires en utilisant l'argument args:

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) 

qui vous donne l'ancienne sortie:

[ 1.30301572 1.30987969 0.51530547 0.48469453] 

Si vous exécutez

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

vous recevez:

[ 1.26670224 1.27158794 0.34096159 0.65903841]