Solving for a specific result to an equation with one variable

Question

If I want to solve for a specific result to an equation with one variable available to increase or decrease to solve for the result I want how could I go about doing it?

I was originally suggested the bisection algorithm but after researching that further I found it to only be applicable for square root equations in python.

For example

rho1 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
H = 1
beta_h1 = .0000000001
def beta_h_func(beta_h):
     return beta_h1*min(Nv/(Qh*Nh), 1) + rho1*H*min(Nv/(Qh*Nh), 1)
while beta_h_func(beta_h1) != 1:
    if beta_h_func(beta_h1) < 1:
        beta_h1 += .000000001
        beta_h_func(beta_h1)
    if beta_h_func(beta_h1) > 1:
        beta_h1 -= .000000001
        beta_h_func(beta_h1)

Where beta_h1 is the variable available for increasing or decreasing. But doing this only results in an infinite back and forth going above and below 1 never ending. What could I change the while function to in order to receive a decisive 1 as my result.

Answer 1

This problem is linear and can easily be solved by hand. The analytical solution is beta_h1 = (1 - rho1*H*min(Nv/(Qh*Nh), 1))/min(Nv/(Qh*Nh), 1). In this case, that value is 0.823.

Assuming this is just an example and your actual problem is more complicated, what you have is a root-finding problem. For that you can use scipy.optimize.root to find the solution. Scipy contains other root-finding methods, including a bisection method that will work for this problem (I'm not sure why you said the bisection method only works for square root problems, whatever that means).

from scipy.optimize import root

rho1 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
H = 1
c = min(Nv/(Qh*Nh), 1)


def beta_h_func(beta_h1):
    return c*(beta_h1 + rho1*H) - 1


sol = root(beta_h_func, 1)
print(sol)

Output:

 message: The solution converged.
 success: True
  status: 1
     fun: [ 0.000e+00]
       x: [ 8.230e-01]
    nfev: 4
    fjac: [[-1.000e+00]]
       r: [-1.000e+00]
     qtf: [ 3.861e-13]

Taking that result, your root is given as sol.x[0] (indexing gives you the float rather than the numpy array).

A couple of comments on your code.

• Your while loop will likely run forever since you are working with floating point numbers and are checking equality. When comparing floating point numbers, check if they're close (with some tolerance), not exactly equal.
• The evaluation of the function after updating beta_h1 is useless since nothing happens with those results.
• The argument to your function is beta_h instead of beta_h1.
• For speed, you can precompute min(...) since the values inside do not change throughout your problem (at least not for the code you shared). I did this in my code above.

Answer 2

Assuming your function will actually equal 1 if you substitute in the correct answer,

>>> def beta_h_func(beta_h):
...      return beta_h*min(Nv/(Qh*Nh), 1) + rho1*H*min(Nv/(Qh*Nh), 1)
...
>>> beta_h_func(.823)==1
True

is your loop logic sound? Let's test it on a simpler function: f(x) = x - 1.

f = lambda x: x - 1

def go(step):
    x = 0
    dir = None
    while (y:=f(x)) != 1:
        if y < 1:
            if dir is None:
                dir = 1
            elif dir == -1:
                return 'looping', x, x+step
            x += step
        elif y > 1:
            if dir is None:
                dir = -1
            elif dir == 1:
                return 'looping', x-step, x
            x -= step
    return(x)

This go function will accept a step and return the value of x at which f(x) is 1. It also has a safety feature to stop if it detects that we are looping between values of x (like you are experiencing). Let's test it:

>>> go(1)
2
>>> go(.5)
2.0
>>> go(.25)
2.0

So the logic is sound. Let's try a smaller step...like 0.1

>>> go(.1)
('looping', 1.9000000000000004, 2.0000000000000004)

Although the logic is sound for exact artihmetic, when you are dealing with floating point numbers you are not assured that the sum of steps will bring you to the single value at which the function will attain the desired value. In my example, stepping by 0.1 from 0 brought us to just over 1.9 and just over 2, but never to 2. Your function (by your own admission) suffers from the same fate. (It also suffers by taking almost 10 billion steps to get to the point at which it loops.)

It is worth your time to see how to implement bisection, e.g. see here.

Advanced Topic

You can even use the built-in bisection routine of Python if you can write your function in terms of integer input. For example, say we wanted to find the solution of x**2 - 70. Let x = i/prec then define a function class that has a __getitem__ method to return the value of the function for a given i and also allows you to specify the value of prec:

>>> class f(object):
...   def __init__(self, prec): self.prec = prec
...   def __getitem__(self, i):
...       x = i/self.prec
...       return x**2-70
...
>>> prec = 1000
>>> f(prec)[1*prec]  # 1**2 - 70
-69.0

Now let the bisect routine find the first value of x that gives a value of f that is >= 0:

>>> from bisect import bisect
>>> prec = 10000; bisect(f(prec), 0, 8*prec, 10*prec)/prec  # nearest 1/1000
8.3667
>>> _**2
70.00166888999999
>>> prec = 2; bisect(f(prec), 0, 8*prec, 10*prec)/prec  # nearest half
8.5
>>> _**2
72.25
>>> prec = 2**10; bisect(f(prec), 0, 8*prec, 10*prec)/prec  # nearest 2**-10
8.3671875
>>> _**2
70.00982666015625