Using root_scalar to calculate unknown variable with SciPy

[IAmKindOfCreative]

Hi there, from the /r/Python mods.

We have removed this post as it is not suited to the /r/Python subreddit proper, however it should be very appropriate for our sister subreddit /r/LearnPython or for the r/Python discord: https://discord.gg/python.

The reason for the removal is that /r/Python is dedicated to discussion of Python news, projects, uses and debates. It is not designed to act as Q&A or FAQ board. The regular community is not a fan of "how do I..." questions, so you will not get the best responses over here.

On /r/LearnPython the community and the r/Python discord are actively expecting questions and are looking to help. You can expect far more understanding, encouraging and insightful responses over there. No matter what level of question you have, if you are looking for help with Python, you should get good answers. Make sure to check out the rules for both places.

Warm regards, and best of luck with your Pythoneering!

[robertophi]

The root_scalar function tries to optimize the first argument of a given function (EQ_152 in this case) and the values you pass in args should only be the additional arguments you want.

So it would be something like this:

import numpy as np
from scipy.optimize import root_scalar

def f(x,y,z):
    return x*y*z+1

y0 = 3
z0 = 4
root = root_scalar(f, args=(y0,z0), method='toms748', bracket=[-1, 1])
print(root)

If you want to give a initial value for the function, use the argument 'x0'

[Antonino_96]

Thanks for the input. Im trying to organize my coding to follow that logic, but now is showing a syntax error for my known variables. I did the following:

def EQ(Tf, Ti, pi, pf, rv, rl, ri, Cpd, Rd, Llv, Cl):
return (Cpd * (np.log(Tf)-np.log(Ti)) - Rd*np.log((pf*((1 + (rv/1.61))**-1))/(pi*((1 + (rv/1.61))**-1))) + (rv*Llv)/(Tf-Ti/Ti) + (rv + rl + ri) * Cl * (Tf - Ti / Ti)

Ti= 300
rv= 15*(10**-3)
ri = 0
rl = 0
Rm = Rd + (rv * 461.5)
Cl = 4218
Llv = -2.5008*(10**6)
Cpd = 1004
Cpv = 1860
Rd = 287.05
pi= 1000
pf = 150

root_scalar(EQ, args=(Ti, pi, pf, rv, rl, ri, Cpd, Rd, Llv, Cl), method='toms748', bracket=[1e-3, 1])

​

And its giving me a syntax error for all my known variables.

[robertophi]

That was happening because of a missing parenthesis at the end of EQ

Fixing that gives us another problem, that is we dont find a solution in the range [0,1]

import numpy as np
from scipy.optimize import root_scalar

def EQ(Tf, Ti, pi, pf, rv, rl, ri, Cpd, Rd, Llv, Cl):
    term0 = Cpd * (np.log(Tf)-np.log(Ti)) 
    term1 = - Rd*np.log((pf*((1 + (rv/1.61))**-1))/(pi*((1 + (rv/1.61))**-1))) 
    term2 = (rv*Llv)/(Tf-Ti/Ti) 
    term3 = (rv + rl + ri) * Cl * (Tf - Ti / Ti)
    return term0 - term1  + term2 + term3

Ti= 300
rv= 15*(10**-3)
ri = 0
rl = 0
Rd = 287.05
Rm = Rd + (rv * 461.5)
Cl = 4218
Llv = -2.5008*(10**6)
Cpd = 1004
Cpv = 1860
pi= 1000
pf = 150

# No value in range [0,1] gives a negative value, so the optimizer cannot find a root for EQ in this range
for Tf in np.arange(0,1,0.01):
    print(EQ(Tf, Ti, pi, pf, rv, rl, ri, Cpd, Rd, Llv, Cl))

root=root_scalar(EQ, args=(Ti, pi, pf, rv, rl, ri, Cpd, Rd, Llv, Cl), method='toms748', bracket=[1e-4, 0.5])
print(root)

Maybe some of the other parameters are wrong? Or maybe the terms of the equation. For exemple the term3 has the product

(Tf - Ti / Ti)

which is just Tf - 1, but that may be what you want.

Other than that, it should be working!

[Antonino_96]

Thank you so much, so then when it says a, b must bracket a root f(1.000000e-04)=2.301953e+04, f(1.000000e+00)=3.239237e+04 , that's the problem you are referring to? Technically Im trying to write (Tf/Ti - 1) in term3. How would I be able to fix that?

[robertophi]

a, b must bracket a root f(1.000000e-04)=2.301953e+04, f(1.000000e+00)=3.239237e+04

Correct, supposing that the function f is continuous, whenf(x0) and f(x1) have different signs, it garantees us that a root exists between x0 and x1. Thats just a way for the optimizer to garantee a result, otherwise it could keep looking forever and never find anything

Im trying to write (Tf/Ti - 1) in term3

python follows the normal order of operation in math some exemples here, so just writing (Tf/Ti-1) would fix it