首页 文章

求解python中的colebrook(非线性)方程

提问于
浏览
3

我想在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)......现在它可以工作,但结果是零

2 回答

  • 0

    这很好用 . 我使用全局变量使用了一个更简单的函数定义,你使用泛型 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
    
  • 1

    你的标签有点偏差:你将它标记为 sympy 这是一个用于符号计算的库,但是说你想用数字来解决它 . 如果后者是你的实际意图,这里是相关的scipy文档:

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

相关问题