求解非线性方程组

Question

我正试图求解这个非线性方程组：

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

在事先不知道的情况下，我可以做些什么来获得好的解决方案

Answer 1

试试这个，它在ini的3个范围内循环，调用solve，如果状态为1，我们会返回，因为状态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