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

发布时间：2020-12-20 13:40:32 所属栏目：Python 来源：网络整理

导读：我试图解决dB的以下等式(为简单起见,我在问题标题中将dB表示为x)：等式中的所有其他项都是已知的.我尝试使用SymPy来象征性地解决dB,但我一直在节省时间.我也尝试过使用scipy.optimize的fminbound,但dB的答案是错误的(参见下面的使用fminbound方法的Python代

我试图解决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

解决方法

以下是四个单调的根方法：

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

还有newton-raphson(割线/哈利)方法：

from scipy.optimize import newton
dbub = newton(dB,dbed)
print 'dbub = ',dB(dbub)

如果你有一个包围间隔,scipy文档建议使用brentq.

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!