3

I am trying to solve a system of 4 exponential equations with two variables. However If I use fsolve python will only allow me two use as many equations as I have variables. But as I have infinitely many pairs of solutions (if only two equations are used) and I need to find the pair of variables that fits not only two but all four equations, fsolve does not seem to work.

So right know my code look something like this:

from scipy.optimize import fsolve
c_1 = w+(y/(x*R))
c_2 = offset - (c_1/(x*R)) #offset(1.05), w(10) and R(0.35) are constants 
def equations(p):
    x,y = p
    eq1 = (c_1*sp.exp(-y*R*1.017))/(y*R)+c_2-(x*1.017)/(y*R)-(5.1138*2*np.pi)
    eq2 = (c_1*sp.exp(-y*R*2.35))/(y*R)+c_2-(x*2.35)/(y*R)-(2.02457*4*np.pi)
    eq3 = (c_1*sp.exp(-y*R*2.683))/(y*R)+c_2-(x*2.683)/(y*R)-(6.0842178*4*np.pi)
    return (eq1,eq2,eq3)

x, y =  fsolve(equations, (1,1))

print (equations((x, y)))

However this will give me wildly different result depending on the initial guesses I give. I now want to so edit this code so that I can put additional equations to guarantee that I get the correct pair of x,y as a solution. Simply adding a third equation here in the into the function will return a TypeError of course so I was wondering how this can be done.

  • 1
    Show your real equations, not placeholders, and show example parameters so that this is reproducible. – Reinderien Mar 02 '23 at 22:58
  • Hi @Reinderien I changed the question to the actual code that I am using. Please note that these three equations are not fixed. Essentially I have about 100 equations (all of the same form as the above), but I think solving for 3-4 equations should be enough so that x and y only have one solution and not infinite. –  Mar 03 '23 at 14:40

1 Answers1

2

To add an arbitrary number of equations, I don't think you can use fsolve at all. Just run a least-squares minimisation, and make sure that you vectorise properly instead of rote repetition. The following does run and generate a result, though I have not assessed its accuracy.

import numpy as np
from scipy.optimize import minimize

offset = 1.05
w = 10
R = 0.35
a = np.array((-1.017, -2.350, -2.683))
b = np.array((
    -5.1138000*2,
    -2.0245700*4,
    -6.0842178*4,
))*np.pi


def equations(x: float, y: float) -> np.ndarray:
    c_1 = w + y/x/R
    c_2 = offset - c_1/x/R

    return (c_1*np.exp(a*y*R) + a*x)/y/R + c_2 + b


def error(xy: np.ndarray) -> np.ndarray:
    eqs = equations(*xy)
    return eqs.dot(eqs)


result = minimize(
    fun=error,
    x0=(1, 1),
    method='nelder-mead',
)
print(result.message)
print('x, y =', result.x)
print('equations:', equations(*result.x))

Optimization terminated successfully.
x, y = [0.19082078 0.13493941]
equations: [ 27.4082168   13.91087443 -42.00388728]
Reinderien
  • 11,755
  • 5
  • 49
  • 77
  • Ok Thank you very much! The equations are set up so that they should all equal 0 and not [ 27.4082168 13.91087443 -42.00388728] for the resulting x and y. How can I change it so that I get the x and y for equations: [ 0,0,0]. Or is [ 27.4082168 13.91087443 -42.00388728] the closest you can get to zero for all 3? Thank you in advance ! –  Mar 05 '23 at 16:50
  • I believe that this is the closest we can get, unless you have evidence otherwise. Do you have manually-calculated roots that do better? – Reinderien Mar 05 '23 at 16:54
  • 1
    II unfortunately do not, so I will just believe you that that is best we can get so thank you very much for your time and help! –  Mar 05 '23 at 17:15