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