英文:
solve a system of nonlinear equations using scipy fsolve (math domain error encountered)
问题
我尝试使用Scipy的fsolve来找到一个包含两个非线性方程的系统的答案。 这两个方程是：
f1 = math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**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。以下是我的代码。
from scipy.optimize import fsolve
import math
# 一些全局变量
m = 46.663
chi = 1.1500799949128826
def binodal_fsolve():
def equations(p):
x, y = p
out = []
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))
out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
return out
initial_guess = [0.3, 0.99]
ans = fsolve(equations, initial_guess)
return ans
def test_answers(phiL, phiR):
def functions(x):
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
return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
print(test_answers(0.2542983070, 0.9999999274))
# (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
res = binodal_fsolve()
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:
f1 = math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**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.
from scipy.optimize import fsolve
import math
# some global variables
m = 46.663
chi = 1.1500799949128826
def binodal_fsolve():
def equations(p):
x, y = p
out = []
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))
out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
return out
initial_guess = [0.3, 0.99]
ans = fsolve(equations, initial_guess)
return ans
def test_answers(phiL, phiR):
def functions(x):
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
return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
print (test_answers(0.2542983070, 0.9999999274))
# (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
res = binodal_fsolve()
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算法的修改来以最小二乘的方式解决非线性方程组。有关更多方法，请参阅文档。
from scipy.optimize import root
import numpy.lib.scimath as math
# 一些全局变量
m = 46.663
chi = 1.1500799949128826
def binodal_fsolve():
def equations(p):
x, y = p
out = []
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))
out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
return out
initial_guess = [0.3, 0.99]
#ans = fsolve(equations, initial_guess)
ans = root(equations, initial_guess, method='lm')
return ans
def test_answers(phiL, phiR):
def functions(x):
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
return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
print (test_answers(0.2542983070, 0.9999999274))
# (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
res = binodal_fsolve()
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.
from scipy.optimize import root
import numpy.lib.scimath as math
# some global variables
m = 46.663
chi = 1.1500799949128826
def binodal_fsolve():
def equations(p):
x, y = p
out = []
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))
out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))
return out
initial_guess = [0.3, 0.99]
#ans = fsolve(equations, initial_guess)
ans = root(equations, initial_guess, method='lm')
return ans
def test_answers(phiL, phiR):
def functions(x):
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
return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]
print (test_answers(0.2542983070, 0.9999999274))
# (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
res = binodal_fsolve()
print (res)
Which gives the following roots x and y: : array([0.25429812, 0.99999993]).
The full output:
(1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
/home/user/.local/lib/python3.6/site-packages/scipy/optimize/minpack.py:401: ComplexWarning: Casting complex values to real discards the imaginary part
gtol, maxfev, epsfcn, factor, diag)
cov_x: array([[6.49303571e-01, 8.37627537e-07],
[8.37627537e-07, 1.08484856e-12]])
fjac: array([[ 1.52933340e+07, -1.00000000e+00],
[-1.97290115e+01, -1.24101235e+00]])
fun: array([-2.22945317e-07, -7.20367503e-04])
ipvt: array([2, 1], dtype=int32)
message: 'The relative error between two consecutive iterates is at most 0.000000'
nfev: 84
qtf: array([-0.00338589, 0.00022828])
status: 2
success: True
x: array([0.25429812, 0.99999993])
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论