用Python中的牛顿拉夫森方法求解colebrook（非线性）方程

Question

我曾尝试解决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的变化值处，并针对Re绘制f。

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。问题是您编写代码的方式只能针对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-但这将是另一个问题。如果您需要修复代码的帮助，建议使用physics标记询问一个新问题，以便获得有用的结果，而不是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()