使用Scipy的fsolve解决非线性方程组（遇到数学域错误）。

huangapple
go评论105阅读模式
英文:

solve a system of nonlinear equations using scipy fsolve (math domain error encountered)

问题

我尝试使用Scipy的fsolve来找到一个包含两个非线性方程的系统的答案。 这两个方程是：

  1. f1 = math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2
  2. f2 = math.log(1 - x) - (m - 1)*x + chi*m*x**2

在这种情况下，m 和 chi 是常数。基本目标是找到同时满足 f1(x) = f1(y) 和 f2(x) = f2(y) 的 x 和 y。我知道 x 和 y 的初始猜测值为 0.3 和 0.99。以下是我的代码。

  1. from scipy.optimize import fsolve
  2. import math
  4. # 一些全局变量
  5. m = 46.663
  6. chi = 1.1500799949128826
  8. def binodal_fsolve():
  9.     def equations(p):
  10.         x, y = p
  11.         out = []
  12.         out.append(math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2 - (math.log(y) + 1. - ((1. + (m - 1)*y) / m) + chi * (1 - y)**2))
  13.         out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
  15.         return out
  17.     initial_guess = [0.3, 0.99]
  18.     ans = fsolve(equations, initial_guess)
  20.     return ans
  22. def test_answers(phiL, phiR):
  23.     def functions(x):
  24.         return math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2, math.log(1 - x) - (m - 1)*x + chi*m*x**2
  26.     return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
  28. print(test_answers(0.2542983070, 0.9999999274))
  29. # (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
  30. res = binodal_fsolve()
  31. print(res)

当我执行代码时，我总是遇到"math domain error"错误。然而，如果我尝试使用MAPLE的fsolve解决它，我可以得到答案(0.2542983070, 0.9999999274)。

通过将这些值代入方程中，我得到(1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)，这表明答案是正确的。
我不知道如何使scipy的fsolve工作。任何建议将不胜感激。

英文:

I tried to use Scipy's fsolve to find the answers to a system of two nonlinear equations.
The two equations are:

  1. f1 = math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2
  2. f2 = math.log(1 - x) - (m - 1)*x + chi*m*x**2

m and chi are constants in this case. The essential goal is to find x, y that satisfies simultaneously f1(x) = f1(y) and f2(x) = f2(y). I know the initial guess for x, y are 0.3 and 0.99. Below is my code.

  1. from scipy.optimize import fsolve
  2. import math
  4. # some global variables
  5. m = 46.663
  6. chi = 1.1500799949128826
  8. def binodal_fsolve():
  9.     def equations(p):
  10.         x, y = p
  11.         out = []
  12.         out.append(math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2 - (math.log(y) + 1. - ((1. + (m - 1)*y) / m) + chi * (1 - y)**2))
  13.         out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
  15.         return out
  16.     
  17.     initial_guess = [0.3, 0.99]
  18.     ans = fsolve(equations, initial_guess)
  20.     return ans
  22. def test_answers(phiL, phiR):
  23.     def functions(x):
  24.         return math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2, math.log(1 - x) - (m - 1)*x + chi*m*x**2
  26.     return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
  28. print (test_answers(0.2542983070, 0.9999999274))
  29. # (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
  30. res = binodal_fsolve()
  31. print (res)

When I executed the code, I always encountered the math domain error.
However, if I tried to solve it using MAPLE fsolve. I can get the answers (0.2542983070, 0.9999999274).

By plugging these back to the equations, I get (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684) which suggests the answers are correct.
I don't know how to make scipy fsolve work. Any suggestions will be greatly appreciated.

答案1

得分: 1

在这种情况下，您可以使用numpy.lib.scimath中的log函数，当其参数为负数时返回一个复数。

不要使用scipy.optimize.fsolve，而是使用scipy.optimize.root，并将方法更改为lm，该方法使用Levenberg-Marquardt算法的修改来以最小二乘的方式解决非线性方程组。有关更多方法，请参阅文档。

  1. from scipy.optimize import root
  3. import numpy.lib.scimath as math
  5. # 一些全局变量
  6. m = 46.663
  7. chi = 1.1500799949128826
  9. def binodal_fsolve():
  10.     def equations(p):
  11.         x, y = p
  12.         out = []
  13.         out.append(math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2 - (math.log(y) + 1. - ((1. + (m - 1)*y) / m) + chi * (1 - y)**2))
  14.         out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
  16.         return out
  18.     initial_guess = [0.3, 0.99]
  19.     #ans = fsolve(equations, initial_guess)
  20.     ans = root(equations, initial_guess, method='lm')
  22.     return ans
  24. def test_answers(phiL, phiR):
  25.     def functions(x):
  26.         return math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2, math.log(1 - x) - (m - 1)*x + chi*m*x**2
  28.     return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
  30. print (test_answers(0.2542983070, 0.9999999274))
  31. # (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
  32. res = binodal_fsolve()
  33. print (res)

这将给出以下根：x和y的值：array([0.25429812, 0.99999993])。

英文:

In this case you can use the log function from numpy.lib.scimath that returns a complex number when its argument is negative.

Instead of using scipy.optimize.fsolve, use scipy.optimize.root and change the method to lm which solves the system of nonlinear equations in a least squares sense using a modification of the Levenberg-Marquardt algorithm. For more methods, see the documentation.

  1. from scipy.optimize import root
  3. import numpy.lib.scimath as math
  5. # some global variables
  6. m = 46.663
  7. chi = 1.1500799949128826
  9. def binodal_fsolve():
  10.     def equations(p):
  11.         x, y = p
  12.         out = []
  13.         out.append(math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2 - (math.log(y) + 1. - ((1. + (m - 1)*y) / m) + chi * (1 - y)**2))
  14.         out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
  16.         return out
  18.     initial_guess = [0.3, 0.99]
  19.     #ans = fsolve(equations, initial_guess)
  20.     ans = root(equations, initial_guess, method='lm')
  22.     return ans
  24. def test_answers(phiL, phiR):
  25.     def functions(x):
  26.         return math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2, math.log(1 - x) - (m - 1)*x + chi*m*x**2
  28.     return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
  30. print (test_answers(0.2542983070, 0.9999999274))
  31. # (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
  32. res = binodal_fsolve()
  33. print (res)

Which gives the following roots x and y: : array([0.25429812, 0.99999993]).

The full output:

  1. (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
  2. /home/user/.local/lib/python3.6/site-packages/scipy/optimize/minpack.py:401: ComplexWarning: Casting complex values to real discards the imaginary part
  3.   gtol, maxfev, epsfcn, factor, diag)
  4.    cov_x: array([[6.49303571e-01, 8.37627537e-07],
  5.        [8.37627537e-07, 1.08484856e-12]])
  6.     fjac: array([[ 1.52933340e+07, -1.00000000e+00],
  7.        [-1.97290115e+01, -1.24101235e+00]])
  8.      fun: array([-2.22945317e-07, -7.20367503e-04])
  9.     ipvt: array([2, 1], dtype=int32)
  10.  message: 'The relative error between two consecutive iterates is at most 0.000000'
  11.     nfev: 84
  12.      qtf: array([-0.00338589,  0.00022828])
  13.   status: 2
  14.  success: True
  15.        x: array([0.25429812, 0.99999993])

huangapple
  • 本文由 发表于 2020年1月7日 01:46:33
  • 转载请务必保留本文链接：https://go.coder-hub.com/59616645.html
  • optimization
  • scipy
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定