为什么sympy nsolve给我的值不同于图形显示的值？

如何解决为什么sympy nsolve给我的值不同于图形显示的值？

我正在尝试找到解决微分方程组的参数值（在费率为零时查找，并且值不再更改）。我用ODEINT迭代了系统，并以图形方式显示了它，以查看每个ODE最终收敛到的值。

奇怪的是，NSOLVE给我的值不同于图形显示为最终稳定点的值。

from scipy.optimize import fsolve
import math
import numpy as np
import sympy as sy
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import scipy.signal
from numpy import linalg as LA



a = [0.9,0.9,0.9] #maximum uptake rate
b = [2,2,2] #half-saturation (concentration of nutrients at half maximum rate)
m = [0.2,0.2,0.2] #mortality rate
I = 0.7 #nutrient input
E = 0.3 #nutrient output
r = [0.3] #maximum recycling rate
l = [0] #deforestation rate (biomass loss)
q = [1] #Hill coefficient - ease of decomposition of the matter. This can vary through the season.
s = [0.9] #half-saturation (concentration of dead organic matter (delta*R) at half maximum rate)


"""
setting up a range for iterations of the system for each moment t in time
"""
val = 1000
y0 = 1,1,1
trange = np.linspace(0,val,val) #equally spaced elements


def system(y,t,a,b,m,I,E,r,l,q,s):
    
    N,R,H,P = y
    dydt = [I - E*N + (r[0]*(m[0]*R)**q[0])/(s[0]**q[0]+(m[0]*R)**q[0]) - a[0]*N*R/(b[0] + N),a[0]*N*R/(b[0] + N) - m[0]*R - l[0]*R - a[1]*H*R/(b[1] + R),a[1]*H*R/(b[1] + R) - m[1]*H - a[1]*H*P/(b[2] + H),a[1]*H*P/(b[2] + H) - m[2]*P]
    
    return dydt



solved = odeint(system,y0,trange,args = (a,s),atol = 1.49012e-10)



N = solved[val-1,0]
R = solved[val-1,1]
H = solved[val-1,2]
P = solved[val-1,3]


print(N,P)

返回1.5119062213983654 0.7448846250807435 0.5724768913289174 0.1295697618640645，这很有意义！这是经过时间迭代的ODE：

N,P = sy.symbols("N R H P")




equilibrium_values = (sy.nsolve([I - E*N + (r[0]*(m[0]*R)**q[0])/(s[0]**q[0]+(m[0]*R)**q[0]) - a[0]*N*R/(b[0] + N),a[1]*H*P/(b[2] + H) - m[2]*P],(N,P),(1,1)) )

print(equilibrium_values)

和balance_values返回：

Matrix（[[1.66676383218065]，[0.571428571428325]，[0.597439764807860]，[-1.76967713603266e-13]]）

不是以图形方式找到的平衡值，这应该是求解方程组的结果...有人知道为什么吗？我的Python是否有问题，还是我误解了其背后的数学原理？

解决方法

非线性系统因其在使用数值求解器时对初始猜测的敏感性而臭名昭著。在这种情况下，您没有做错任何事情。如果您使用更接近您在屏幕上看到的初始猜测，您可能会得到 P 不为零的解决方案：

guess N,R,H,P= (1.5,0.7,0.6,0.1)
N = 1.50905265512538
R = 0.749587544253425
H = 0.571428571428571
P = 0.129589598408824

guess N,P= (1.51,0.1)
N = 1.66676383218043
R = 0.571428571428571
H = 0.597439764807826
P = -1.21722045507780e-18

在这种情况下，能够对输入变量进行系统调整是很好的。为此，可以使用如下所示的 tweak 之类的东西：

def tweak(v,p=.1,d=None):
    """Given a vector of initial guesses,return
    vectors with values that are modified by a factor of +/-p
    along with the unmodified values until all possibilities
    are given. d controls number of decimals shown in returned
    values; all are shown when d is None (default).

    EXAMPLES
    ========

    >>> for i in tweak((1,2),.1):
    ...     print(i)
    [1.1,2.2]
    [1.1,2]
    [1.1,1.8]
    [1,2.2]
    [1,2]
    [1,1.8]
    [0.9,2.2]
    [0.9,2]
    [0.9,1.8]
    """"
    from sympy.utilities.iterables import cartes
    if d is None:
        d = 15
    c = []
    for i in v:
        c.append([i*(1+p),i,i*(1-p)])
    for i in cartes(*c):
        yield [round(i,d) for i in i]

eqs = [
    I - E*N + (r[0]*(m[0]*R)**q[0])/(s[0]**q[0]+(m[0]*R)**q[0]) - a[0]*N*R/(b[0] + N),a[0]*N*R/(b[0] + N) - m[0]*R - l[0]*R - a[1]*H*R/(b[1] + R),a[1]*H*R/(b[1] + R) - m[1]*H - a[1]*H*P/(b[2] + H),a[1]*H*P/(b[2] + H) - m[2]*P]

saw = set()
vv = 1.5,0.1
for v0 in tweak(vv,d=3):
    equilibrium_values = nsolve(eqs,(N,P),v0)
    ok = round(equilibrium_values[-1],3)
    if ok not in saw:
        saw.add(ok)
        print('guess N,P=',v0,ok)

给出

guess N,P= [1.65,0.77,0.66,0.11] 0.130
guess N,0.1] 0.0