找到指数总和的解 [英] Finding the solution of a summation of exponentials
问题描述
我正在尝试用 Python 对这个方程进行数值求解(numpy/scipy,一切都可用)
I'm trying to solve numerically this equation in Python (numpy/scipy, everything is available)
在这个公式中,K 是一个常数,f 和 g 是两个依赖于 E 计数器的术语(这是一个离散表示一个积分),其中 x 是我正在寻找的变量.
In this formula K is a constant, f and g are two terms those depends on the E counter (this is a discrete representation of an integral) where x is the variable I'm looking for.
例如,E 是 3 个术语:
As an example, with E being 3 terms that'd be:
还有 f(E) 和 g(E) 是已知的.
also f(E) and g(E) are known.
我从 numpy 中了解到使用fsolve",但我不明白如何自动生成一个函数,该函数是术语的总和.我可以手动完成,但 50 个术语需要一段时间,而且我很想学习新的东西.
I read about using "fsolve" from numpy, although I can't understand how to automatically generate a function that's a summation of terms. I may do it manually but being 50 terms that'd take a while, also I'd love to learn something new.
推荐答案
您可以使用 scipy.optimize.fsolve
其中函数使用 numpy.sum
:
You can use scipy.optimize.fsolve
where the function is constructed using 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))
请注意,对于您的特定函数,您还可以通过提供函数的导数来辅助算法:
Note that for your specific function you can also assist the algorithm by providing the derivative of the function:
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 个输出(对于每一行再次产生一个).因此输入和输出是相互独立的,对应的雅可比是一个对角矩阵.下面是一些示例代码:
Finding multiple roots
You can vectorize the procedure in a sense that the function accepts N inputs (for example for each of the rows) and produces N outputs (again one for each row). Hence the inputs and outputs are independent of each other and the corresponding jacobian is a diagonal matrix. Here's some sample code:
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)
这篇关于找到指数总和的解的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!