ZeroDivisionError: деление с плавающей точкой на ноль (решение нелинейного уравнения Колебрука методом Ньютона-Рафсона в python) - PullRequest
Новая
       66

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

1 голос
/ 22 января 2020

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

ZeroDivisionError: деление с плавающей точкой на ноль

вот полный возврат : 

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

Я пытаюсь найти коэффициент трения (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()

1 Ответ

2 голосов
/ 22 января 2020

Ваша проблема в этой строке: 

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)

Я думаю, что вместо этого вы будете sh: 

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

Однако это не сработает, потому что вы sh сюжет f против Re. Поэтому вместо этого вы должны составить список fa и добавить результаты: 

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

Final result

...