Accolade l'une des deux racines de l'algorithme de recherche de racine pour les racines d'une fonction multivariée

Toutes mes excuses pour le titre (peut-être trompeur) et la question probablement déroutante elle-même, j'ai beaucoup de mal à formuler mon problème et surtout à le compresser en une seule phrase pour le titre. Je veux trouver les racines d'une fonction f (w, t, some_other_args) avec deux variables, w et t , en utilisant python. La vraie structure des fonctions est vraiment longue et compliquée, vous pouvez la trouver à la fin de ce post. L'important est qu'il contienne la ligne suivante:

import math as m
from scipy import optimize
import numpy as np


def func(t, w, r_1, r_2, r_3):

    k = 1.5 * m.sqrt((1.0 - w) / (1.0 - 0.25 * w))

    k23 = 2 * k / 3

    z1 = 1 / (1 + k23)
    z2 = 1 / (1 - k23)
    z3 = 3 * ((1 / 5 + r_1 - r_2 - 1 / 5 * r_1 * r_2) / (z1 - r_2 * z2)) * m.exp(t * (k - 1))
    z4 = -(z2 - r_2 * z1) / (z1 - r_2 * z2) * m.exp(2 * k * t)
    z5 = -(z1 - r_2 * z2) / (z2 - r_2 * z1)
    z6 = 3 * (1 - r_2 / 5) / (z2 - r_2 * z1)

    beta_t = r_3 / (z2 / z1 * m.exp(2 * k * t) + z5) * (z6 - 3 / (5 * z1) * m.exp(t * (k - 1)))
    alpha_t = beta_t * z5 - r_3 * z6

    beta_r = (z3 - r_1 / 5 / z2 * m.exp(-2 * t) * 3 - 3 / z2) / (z1 / z2 + z4)
    alpha_r = -z1 / z2 * beta_r - 3 / z2 - 3 / 5 * r_1 / z2 * m.exp(-2 * t)

    It_1 = 1 / 4 * w / (1 - 8 / 5 * w) * (alpha_t * z2 * m.exp(-k * t) + beta_t * z1 * m.exp(k * t) + 3 * r_3 * m.exp(-t))

    Ir_1 = (1 / 4 * w / (1 - 8 / 5 * w)) * (z1 * alpha_r + z2 * beta_r + 3 / 5 + 3 * r_1 * m.exp(-2 * t))

    T = It_1 + m.exp(-t) * r_3
    R = Ir_1 + m.exp(-2 * t) * r_1

    return [T, R]


def calc_1(t, w, T, R, r_1, r_2, r_3):
    t_begin = float(t[0])
    T_new, R_new = func(t_begin, w, r_1, r_2, r_3)
    a = abs(-1 + T_new/T)
    b = abs(-1 + R_new/R)
    return np.array([a, b])


def calc_2(x, T, R, r_1, r_2, r_3):
    t = x[0]
    w = x[1]
    T_new, R_new = func(t, w, r_1, r_2, r_3)
    a = abs(T - T_new)
    b = abs(R - R_new)
    return np.array([a, b])


def approximate_w(R):
    k = (1 - R) / (R + 2 / 3)
    w_approx = (1 - ((2 / 3 * k) ** 2)) / (1 - ((1 / 3 * k) ** 2))
    return w_approx


def approximate_t(w, T, R, r_1, r_2, r_3):
    t = optimize.root(calc_1, x0=np.array([10, 0]), args=(w, T, R, r_1, r_2, r_3))
    return t.x[0]


def solve(T, R, r_1, r_2, r_3):
    w_x = approximate_w(R)
    t_x = approximate_t(w_x, T, R, r_1, r_2, r_3)
    sol = optimize.fsolve(calc_2, x0=np.array([t_x, w_x]), args=(T, R, r_1, r_2, r_3))
    return sol


# Values for testing:
T = 0.09986490557943692
R = 0.8918728343037964
r_1 = 0
r_2 = 0
r_3 = 1

print(solve(T, R, r_1, r_2, r_3))

Cela signifie que w ne peut pas dépasser 1, car cela conduirait à calculer la racine carrée d'un nombre négatif, ce qui, bien sûr, est impossible. J'ai des algorithmes pour calculer les valeurs approximatives de w et t en utilisant d'autres valeurs dans ma fonction, mais elles sont très inexactes.

Donc, j'essaye pour calculer les racines avec scipy.optimize.fsolve (après avoir essayé littéralement tous les algorithmes de recherche de racines que j'ai pu trouver en ligne, j'ai trouvé que celui-ci était le meilleur pour ma fonction) en utilisant ces valeurs approximatives comme points de départ, qui ressemblerait à ceci:

w_approx: 0.9960090844989311
t_approx: 24.26777844720981
Values: t:24.26777844720981, w:0.9960090844989311
Values: t:24.26777844720981, w:0.9960090844989311
Values: t:24.26777844720981, w:0.9960090844989311
Values: t:24.267778808827888, w:0.9960090844989311
Values: t:24.26777844720981, w:0.996009099340623
Values: t:16.319554685876746, w:1.0096680915775516
      solution = optimize.fsolve(f, x0=np.array([t_approx, w_approx]), args=(some_other_args))
    File "C:\Users\...\venv\lib\site-packages\scipy\optimize\minpack.py", line 148, in fsolve
      res = _root_hybr(func, x0, args, jac=fprime, **options)
    File "C:\Users\...\venv\lib\site-packages\scipy\optimize\minpack.py", line 227, in _root_hybr
      ml, mu, epsfcn, factor, diag)
    File "C:\Users\...\algorithm.py", line 9, in f
      k = 1.5 * m.sqrt((1.0 - w) / (1.0 - 0.25 * w))
ValueError: math domain error

Pour la plupart des valeurs, cela fonctionne parfaitement. Si w est trop proche de 1, cependant, il arrive toujours un moment où fsolve essaie une valeur supérieure à 1 pour w , qui, dans turn, lève une ValueError (car calculer la racine d'un nombre négatif est mathématiquement impossible). Voici un exemple d'impression des valeurs que fsolve utilise, où w devrait être quelque part autour de 0.997:

solution = optimize.fsolve(f, x0=np.array([t_approx, w_approx]), args=(some_other_args))

Alors, comment puis-je dire à Optimize.fsolve que w ne peut pas être plus grand que 1? Ou quels sont les algorithmes alternatifs pour faire quelque chose comme ça (je connais brentq et ainsi de suite, mais tous nécessitent de donner un intervalle pour les deux racines, ce que je ne sais pas que vous voulez faire.)?

Code de test (Ce qui est important à noter ici: même si func est théoriquement censé calculer R et T étant donné t et w , je dois l'utiliser dans l'autre sens. C'est un peu maladroit, mais je n'y parviens tout simplement pas réécrire la fonction pour qu'elle accepte T, R pour calculer t, w - c'est un peu trop pour mon expertise mathématique médiocre;)):

k = 1.5 * m.sqrt((1.0 - w) / (1.0 - 0.25 * w))