Python Newton Method - Blowing up

Question

I am running into some issues with Newton method. I tried to debug but dont know what seems to be the root cause. Thought I will post it out here to see if somebody can help.

I have a simple function to compute oil density, rho_p0, which takes an initial guess and then computes rho_a and then iterates to make the final value for rho_p0.

from scipy.optimize import newton
def mccain_hill_rhop0(rho_p0, Rs,Sg,So):
        try:
            a0 = -49.8930
            a1 = 85.0149
            a2 = - 3.70373
            a3 = 0.0479818
            a4 = 2.98914
            a5 = - 0.0356888
            print(rho_p0)
            rho_a=(a0+a1*Sg+a2*Sg*rho_p0+a3*Sg*rho_p0**2+a4*rho_p0+a5*rho_p0**2)
            rho_p0=(Rs*Sg+4600*So)/(73.71+Rs*Sg/rho_a)
            return rho_p0
        except RuntimeError:
            return None
Rs=1000
Sg=0.55
So=0.8
rho_p01=52.8-0.01*Rs
rho_p0= newton(mccain_hill_rhop0, rho_p01, args=(Rs,Sg,So,), tol=10**(-3) , maxiter=100)
print(rho_p0)

Here is the error message. I dont know why it blows after the second iteration. Any idea? Thanks for your help

RuntimeError: Tolerance of 2.5470966514398777e+33 reached. Failed to converge after 7 iterations, value is 2.5470966514398786e+33.
42.8
42.80438
-452.7973894722485
1393.5747708318825
672642.4222763113
1505488387.2917664
8.483536060177454e+17
2.5470966514398786e+33

No idea what you are trying to compute, but you do know that Newton's Method is not guaranteed to converge, right? — AirSquid, Jun 01 '22 at 18:43
You might be getting underflow in rho_a: Try rearranging the main "return" value so you're multiplying by rho_a rather than dividing by a number with it in the denominator. Also try printing the value to be returned and along-side rho_po and the inputs: `print((rho_p0, rho_a, mccain, Rs, Sg, So))`. Note the doubled parens. — Sarah Messer, Jun 01 '22 at 18:44
I am trying to compute rho_p0, by providing an initial value rho_p01 to the function — Dark Knight, Jun 01 '22 at 18:45
what are you trying to compute? Newton's method is designed to find 0 of a function, also, this function has some discontinuities that would make method crash. Maybe you should search another method with restrictions — Ulises Bussi, Jun 01 '22 at 18:51
I did rearrange the return part to return rho_a*(Rs*Sg+4600*So)/(73.71*rho_a+Rs*Sg) . and printed print((rho_p0, rho_a, Rs, Sg, So)) It still blows up . (42.8, 20.580650689599977, 1000, 0.55, 0.8) (42.80438, 20.581334274602213, 1000, 0.55, 0.8) (-452.79738946006813, -2340.730937905859, 1000, 0.55, 0.8) — Dark Knight, Jun 01 '22 at 18:54
@UlisesBussi I have been using newton optimization instead of writing the while loops until the difference between iterations are within tolerance on many cases. This one seems to crash... — Dark Knight, Jun 01 '22 at 19:00
I don't understand the new return, but, the first problem is: it's not clear which one is the rho_p0 value. Is the value that makes the function image 0? is the maximun value ? — Ulises Bussi, Jun 01 '22 at 19:06
@UlisesBussi the computation is like this. rho_p0 = (Rs*Sg+4600*So)/(73.71+Rs*Sg/rho_a) , where rho_a is a f(rho_p0) . So I am trying to iterate to solve for rho_p0. — Dark Knight, Jun 01 '22 at 21:39
@SarahMesser .. I rearranged the equation to keep rho_a in the numerator and it work... Thanks. — Dark Knight, Jun 01 '22 at 21:50

Dark Knight · Answer 1 · 2022-06-01T22:31:22.320

I rearranged the return function to have the variable I am solving for in the numerator and it worked.

from scipy.optimize import newton
def mccain_hill_rhop0(rho_p0, Rs,Sg,So):
        try:
            a0 = -49.8930
            a1 = 85.0149
            a2 = - 3.70373
            a3 = 0.0479818
            a4 = 2.98914
            a5 = - 0.0356888
            rho_a=(a0+a1*Sg+a2*Sg*rho_p0+a3*Sg*rho_p0**2+a4*rho_p0+a5*rho_p0**2)
            print((rho_p0, rho_a))
            rho_p0=rho_a*(Rs*Sg+4600*So-73.71*rho_p0)/(Rs*Sg)
            return rho_p0
        except RuntimeError:
            return None

output

(42.8, 20.580650689599977)
(42.80438, 20.581334274602213)
(59.200734064755764, 20.63974342124446)
(57.3824271698717, 20.879752959108018)
(57.3871110909038, 20.879213691644452)
57.38705738842851

Python Newton Method - Blowing up

1 Answers1