如何求解具有两个变量的指数方程?
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 获得。目标是解决 x
和 y
.
尝试
为了求解有边界的方程,我想出了 optimize.minimize
来得到 x
和 y
,其中残差的平方最小。
完整示例如下:
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
作为响应,所以我想说有无限多的解决方案。
背景
我有一个这样的指数方程:
a - b * np.exp(-c/x) - y * np.exp(-delta/x).sum() * 2
其中 a、b 和 c 是常量,delta 是一维数组,可从 Google Drive 获得。目标是解决 x
和 y
.
尝试
为了求解有边界的方程,我想出了 optimize.minimize
来得到 x
和 y
,其中残差的平方最小。
完整示例如下:
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
作为响应,所以我想说有无限多的解决方案。