Решите уравнение, используя численный решатель python в numpy

У меня есть уравнение:

R - ((1.0 - np.exp(-tau))/(1.0 - np.exp(-a*tau))) = 0.

Я хочу решить для tau в этом уравнении, используя численный решатель, доступный в numpy. Каков наилучший способ сделать это?

Значения для R и a в этом уравнении изменяются для разных реализаций этой формулы, но фиксируются при определенных значениях, когда они должны быть решены для tau.

Ответ 1

В обычных математических обозначениях ваше уравнение

$R = \frac{1 - e^{-\tau}}{1 - e^{-a\cdot\tau}}$

Функция SciPy fsolve выполняет поиск точки, в которой данное выражение равно нулю ( "нуль" или "корень" выражения). Вам нужно предоставить fsolve первоначальное предположение, что "рядом" с вашим желаемым решением. Хорошим способом найти такое первоначальное предположение является просто построение выражения и поиск нулевого пересечения.

#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# Define the expression whose roots we want to find

a = 0.5
R = 1.6

func = lambda tau : R - ((1.0 - np.exp(-tau))/(1.0 - np.exp(-a*tau))) 

# Plot it

tau = np.linspace(-0.5, 1.5, 201)

plt.plot(tau, func(tau))
plt.xlabel("tau")
plt.ylabel("expression value")
plt.grid()
plt.show()

# Use the numerical solver to find the roots

tau_initial_guess = 0.5
tau_solution = fsolve(func, tau_initial_guess)

print "The solution is tau = %f" % tau_solution
print "at which the value of the expression is %f" % func(tau_solution)

Ответ 2

Вы можете переписать уравнение как

Для целых a и ненулевых R вы получите решения a в сложном пространстве;
Существуют аналитические решения для a=0,1,...4 (см. здесь);

Таким образом, в общем случае у вас может быть одно, многократное или никакое решение, а некоторые или все из них могут быть сложными. Вы можете легко бросить scipy.root в этом уравнении, но ни один численный метод не позволит найти все решения.

Решить в сложном пространстве:

import numpy as np
from scipy.optimize import root

def poly(xs, R, a):
    x = complex(*xs)
    err = R * x - x + 1 - R
    return [err.real, err.imag]

root(poly, x0=[0, 0], args=(1.2, 6))