在Python中求解colebrook（非线性）方程

Question

我想用 python 做这家伙在 MATLAB 中做了什么。

我已经安装了 anaconda，所以我有 numpy 和 sympy 库。到目前为止，我已经尝试使用 numpy nsolve，但这不起作用。我应该说我是 python 新手，而且我知道如何在 MATLAB 中执行此操作：P。

等式：

-2*log(( 2.51/(331428*sqrt(x)) ) + ( 0.0002 /(3.71*0.26)) ) = 1/sqrt(x)

通常，我会迭代地解决这个问题，只需猜测左侧的 x，然后求解右侧的 x。将解放在左边，再次求解。重复直到左 x 接近右。我知道应该采取什么解决方案。

所以我可以这样做，但这不是很酷。我想用数字来做。 我的15欧元卡西欧计算器可以解决这个问题，所以我认为它应该不会太复杂？

谢谢您的帮助，

编辑：所以我尝试了以下方法：

from scipy.optimize import brentq

w=10;
d=0.22;
rho=1.18;
ni=18.2e-6;

Re=(w*d*rho)/ni
k=0.2e-3;
d=0.26;

def f(x,Re,k,d):
    return (
        -2*log((2.51/(Re*sqrt(x)))+(k/(3.71*d)),10)*sqrt(x)+1
            );

print(
    scipy.optimize.brentq
        (
        f,0.0,1.0,xtol=4.44e-12,maxiter=100,args=(),full_output=True,disp=True
        )
    );

我得到了这个结果：

    r = _zeros._brentq(f,a,b,xtol,maxiter,args,full_output,disp)
TypeError: f() takes exactly 4 arguments (1 given)

是因为我在求解的同时也在求解常数吗？

编辑2： 所以我想我必须通过 args=() 关键字分配常量，所以我改变了：

f,0.0,1.0,xtol=4.44e-12,maxiter=100,args=(Re,k,d),full_output=True,disp=True

但现在我明白了：

-2*log((2.51/(Re*sqrt(x)))+(k/(3.71*d)),10)*sqrt(x)+1
TypeError: return arrays must be of ArrayType

无论如何，当我输入不同的方程时；可以说

2*x*Re+(k*d)/(x+5)

它有效，所以我想我必须改变方程。

所以它死在这里：log(x,10)..

edit4：正确的语法是 log10(x)...现在它可以工作，但结果我得到零

Answer 1

这个效果很好。我在这里做了一些事情。首先，我使用您定义的全局变量对函数进行了更简单的定义。我发现这比将 args= 传递给求解器要好一些，如果您需要类似的东西，它还可以更轻松地使用您自己的自定义求解器。我使用通用

root

函数作为入口点，而不是使用特定的算法 - 这很好，因为您可以稍后简单地传递不同的方法。我还按照 PEP 8 的建议修复了您的间距，并修复了您对方程式的错误重写。我发现简单地写 LHS - RHS 比像你那样操作更直观。另外，请注意，我已将所有整数文字替换为 1.0 或其他值，以避免整数除法出现问题。 0.02 被认为是摩擦系数的一个相当标准的起点。

import numpy
from scipy.optimize import root

w = 10.0
d = 0.22
rho = 1.18
ni = 18.2e-6

Re = w*d*rho/ni
k = 0.2e-3

def f(x):
    return (-2*numpy.log10((2.51/(Re*numpy.sqrt(x))) + (k/(3.71*d))) - 1.0/numpy.sqrt(x))

print root(f, 0.02)

我还必须提到，对于这个问题，定点迭代实际上比牛顿法还要快。您可以通过定义

f2

来使用内置定点迭代例程，如下所示：

def f2(x):
    LHS = -2*numpy.log10((2.51/(Re*numpy.sqrt(x))) + (k/(3.71*d)))
    return 1/LHS**2

时序（从根开始更远以显示收敛速度）：

%timeit root(f, 0.2)
1000 loops, best of 3: 428 µs per loop

%timeit fixed_point(f2, 0.2)
10000 loops, best of 3: 148 µs per loop

Answer 2

您的标签有点偏离：您将其标记为

sympy

，这是一个用于 symbolic 计算的库，但说您想以数字方式解决它。如果后者是您的实际意图，这里是相关的 scipy 文档：

http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#root-finding

Answer 3

带有

fixed_point

的 Scipy 是首选，因为

root

对于远处的猜测值不会收敛，就像 @chthonicdaemon %timeit 示例中的 0.2。

Answer 4

上面答案中 numpy 的求解器解决方案对于这个问题来说相当慢。正常的减半方法至少比 root 和fixed_point 方法快 5 倍且准确。下面的代码在factor()函数中有简单的减半方法。 Factor()函数和fixed_point方法都是定时的。

import numpy
import time

from scipy.optimize import fixed_point

w = 10.0
d = 0.5
rho = 1.18
ni = 18.2e-6

Re = w*d*rho/ni
k = 0.2e-3

def f2(x):
    LHS = -2*numpy.log10((2.51/(Re*numpy.sqrt(x))) + (k/(3.71*d)))
    return 1/LHS**2

def f(x):
    return (-2*numpy.log10((2.51/(Re*numpy.sqrt(x))) + (k/(3.71*d))) - 1.0/numpy.sqrt(x))

def factor():
    min_value = 0.001
    max_value = 0.2
    tolerance = 0.000001
    iterations = 0
    max_iterations = 1000

    while iterations < max_iterations:
        factor = (min_value + max_value) / 2
        result = f(factor)
        if abs(result) > tolerance:
            iterations += 1
            if result < 0:
                min_value = factor
            else:
                max_value = factor
        else:
            return factor
    return 0

tic = time.perf_counter()
a = factor()
toc = time.perf_counter()
print(f"*** Solution reached in {toc - tic:0.6f} seconds")
print(f"Factor found: {a}") 

tic = time.perf_counter()
a = fixed_point(f2, 0.2)
toc = time.perf_counter()
print(f"*** Solution reached in {toc - tic:0.6f} seconds")
print(f"Factor found: {a}")

*** 0.000098 秒内得出解决方案 找到的因子：0.01750585120916367 *** 0.000429 秒内得出解决方案 找到的因子：0.017505852375130332