在 python 中用 Newton Raphson 方法求解 colebrook(非线性)方程

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

我已经尝试在 python 中求解摩擦系数的 colebrook(非线性)方程,但我一直收到此错误:

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

我知道出了点问题,但不确定是什么。

我正试图找到这个 equation

的摩擦系数 (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()

函数func return是一个数组对象,所以

float(func(f))

无效 Python。问题是您编写代码的方式只能为每次调用 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:
            # ...

将 运行 没有错误,但它将 return 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()