Нахождение решения суммирования экспонент - PullRequest
Новая
       35

Нахождение решения суммирования экспонент

0 голосов
/ 18 января 2019

Я пытаюсь численно решить это уравнение в Python (numpy / scipy, все доступно)

image

В этой формуле K - это константа, f и g - два члена, которые зависят от счетчика E (это дискретное представление интеграла ) где x - переменная, которую я ищу.

Например, с E 3 терминами, которые будут:

image

также f (E) и g (E) .

Я читал об использовании "fsolve" из numpy, хотя не могу понять, как автоматически генерировать функцию, которая является суммированием терминов. Я могу сделать это вручную, но из 50 терминов, которые потребуют времени, я также хотел бы узнать что-то новое.

1 Ответ

0 голосов
/ 18 января 2019

Вы можете использовать scipy.optimize.fsolve, если функция построена с использованием numpy.sum: 

import numpy as np
import scipy.optimize

np.random.seed(123)

f = np.random.uniform(size=50)
g = np.random.uniform(size=f.size)
K = np.sum(f * np.exp(-g*np.pi))

def func(x, f, g, K):
    return np.sum(f * np.exp(-g*x), axis=0) - K

# The argument to the function is an array itself,
# so we need to introduce extra dimensions for f, g.
res = scipy.optimize.fsolve(func, x0=1, args=(f[:, None], g[:, None], K))

Обратите внимание, что для вашей конкретной функции вы также можете помочь алгоритму, предоставив производную функции: 

def derivative(x, f, g, K):
    return np.sum(-g*f * np.exp(-g*x), axis=0)

res = scipy.optimize.fsolve(func, fprime=derivative,
                            x0=1, args=(f[:, None], g[:, None], K))

Нахождение нескольких корней

Вы можете векторизовать процедуру в том смысле, что функция принимает N входных данных (например, для каждой строки) и производит N выходных данных (опять же по одному для каждой строки). Следовательно, входы и выходы не зависят друг от друга, и соответствующий якобиан является диагональной матрицей. Вот пример кода: 

import numpy as np
import scipy.optimize

np.random.seed(123)

image = np.random.uniform(size=(4000, 3000, 2))
f, g = image[:, :, 0], image[:, :, 1]
x_original = np.random.uniform(size=image.shape[0])  # Compute for each of the rows.
K = np.sum(f * np.exp(-g * x_original[:, None]), axis=1)

def func(x, f, g, K):
    return np.sum(f * np.exp(-g*x[:, None]), axis=1) - K

def derivative(x, f, g, K):
    return np.diag(np.sum(-g*f * np.exp(-g*x[:, None]), axis=1))

res = scipy.optimize.fsolve(func, fprime=derivative,
                            x0=0.5*np.ones(x_original.shape), args=(f, g, K))
assert np.allclose(res, x_original)
...