Скобка один из двух корней в алгоритме поиска корней для корней многомерной функции

Question

Извиняюсь за (возможно, вводящий в заблуждение) заголовок и, возможно, сам вопрос, который вводит в заблуждение, я много борюсь с формулировкой своей проблемы и особенно сжимаю ее в одно предложение для заголовка. Я хочу найти корни функции f(w, t, some_other_args) с двумя переменными, w и t, используя python. Реальная структура функций действительно длинная и сложная, вы можете найти ее в конце этого поста. Важно, что он содержит следующую строку:

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

Это означает, что w не может превышать 1, потому что это привело бы к вычислению квадратного корня отрицательного числа, что, конечно, невозможно. У меня есть алгоритмы для вычисления приблизительных значений w и t с использованием других значений в моей функции, но они очень неточны.

Итак, я пытаюсь вычислить корни с помощью scipy.optimize.fsolve (после попытки буквально каждого алгоритма поиска корней, который я смог найти в Интернете, я нашел этот алгоритм наилучшим для моей функции), используя эти приблизительные значения в качестве отправных точек, что выглядеть так:

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

Для большинства значений это работает отлично. Однако, если w слишком близко к 1, всегда наступает момент, когда fsolve пытается получить значение, большее 1, для w, что, в свою очередь, повышает ValueError (поскольку вычисляется корень отрицательного значения). число математически невозможно). Это пример распечатки значений, которые использует fsolve, где w должно быть где-то около 0,997:

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

Итак, как я могу сказать optimize.fsolve, что w не может быть больше 1? Или каковы альтернативные алгоритмы для выполнения чего-то подобного (я знаю о brentq и т. Д., Но все они требуют дать интервал для обоих корней, что я не хочу делать)?

---

Код для тестирования (Здесь важно отметить: хотя теоретически func должен вычислять R и T с учетом t и w, я должен использовать его наоборот. немного неуклюже, но мне просто не удается переписать функцию, чтобы она принимала T, R для вычисления t, w - это слишком много для моей посредственной математической экспертизы;)):

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

Answer 1

Вы должны попытаться явно определить свою функцию, прежде чем оптимизировать ее, чтобы вам было проще проверять домен.

По сути, у вас есть функция T и R. это сработало для меня:

def func_to_solve(TR_vector, r_1, r_2, r_3):
    T, R = TR_vector   # what you are trying to find
    w_x = approximate_w(R)
    t_x = approximate_t(w_x, T, R, r_1, r_2, r_3)
    return (calc_2([t_x, w_x], T, R, r_1, r_2, r_3))

def solve(TR, r_1, r_2, r_3):
    sol = optimize.fsolve(func_to_solve, x0=TR, args=(r_1, r_2, r_3))
    return sol

Также замените m.exp на np.exp

Answer 2

Вы сообщили об ошибке, так как fsolve не может справиться с неявными ограничениями в преобразовании w в k. Эту проблему можно решить радикально, инвертировав эту зависимость, сделав func зависимым от t и k.

def w2k(w): return 3 * m.sqrt((1.0 - w) / (4.0 - w))
    #k = 1.5 * m.sqrt((1.0 - w) / (1.0 - 0.25 * w))
    # (k/3)**2 * (4-w)= 1-w 
def k2w(k): return 4 - 3/(1-(k/3)**2)

def func(t, k, r_1, r_2, r_3):
    w = k2w(k)
    print "t=%20.15f, k=%20.15f, w=%20.15f"%(t,k,w)
    ...

Затем удалите абсолютные значения из значений функций в calc1 и calc2. Это только делает ваши решения недифференцируемыми точками, что плохо для любого алгоритма поиска корня. Изменения знака и гладкие корни хороши для методов, подобных Ньютону.

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

Нет особого смысла находить значение для t, решая уравнение, сохраняя w соотв. k исправлено, это только удваивает вычислительные усилия.

def approximate_k(R):
    k = (1 - R) / (R + 2 / 3)
    return k

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

t,k = solve(T, R, r_1, r_2, r_3)
print "t=%20.15f, k=%20.15f, w=%20.15f"%(t, k, k2w(k))

С этими модификациями решение

t=  14.860121342410327, k=   0.026653140486605, w=   0.999763184675043

найдено в 15 оценках функций.

Answer 3

Как насчет логистика аргумент, который вы хотите ограничить? Я имею в виду, внутри f вы могли бы сделать

import numpy as np

def f(free_w, ...):
    w = 1/(1 + np.exp(-free_w)) # w will always lie between 0 and 1
    ...
    return zeros

И тогда вам просто нужно применить то же логистическое преобразование к значению решения free_w, чтобы получить w*. См

solution = optimize.fsolve(f, x0=np.array([t_approx, w_approx]), args=(some_other_args))
free_w   = solution[0]
w        = 1/(1 + np.exp(-free_w))