为什么 sympy nsolve 给我的值与图形显示的值不同? [英] Why is sympy nsolve giving me different values than what is graphically shown?

问题描述

我试图找出求解微分方程组的参数值(找出比率为零且值不再变化的时间).我使用 ODEINT 迭代系统并以图形方式显示它以查看每个 ODE 最终收敛到的值.

I am trying to find the values of parameters that solve a system of differential equations (find when the rates are zero, and the values no longer change). I iterated the system with ODEINT and graphically shown it to see what are the values that each ODE eventually converges to.

奇怪的是，NSOLVE 给我的值与图表显示的最终稳定点不同.

The weird part is, NSOLVE is giving me different values than what the graph shows as the final stable point.

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, 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, b, m, I, E, r, l, q, 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, R, H, P)

返回 1.5119062213983654 0.7448846250807435 0.5724768913289174 0.1295697618640645，这是有道理的！这是通过时间迭代的 ODE:

returns 1.5119062213983654 0.7448846250807435 0.5724768913289174 0.1295697618640645, which makes sense! Here's the ODE iterated through time:

N, R, H, 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[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],(N, R, H, P),(1,1, 1, 1)) )

print(equilibrium_values)

和 equilibrium_values 返回:

and equilibrium_values returns:

矩阵([[1.66676383218065], [0.571428571428325], [0.597439764807860], [-1.76967713603266e-13]])

Matrix([[1.66676383218065], [0.571428571428325], [0.597439764807860], [-1.76967713603266e-13]])

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

Which are not the graphically found equilibrium values, which should be the result of solving the system of equations... Anyone has a clue why? Is it an issue with my Python, or have I misunderstood the math behind it?

推荐答案

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

Nonlinear systems are notorious for their sensitivity to initial guesses when using a numeric solver. In this case you have done nothing wrong. If you use initial guesses closer to what you see on the screen you might get a solution where P does not go to zero:

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,R,H,P= (1.51, 0.7, 0.6, 0.1)
N = 1.66676383218043
R = 0.571428571428571
H = 0.597439764807826
P = -1.21722045507780e-18

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

In cases like these it is nice to be able to make a systematic tweaking of the input variables. For that, something like tweak shown below can be used:

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.7, 0.6, 0.1
for v0 in tweak(vv, d=3):
    equilibrium_values = nsolve(eqs, (N, R, H, P), v0)
    ok = round(equilibrium_values[-1], 3)
    if ok not in saw:
        saw.add(ok)
        print('guess N,R,H,P=',v0,ok)

给出

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

这篇关于为什么 sympy nsolve 给我的值与图形显示的值不同?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！