如何求解具有两个变量的指数方程?

How to solve exponential equation with two variables?

背景

我有一个这样的指数方程:

a - b * np.exp(-c/x) - y * np.exp(-delta/x).sum() * 2

其中 a、b 和 c 是常量,delta 是一维数组,可从 Google Drive 获得。目标是解决 xy.

尝试

为了求解有边界的方程,我想出了 optimize.minimize 来得到 xy,其中残差的平方最小。

完整示例如下:

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

a = 35167.7
b = 11919.5
c = 1.68
delta = np.load('delta.npy')

def residual_sum(x0, a, b, c, delta):
    x = x0[0]
    y = x0[1]

    residual = a - b * np.exp(-c/x) - y * np.exp(-delta/x).sum() * 2

    return residual**2

bnds = ((1, 24), (30, 1500))
x_init = 3
y_init = 300
x0 = [x_init, y_init]

solution = opt.minimize(residual_sum, x0, bounds=bnds, args=(a, b, c, delta))
print(solution['x'][0], solution['x'][1])

但是,当我更改初始值时,结果却大不相同:

solution = opt.minimize(residual_sum, [3, 300], bounds=bnds, args=(a, b, c, delta))
print(solution['x'][0], solution['x'][1])

solution = opt.minimize(residual_sum, [5, 500], bounds=bnds, args=(a, b, c, delta))
print(solution['x'][0], solution['x'][1])

Output:
10.110838104427442 87.90413009609203
24.0 80.08308172149127

问题

所以,我用手动输入 x 和 y 检查方程:

x = np.linspace(0, 24, 200)
y = np.linspace(0, 1500, 500)

res = (a - b * np.exp(-c/x))[:, None] - np.exp(-delta[:,np.newaxis]/x).sum(axis=0)[:, None] * y * 2

X, Y = np.meshgrid(x, y)

fig, ax = plt.subplots()

m = ax.pcolormesh(X, Y, (res**2).T, cmap='viridis', norm=matplotlib.colors.LogNorm())

bnds = ((1, 24), (30, 1500))
x_init = 3
y_init = 300
x0 = [x_init, y_init]

solution = opt.minimize(residual_sum, x0, bounds=bnds, args=(a, b, c, delta))
ax.scatter(solution['x'][0], solution['x'][1], marker='*', facecolors='none', edgecolors='r')

plt.colorbar(m)

图中有明显的低值线。知道如何获得“正确”的解决方案吗?

这是具有较低 vmax 的相同数据:

m = ax.pcolormesh(X, Y, (res**2).T, cmap='viridis', vmax=2e3)

您的等式如下所示:

a - b * np.exp(-c/x) - y * np.exp(-delta/x).sum() * 2 = 0

也就是说:

y = (a - b * np.exp(-c/x)) / np.exp(-delta/x).sum() / 2

您可以为 x 插入任何值并获得相应的 y:

>>> import numpy as np; import matplotlib.pyplot as plt
>>> a,b,c = 35167.7, 11919.5, 1.68
>>> delta  = np.load('delta.npy')
>>> def the_y(x): return (a - b * np.exp(-c/x)) / np.exp(-delta/x).sum() / 2
... 
>>> import matplotlib.pyplot as plt
>>> X = np.linspace(1, 24, 1000)
>>> plt.plot(X, [the_y(x) for x in X]); plt.show()

情节如下:

确实,对于 X = np.linspace(1, 24, 1000) 中的每个 x,我得到了相应的 y 值。您可以生成无数个 x 并得到无数个 y 作为响应,所以我想说有无限多的解决方案。