在 Python 中求解 x 的高度非线性方程

Question

我正在尝试求解以下 dB 等式（为简单起见，我在问题标题中将 dB 表示为 x）：

等式中的所有其他项都是已知的。我尝试使用 SymPy 来象征性地求解 dB，但我不断收到超时错误。我也尝试使用scipy.optimize 中的fminbound，但是dB 的答案是错误的（请参阅下面的使用fminbound 方法的Python 代码）。

有人知道用 Python 求解 dB 方程的方法吗？

import numpy as np

from scipy.optimize import fminbound

#------------------------------------------------------------------------------
# parameters

umf = 0.063         # minimum fluidization velocity, m/s
dbed = 0.055        # bed diameter, m
z0 = 0              # position bubbles are generated, m
z = 0.117           # bed vertical position, m
g = 9.81            # gravity, m/s^2

#------------------------------------------------------------------------------
# calculations

m = 3                       # multiplier for Umf
u = m*umf                   # gas superficial velocity, m/s

abed = (np.pi*dbed**2)/4.0  # bed cross-sectional area, m^2

# calculate parameters used in equation

dbmax = 2.59*(g**-0.2)*(abed*(u-umf))**0.4
dbmin = 3.77*(u-umf)**2/g

c1 = 2.56*10**-2*((dbed / g)**0.5/umf)

c2 = (c1**2 + (4*dbmax)/dbed)**0.5

c3 = 0.25*dbed*(c1 + c2)**2

dbeq = 0.25*dbed*(-c1 + (c1**2 + 4*(dbmax/dbed))**0.5 )**2

# general form of equation ... (term1)^power1 * (term2)^power2 = term3

power1 = 1 - c1/c2

power2 = 1 + c1/c2

term3 = np.exp(-0.3*(z - z0)/dbed)

def dB(d):
    term1 = (np.sqrt(d) - np.sqrt(dbeq)) / (np.sqrt(dbmin) - np.sqrt(dbeq))
    term2 = (np.sqrt(d) + np.sqrt(c3)) / (np.sqrt(dbmin) + np.sqrt(c3))
    return term1**power1 * term2**power2 - term3

# solve main equation for dB

dbub = fminbound(dB, 0.01, dbed)

print 'dbub = ', dbub

Answer 1

解决标题中的问题很容易：

In [9]:

import numpy as np

import scipy.optimize as so

In [10]:

def f(x):

    return ((x-0.32)**0.8+(x+1.45)**1.1-np.exp(0.8))**2

In [11]:

so.fmin(f, x0=5)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 40

Out[11]:

array([ 0.45172119])

In [12]:

f(0.45172119)

Out[12]:

4.7663411535618792e-13

其他参数都固定了吗？

Answer 2

以下是四种单维度根方法：

from scipy.optimize import brentq, brenth, ridder, bisect
for rootMth in [brentq, brenth, ridder, bisect]:
    dbub = rootMth(dB, 0.01, dbed)
    print 'dbub = ', dbub, '; sanity check (is it a root?):', dB(dbub)

还有 newton-raphson (secant / haley) 方法：

from scipy.optimize import newton
dbub = newton(dB, dbed)
print 'dbub = ', dbub, '; sanity check (is it a root?):', dB(dbub)

如果您有括号间隔，scipy 文档建议使用 brentq。