找到指数求和的解
Finding the solution of a summation of exponentials
我正在尝试对 Python 中的这个方程进行数值求解(numpy/scipy,一切都可用)
这个公式中K是常数,f和g是两个项取决于 E 计数器(这是积分的离散表示),其中 x 是我正在寻找的变量。
例如,E 是 3 个术语:
也知道 f(E) 和 g(E)。
我从 numpy 中读到了有关使用 "fsolve" 的信息,尽管我不明白如何自动生成作为项求和的函数。我可能会手动完成,但需要一段时间的 50 个术语,我也很想学习一些新东西。
您可以使用 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))
请注意,对于您的特定函数,您还可以通过提供函数的导数来辅助算法:
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)
我正在尝试对 Python 中的这个方程进行数值求解(numpy/scipy,一切都可用)
这个公式中K是常数,f和g是两个项取决于 E 计数器(这是积分的离散表示),其中 x 是我正在寻找的变量。
例如,E 是 3 个术语:
也知道 f(E) 和 g(E)。
我从 numpy 中读到了有关使用 "fsolve" 的信息,尽管我不明白如何自动生成作为项求和的函数。我可能会手动完成,但需要一段时间的 50 个术语,我也很想学习一些新东西。
您可以使用 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))
请注意,对于您的特定函数,您还可以通过提供函数的导数来辅助算法:
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)