0

我正在尝试解决这个非线性方程组：

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

在事先不知道的情况下，我该怎么做才能获得好的解决方案？

4

1 回答 1

1

试试这个，它会循环 3 个范围的 ini，调用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
于 2021-10-29T11:21:37.370 回答