numpy 在Python中求解Colebrook(非线性)方程

8hhllhi2  于 12个月前  发布在  Python
关注(0)|答案(4)|浏览(143)

我想在MATLAB中用python what this guy did来做。
我已经安装了anaconda,所以我有numpy和sympy库。到目前为止,我已经尝试了numpy nsolve,但这不起作用。我应该说我是python的新手,而且我知道如何在MATLAB中实现它:P。
The equation
-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)

是不是因为我同时也在解常数?
edit2:所以我想我必须通过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).现在它工作,但我得到零作为一个结果,

ozxc1zmp

ozxc1zmp1#

这个很好用。我在这里做了一些事情。首先,我使用了一个更简单的函数定义,它使用了您已经定义的全局变量。我发现这比把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
8gsdolmq

8gsdolmq2#

你的标签有点不对:你把它标记为sympy,这是一个用于符号计算的库,但是你想用数字来解决它。如果后者是你的实际意图,这里有相关的scipy文档:
http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#root-finding

ftf50wuq

ftf50wuq3#

Scipy与fixed_point是首选的,因为root不会收敛于远处的猜测值,比如@chthonicdaemon %timeit示例中的0.2。

nbysray5

nbysray54#

在上面的答案中,numpy的求解器解决方案对于这个问题来说相当慢。在这里,一个普通的减半方法将比root和fixed_point方法快大约5倍,并且与root和fixed_point方法一样精确。下面的代码在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 = 100

    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

相关问题