求解非线性方程组
Solve a non-linear system of equations
我正在尝试求解这个非线性方程组:
from scipy.optimize import fsolve
import numpy as np
def equations(vars):
x, y, z = vars
eq1 = x - 0.095*np.log(z) - 1.2022
eq2 = y - 100000/x
eq3 = z - y
return [eq1, eq2, eq3]
ini = (1, 1, 1)
x, y, z = fsolve(equations, ini)
print(x, y, z)
系统给了我解决方案但不是问题的解决方案。求解器给了我这个解决方案:x=2.22015, y=14373.01967, z=14372.9181
但真正的解决方案 是 x=2.220157, y=45041.83986, z=45041.83986
.
看来问题出在值的初始化上。如果我将这些值用于初始化:
ini = (2, 40000, 40000)
x, y, z = fsolve(equations, in)
系统给出了真正的解决方案:x=2.220157, y=45041.83986, z=45041.83986
如何在事先不知道的情况下获得好的解决方案?
试试这个,它循环遍历 ini 的 3 个范围,调用 solve,如果状态为 1,我们 return 因为状态 1 是成功或通过状态。我们在 fsolve() 中将 full_output 参数设置为 true 以获取状态信息。
代码
import time
from scipy.optimize import fsolve
import numpy as np
def equations(vars):
x, y, z = vars
eq1 = x - 0.095*np.log(z) - 1.2022
eq2 = y - 100000/x
eq3 = z - y
return [eq1, eq2, eq3]
def sol():
ret = None
for i in range(1, 1000):
print(f'processing ... {i}')
for j in range(1, 1000):
for k in range(1, 1000):
ini = (i, j, k)
res = fsolve(equations, ini, full_output=True)
if res[2] == 1: # if status is 1 then it is solved.
return res
ret = res
return ret
# Test
t1 = time.perf_counter()
res = sol()
print(f'status: {res[2]}, sol: {res[0]}')
print(f'elapse: {time.perf_counter() - t1:0.1f}s')
输出
processing ... 1
status: 1, sol: [2.22015798e+00 4.50418399e+04 4.50418399e+04]
elapse: 2.9s
我正在尝试求解这个非线性方程组:
from scipy.optimize import fsolve
import numpy as np
def equations(vars):
x, y, z = vars
eq1 = x - 0.095*np.log(z) - 1.2022
eq2 = y - 100000/x
eq3 = z - y
return [eq1, eq2, eq3]
ini = (1, 1, 1)
x, y, z = fsolve(equations, ini)
print(x, y, z)
系统给了我解决方案但不是问题的解决方案。求解器给了我这个解决方案:x=2.22015, y=14373.01967, z=14372.9181
但真正的解决方案 是 x=2.220157, y=45041.83986, z=45041.83986
.
看来问题出在值的初始化上。如果我将这些值用于初始化:
ini = (2, 40000, 40000)
x, y, z = fsolve(equations, in)
系统给出了真正的解决方案:x=2.220157, y=45041.83986, z=45041.83986
如何在事先不知道的情况下获得好的解决方案?
试试这个,它循环遍历 ini 的 3 个范围,调用 solve,如果状态为 1,我们 return 因为状态 1 是成功或通过状态。我们在 fsolve() 中将 full_output 参数设置为 true 以获取状态信息。
代码
import time
from scipy.optimize import fsolve
import numpy as np
def equations(vars):
x, y, z = vars
eq1 = x - 0.095*np.log(z) - 1.2022
eq2 = y - 100000/x
eq3 = z - y
return [eq1, eq2, eq3]
def sol():
ret = None
for i in range(1, 1000):
print(f'processing ... {i}')
for j in range(1, 1000):
for k in range(1, 1000):
ini = (i, j, k)
res = fsolve(equations, ini, full_output=True)
if res[2] == 1: # if status is 1 then it is solved.
return res
ret = res
return ret
# Test
t1 = time.perf_counter()
res = sol()
print(f'status: {res[2]}, sol: {res[0]}')
print(f'elapse: {time.perf_counter() - t1:0.1f}s')
输出
processing ... 1
status: 1, sol: [2.22015798e+00 4.50418399e+04 4.50418399e+04]
elapse: 2.9s