How to solve exponential equation with two variables?

Question

Background

I have one exponential equation like this:

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

where a, b, and c are constants, the delta is a 1D array that is available from Google Drive. The goal is to solve the x and y.

Attempt

To solve the equation with bounds, I come up with optimize.minimize to get the x and y where the square of residue is minimum.

Here's the full example:

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])

However, when I change the init values, the results are quite different:

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

Problem

So, I check the equation with manual inputs of x and 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)

There's an obvious low-value line in the figure. Any idea how to get the "correct" solution?

Here's the same data with lower vmax:

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

Answer 1

Your equation looks like this:

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

Which is to say:

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

You can plug any value for x and get the corresponding 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()

Here's the plot:

Indeed, for each x in X = np.linspace(1, 24, 1000) I got the corresponding y value. You could generate a bazillion of xs and get a bazillion of ys in response, so I'd say there are infinitely many solutions.