python - 零除法错误: float division by zero (Solving colebrook (nonlinear) equation with Newton Raphson method in python)

Question

我尝试在 python 中求解摩擦系数的科尔布鲁克(非线性)方程，但我不断收到此错误:

ZeroDivisionError: float 除以零

这是完整的回溯:

Traceback (most recent call last):
  File "c:/Users/BDG/Desktop/kkk/www/Plots/jjj/Code.py", line 49, in <module>
    f = Newton(f0,re)
  File "c:/Users/BDG/Desktop/kkk/www/Plots/jjj/Code.py", line 20, in Newton
    eps_new = func(f, Re)/dydf(f, Re)
  File "c:/Users/BDG/Desktop/kkk/www/Plots/jjj/Code.py", line 13, in func
    return -0.86*np.log((e_D/3.7)+((2.51/Re))*f**(-0.5))-f**(-0.5)
ZeroDivisionError: float division by zero

我正在尝试找到此 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

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.0793/(Re*((251*f**-0.5)/(100*Re)+(10*e_D)/37)*(f**1.5)))+(1/(2*(f**1.5)))

def Newton(f0, Re, conv_hist=True):
    f = f0
    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, func(f,Re), eps_new])

        if abs(dydf(f, Re)) <= eps:
            print('derivative near zero!, dydf =', dydf(f,re))
            print(dydf(f,re), 'iter# =', iteration_counter, 'eps =', eps_new)
            break
        if iteration_counter == 99:
            print('maximum iterations reached!')
            print(f, 'iter# = ', iteration_counter)
            break 

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

    return f
startTime = time.time()

Re = np.linspace(10**4,10**7,100)
f0 = 0.001
for re in range(len(Re)):
    f = Newton(f0,re)
endTime = time.time()

print('Total process took %f seconds!' % (endTime - startTime))
plt.loglog(Re, f,  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

您的问题在这一行:

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

当Re为0时，此操作失败。发生这种情况的原因是:

for re in range(len(Re)):
    f = Newton(f0,re)

我认为你想做的是:

for re in Re:
    f = Newton(f0,re)

但是，这不起作用，因为您希望绘制 f 与 Re 的图。因此，您应该创建一个列表并附加结果:

f = []
for re in Re:
    f.append(Newton(f0,re))