Solving colebrook (nonlinear) equation with Newton Raphson method in python

Question

I have tried solving the colebrook (nonlinear) equation for frictional factor in python but I keep getting this error:

Traceback (most recent call last):
  File "c:/Users/WWW/Desktop/WWWWWWWW/WWWWWWWWWWW/DDD/WWWWW.py", line 46, in <module>
    F_factor = Newton(f=0.1)
  File "c:/Users/WWW/Desktop/WWWWWWWW/WWWWWWWWWWW/DDD/WWWWW.py", line 22, in Newton
    eps_new = float(func(f))/dydf(f)
TypeError: only size-1 arrays can be converted to Python scalars

I know something is wrong but not sure what is.

I am trying to find the friction factor (f) for this equation

-0.86 * log(2.51 / (Re * sqrt(f)) + e / D / 3.7) = 1 / sqrt(f)

at varying values of Reynold's number Re and plotting f against Re.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import time

#Parameters
e_D = 1e-4
Re = np.linspace(10**4,10**7,10)
eps=1e-7

def func(f):
    return -0.86*np.log((e_D/3.7)+(2.51/Re)*f**0.5)-f**0.5

def dydf(f):
    return (((-1.255/Re)*f**-1.5)/((e_D/3.7)+((2.51/Re)*f**-0.5)))-((f**-1.5)/1.72)

def Newton(f,conv_hist=True):
    eps_new = float(func(f))/dydf(f)
    #y_value = func(f)
    iteration_counter = 0
    history = []
    while abs(eps_new) >= eps and iteration_counter < 100:
        eps_new = float(func(f))/dydf(f)
        f = f - eps_new
        #y_value = func(f)
        iteration_counter += 1
        history.append([iteration_counter, f, eps_new])

    # Here, either a solution is found, or too many iterations
        if abs(dydf(f)) <= eps:
            print('derivative near zero!, dydf =', dydf)
            print(func(f), 'iter# =', iteration_counter, 'eps =', eps_new)
            break

    if conv_hist:
            hist_dataframe = pd.DataFrame(history, columns=['Iteration #', 'f', 'eps'])
            print(hist_dataframe.style.hide_index())

    return f, iteration_counter
startTime = time.time()
F_factor = Newton(f=0.1)
endTime = time.time()
print(F_factor)

print('Total process took %f seconds!' % (endTime - startTime))
plt.loglog(F_factor, Re, marker='o')
plt.title('f vs Re')
plt.grid(b=True, which='minor')
plt.grid(b=True, which='major')
plt.xlabel('Re')
plt.ylabel('f')
plt.savefig('fvsRe.png')
plt.show()

Can you include the whole error traceback? – William Miller Jan 20 '20 at 05:38 — William Miller, Jan 20 '20 at 05:38
I just added the whole error traceback. Thanks! – nebula1 Jan 20 '20 at 07:08 — nebula1, Jan 20 '20 at 07:08

score 0 · Accepted Answer · answered Jan 20 '20 at 21:51

The function func returns an array object, so

float(func(f))

is not valid Python. The issue is that the way you have written the code you can only use a specific Reynold's value for each call to Newton i.e.

def func(f, re):
    return -0.86*np.log((e_D/3.7)+(2.51/re)*f**0.5)-f**0.5

def dydf(f, re):
    return (((-1.255/re)*f**-1.5)/((e_D/3.7)+((2.51/re)*f**-0.5)))-((f**-1.5)/1.72)

def Newton(f, re, conv_hist=True):
    eps_new = func(f, re)/dydf(f, re)
    # ...
    while abs(eps_new) >= eps and iteration_counter < 100:
        eps_new = func(f, re)/dydf(f, re)
        # ...
        if abs(dydf(f, re)) <= eps:
            # ...

Will run without the error, but it will return nan for every Reynold's number - but that would be the subject of another question. I would recommend asking a new question with the physics tag if you need help fixing the code so it gets a useful result instead of nan.

Complete runnable code:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import time

#Parameters
e_D = 1e-4
Re = np.linspace(10**4,10**7,10)
eps=1e-7

def func(f, re):
    return -0.86*np.log((e_D/3.7)+(2.51/re)*f**0.5)-f**0.5

def dydf(f, re):
    return (-1.255/re*f**-1.5) / (e_D/3.7 + 2.51/re*f**-0.5) - (f**-1.5) / 1.72

def Newton(f, re, conv_hist=True):
    eps_new = func(f, re)/dydf(f, re)
    iteration_counter = 0
    history = []
    while abs(eps_new) >= eps and iteration_counter < 100:
        eps_new = func(f, re)/dydf(f, re)
        f = f - eps_new
        iteration_counter += 1
        history.append([iteration_counter, f, eps_new])

        if abs(dydf(f, re)) <= eps:
            print('derivative near zero!, dydf =', dydf)
            print(func(f), 'iter# =', iteration_counter, 'eps =', eps_new)
            break

    if conv_hist:
            hist_dataframe = pd.DataFrame(history, columns=['Iteration #', 'f', 'eps'])

    return f, iteration_counter
startTime = time.time()
F_factors = []
for re in Re:
    F_factors.append(Newton(0.1, re)[0])
endTime = time.time()

print('Total process took %f seconds!' % (endTime - startTime))
plt.loglog(F_factors, Re, marker='o')
plt.title('f vs Re')
plt.grid(b=True, which='minor')
plt.grid(b=True, which='major')
plt.xlabel('Re')
plt.ylabel('f')
plt.savefig('fvsRe.png')
plt.show()

Solving colebrook (nonlinear) equation with Newton Raphson method in python

1 Answers1