pythonfsolve

I need to solve an equation numerically, but fsolve gives me a seemingly incorrect answer

I need to solve a single variable in an equation numerically. I tried using fsolve on two different functions that are, according to my understanding, equivalent. Call these functions func1 and func2. If I specify the variable I am solving for, both functions return the same value (the residual of the equation). However when I don't specify the variable and use fsolve to find it, I get different answers depending on whether I use func1 or func2. What am I doing wrong?

data for my question

dHi=array([-125790,49080,4.2])     #  [n butane :  1,3 butadiene  :   H2]
dGi=array([-16570,124520,17.6])
V=array([-1,1,2])
No=array([1,0,0])

dH=sum(V*dHi)
dG=sum(V*dGi)

now function 1

def func1(e):
    R=8.314
    T1=298
    T2=925
    Nt=1+2*e
    Ni=array([1-e,e,2*e])
    
    lnk1=(-dG/(R*T1))
    lnk2=-(dH/R)*(1/T2 - 1/T1)+lnk1
    k2=exp(lnk2)
    
    A1=prod((Ni/Nt)**V)-k2
    
    
    return A1

for function 2 I wrote a separate function that does not require me to specify Ni, but calculates it as a function of e.

def N(e):
    return No+e*V

def func2(e):
    R=8.314
    T1=298
    T2=925
    Nt=1+2*e
     
    lnk1=(-dG/(R*T1))
    lnk2=-(dH/R)*(1/T2 - 1/T1)+lnk1
    k2=exp(lnk2)
    
    A1=prod(((N(e))/Nt)**V)-k2
    
    
    return A1

to prove N(e) and Ni is equivalent

e=0.1
Ni=array([1-e,e,2*e])
print(Ni,N(e))

I get

[0.9 0.1 0.2] [0.9 0.1 0.2]

Now to compare func1 and func2

print(fsolve(func1,0.03), fsolve(func2,0.03))

[0.10045184] [0.03108138]

If I check the second answer with both functions..

print(func1(0.03108138),func2(0.03108138))

1.2794325793047054e-11 1.2794325793047054e-11

So e = 0.03108138, and both functions can confirm this, but why does fsolve give the wrong answer for func1 ?

Solution

  • The function you pass to scipy.optimize.fsolve is supposed to accept a 1-dimensional array, and return a 1-dimensional array of the same length. (This doesn't mean it should broadcast - the function is supposed to represent a system of N nonlinear equations in N variables for some N, so the input represents N input variables and the output represents the N output values.)

    Whenever scipy.optimize.fsolve calls your function, it will pass the function an array. Even if you give an initial estimated root as a scalar, scipy.optimize.fsolve will pass your function 1-element arrays.

    Neither of your functions are designed to handle array input. By sheer luck, your func2 happens to do the right thing when passed a single-element array, but func1 builds a Ni array of the wrong shape, then broadcasts the (Ni/Nt)**V computation in a way you don't want, and ends up computing the wrong value.

    Write your functions to take and return 1-dimensional 1-element arrays instead of scalars, and you will get the right output.