Float division by 0 error in optimize.newton

Question

I'm a beginning programmer, and I'm trying to find the root of an equation by using Newton's method. The value keeps returning 0 however, which results in a 'ZeroDivisionError: float division by zero' in 'self.β_func'. I'm trying to find the x-value in which p_max_T = p_max, which is why I've subtracted them from each other in beta_func.

It concerns the following piece of code:

import math
import scipy.optimize

P_1 = 1.22E-3
P_2 = 1.3E-3
T_P1 = 285
T_P2 = 288
T_AP = (1 / T_P1 - 1 / T_P2)**(-1) * math.log(P_2 / P_1)
T_APH = 25924
T_APL = 27774
T_PL = 271
T_PH = 296
I_sat = 200
temp = 13 #Example of a temperature in .csv file
α = 3.75e-5


p_max_T = (P_1 * math.exp(T_AP / T_P1 - T_AP / (temp + 273.15))) / (1 + math.exp(T_APL / (temp + 273.15) - T_APL / T_PL) + math.exp(T_APH / T_PH - T_APH / (temp + 273.15)))

#p_max = (α * I_sat / math.log(1 + α / x)) * (α / (α + x)) * (x / (α + x))**(x / α) 

beta_func = lambda x: p_max_T - (α * I_sat / math.log(1 + α / x)) * (α / (α + x)) * (x / (α + x))**(x / α) 

beta_test = scipy.optimize.newton(beta_func, 1E-9)

For some reason, optimize.newton finds x=0, which results in a 'division by 0' error.

The initial guess of x is 1E-9. I have tried solving the equation for every temp that is in my .csv file, but none of them should result in x=0.

Any advice or insights would be very welcome.

Thanks!

UPDATE

---

I've tried applying the bisection method. However, it states that no root is present.

def f(x):
    return((P_1 * math.exp(T_AP / T_P1 - T_AP / (temp + 273.15))) / (1 + math.exp(T_APL / (temp + 273.15) - T_APL / T_PL) + math.exp(T_APH / T_PH - T_APH / (temp + 273.15)))-(α * I_sat / math.log(1 + α / temp)) * (α / (α + temp)) * (temp / (α + temp))**(temp / α))

def bisection_method(a, b, tol):
    if f(a)*f(b) > 0:
        #end function, no root.
        print("No root found.")
    else:
        while (b - a)/2.0 > tol:
            midpoint = (a + b)/2.0
            if f(midpoint) == 0:
                return(midpoint) 
            elif f(a)*f(midpoint) < 0: 
                b = midpoint
            else:
                a = midpoint
        return(midpoint)

answer = bisection_method(0, 1, 1E-10)
print("Answer:", answer)

Answer 1

So, I have no idea what equations you are using, but let me show you a couple things that might help you diagnose your issue here, especially if you are new to programming and working with Newton's method or similar things.

A couple pointers:

1. Make some functions that can be repeatedly called with values. I've shown you some ideas below. You can add parameters to the functions if other items are variables. Right now these rely on "global" variables from the stuff you define above, but that's ok for a small program. You can then use these to send in test values to ensure they are working and producing expected results. You can take what I started below and make a "beta function" and use that later

2. Break up super long equations into components that can be individually inspected. See what I did with your first equation. I could inspect a and b separately inside the function, if needed. That and it is much more readable.

3. Make some plots to see if things make sense. Using pyplot is pretty easy and there are a ton of examples out there. There are several other plotting packages that are also very straightforward. So the basic recipe is: Make a function, make a vector (list) of "x" values, use a list comprehension to make a corresponding list of "y" values using your function, pass the x and y vectors to the plotter. I did that for 2 of your functions (below). I suspect there is something wrong with your p_max formula. It looks really odd.

You are working with very small numbers right near zero and the slopes are not clear yet, so Newton may be sending things toward zero causing your error. Tinker with your functions by sending them test values with assertion statements or such until you know they are correct, then proceed.

Code:

import math
import scipy.optimize
from matplotlib import pyplot as plt
import numpy as np

P_1 = 1.22E-3
P_2 = 1.3E-3
T_P1 = 285
T_P2 = 288
T_AP = (1 / T_P1 - 1 / T_P2)**(-1) * math.log(P_2 / P_1)
T_APH = 25924
T_APL = 27774
T_PL = 271
T_PH = 296
I_sat = 200
temp = 13 #Example of a temperature in .csv file
α = 3.75e-5

def p_max_T(temp):
    a = (P_1 * math.exp(T_AP / T_P1 - T_AP / (temp + 273.15)))
    b = (1 + math.exp(T_APL / (temp + 273.15) - T_APL / T_PL) + math.exp(T_APH / T_PH - T_APH / (temp + 273.15)))
    result = a/b
    return result

def p_max(x):
    result = (α * I_sat / math.log(1 + α / x)) * (α / (α + x)) * (x / (α + x))**(x / α)
    return result



temps = np.arange(start=0, stop=5, step=1e-4)
pmaxt = [p_max_T(temp) for temp in temps]

x_vals = np.arange(start=1e-10, stop=2, step=1e-5)
pmax = [p_max(x) for x in x_vals]

plt.plot(temps, pmaxt, label='p max t')
plt.plot(x_vals, pmax, label='p max')
plt.xlabel('t or x')
plt.legend()
plt.show()

Plot

Update (so this is turning into a long post... :))

Perhaps your functions are correct... IDK. I made a beta function and plotted it (for fixed value of T=13) in vicinity of your solution. There is a root near there, but Newton will likely fail here due to asymptotic behavior. You might want to consider bisection or another method.

new code segment and plot...

def beta(x, temp=13):
    return p_max_T(temp) - p_max(x)


start = 10e-10
stop = 1e-6
steps = 1000
step = (stop-start)/steps
x_vals = np.arange(start, stop, step)

pmax = [p_max(x) for x in x_vals]
beta_val = [beta(x, temp) for x in x_vals]

# plt.plot(temps, pmaxt, label='p max t')
plt.plot(x_vals, pmax, label='p max')
plt.plot(x_vals, beta_val, label='beta')
plt.axhline(color='k', alpha=0.5, ls='--')
plt.xlabel('x')
plt.title(f'Beta for t={temp} and p_max_T={p_max_T(temp):0.2e}')
plt.legend()
plt.show()

beta_test = scipy.optimize.newton(beta, 1E-8)