Решение уравнения Колебрука (нелинейного) методом Ньютона Рафсона в python

Question

Я пытался решить уравнение Колебрука (нелинейное) для коэффициента трения в python, но продолжаю получать эту ошибку:

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

Я знаю, что что-то не так, но не уверен, что именно.

Я пытаюсь найти коэффициент трения (f) для этого уравнения

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

при различных значениях числа Рейнольдса Re и построении графика f относительно 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()

Answer 1

Функция func возвращает объект массива, поэтому

float(func(f))

недопустимо Python. Проблема в том, что при написании кода вы можете использовать только определенное c значение Рейнольдса для каждого вызова Newton, т.е.

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:
            # ...

будет работать без ошибки, но вернется nan для каждого числа Рейнольда - но это будет предметом другого вопроса. Я бы порекомендовал задать новый вопрос с тегом физика , если вам нужна помощь в исправлении кода, чтобы он получил полезный результат вместо nan.

Полный исполняемый код:

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()