我在python中用数值方法求解非线性方程组时遇到了一个问题



我目前正试图用Python解决我的流体力学问题,该问题导致了三个不同的方程。经过一些研究,我试图使它尽可能简单,这产生了以下代码:

from scipy.optimize import fsolve 
import math
def equations(p):
d, Re, f = p
return (d**5 - 0.0165*f, 97200 - Re*d, f**(-0.5) + 2*math.log10(0.00015/(3.7*d) + 2.51/(Re*f**(0.5))))

d, Re, f =  fsolve(equations, (0.22, 0.0213, 490000))
print(equations((d, Re, f)))

这给出了一个长错误:

> ValueError                                Traceback (most recent call
> last) <ipython-input-19-dcc7fc844bb6> in <module>()
>       6     return (d**5 - 0.0165*f, 97200 - Re*d, f**(-0.5) + 2*math.log10(0.00015/(3.7*d) + 2.51/(Re*f**(0.5))))
>       7 
> ----> 8 d, Re, f =  fsolve(equations, (0.22, 0.0213, 490000))
>       9 
>      10 print(equations((d, Re, f)))
> 
> ~Anaconda3libsite-packagesscipyoptimizeminpack.py in
> fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev,
> band, epsfcn, factor, diag)
>     146                'diag': diag}
>     147 
> --> 148     res = _root_hybr(func, x0, args, jac=fprime, **options)
>     149     if full_output:
>     150         x = res['x']
> 
> ~Anaconda3libsite-packagesscipyoptimizeminpack.py in
> _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
>     225         with _MINPACK_LOCK:
>     226             retval = _minpack._hybrd(func, x0, args, 1, xtol, maxfev,
> --> 227                                      ml, mu, epsfcn, factor, diag)
>     228     else:
>     229         _check_func('fsolve', 'fprime', Dfun, x0, args, n, (n, n))
> 
> <ipython-input-19-dcc7fc844bb6> in equations(p)
>       4 def equations(p):
>       5     d, Re, f = p
> ----> 6     return (d**5 - 0.0165*f, 97200 - Re*d, f**(-0.5) + 2*math.log10(0.00015/(3.7*d) + 2.51/(Re*f**(0.5))))
>       7 
>       8 d, Re, f =  fsolve(equations, (0.22, 0.0213, 490000))
> 
> ValueError: math domain error

现在,我知道它一定与对数值有关,但我认为方程的复杂性也可能在这里发挥作用。

您可以在equations(p)函数中添加print(d, Re, f),并查看值在每次迭代中的变化情况。

我对它进行了测试,这些是错误之前的结果

0.22 0.0213 490000.0
0.22 0.0213 490000.0
0.22 0.0213 490000.0
0.22000000327825547 0.0213 490000.0
0.22 0.02130000031739473 490000.0
0.22 0.0213 490000.007301569
4384243.368454826 -464.7214227491007 3096116.988633934

正如您所看到的,在计算时,最后一行值实际上对f**(-0.5) + 2*math.log10(0.00015/(3.7*d) + 2.51/(Re*f**(0.5)))操作无效,因为0.00015/(3.7*d) + 2.51/(Re*f**(0.5))的计算结果为-3.069524039032724e-06,接近0,因此被识别为0。因此,由于log(0)本身无效,它返回Domain Error

要解决此问题,您可能需要重新访问您的公式或扩展有效位数字的数量,以容纳更多数字并使其不再为0。

最新更新