python迭代法求解非线性方程_在Python中求解x的高度非线性方程

I am trying to solve the following equation for dB (for simplicity, I stated dB as x in the question title):

All of the other terms in the equation are known. I tried using SymPy to symbolically solve for dB but I kept getting time out errors. I also tried using fminbound from scipy.optimize but the answer for dB is wrong (see below for Python code using the fminbound approach).

Does anyone know of a way to solve the equation for dB using 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

解决方案

Here are the four single-dim root-methods:

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)

Also the newton-raphson (secant / haley) method:

from scipy.optimize import newton

dbub = newton(dB, dbed)

print 'dbub = ', dbub, '; sanity check (is it a root?):', dB(dbub)

The scipy documentation recommends brentq if you have a bracketing interval.