python中的这个欧拉方法代码有什么问题？

计算科学 Python 一体化

2021-12-20 18:00:01

我正在测试我的集成方法重现分析获得的结果（对于参数的各种值）的能力，但结果证明是完全错误的，我找不到我的代码有什么问题。

例如，对于下面的参数值， $x_1$ 应该倾向于 $40$ 随着时间（t）。然而，当我在 PyCharm 中运行这段代码时， $x_1$ 停在 $0.1227$ .

使用的动态方程是

\frac{d x_{1}}{d t} = (N - x_{1}) μ - \frac{θ x_{1}}{k + x_{1}^{2}}

$\frac{dx_{1}}{dt} = (N - x_1) \mu - \frac{\theta x_1}{k + x_1^2}$

任何帮助将不胜感激，谢谢！

代码：

#Parameters:
N = 50  
theta= 0.6  
µ = 0.6/412.5  
k = 1  
n = 2  
x1 = 0   #x1(0)  
listx1 = []  #vector of all x1-values  
dt = 0.01    #delta t  

from math import pow  
import matplotlib.pyplot as plt  

t= 3600*24  # time is 1 day  
while t>0:  
    listx1.append(x1)  
    dx1 = (µ * (N-x1) - (theta * x1) / (k + pow(x1, 2))) * dt  
    x1 = x1 + dx1  
    t -= dt  
    
for i in range(500):
    print("x1=", listx1[i])

print("x1 final:", listx1[-1])

t = []
for i in range(len(listx1)): 
    t.extend([i])

plt.plot(t, listx1, 'r')
plt.autoscale(enable=True, axis='both', tight=False)
v = [0, 3600 * 24, 0, 60] 
plt.axis(v)
plt.show()

1个回答

我没有检查您的代码，但是，您得到的结果也由scipy.integrate.odeint

from scipy.integrate import odeint

def ode(x, t, N, theta, mu, k):
    dxdt = (N - x) * mu - theta * x / (k + x**2)

    return dxdt

N = 50  
theta= 0.6  
mu = 0.6/412.5  
k = 1  

x0 = 0
t = np.linspace(0, 30, 101)
sol = odeint(ode, x0, t, args=(N, theta, mu, k))

这是解决方案的图（请注意，对于 $t=30$ )

稳态值为sol[-1]，等于 $0.1227$ .

但是，您说您期望的值约为 $40$ . 让我们看看发生了什么。在稳定状态下，它必须保持（与 $x_{ss}\triangleq \lim_{t\rightarrow \infty} x_1(t)$ )

0 = (N - x_{s s}) μ - \frac{θ x_{s s}}{k + x_{s s}^{2}}

$0 = (N - x_{ss}) \mu - \frac{\theta x_{ss}}{k+x_{ss}^2}$

因此，稳态解是三阶多项式的根，即系统应该有三个可能的稳态值。这些可以使用 sympy 找到，如下所示。

import sympy as sp

t, N, theta, mu, k = sp.symbols('t, N, theta, mu, k', real = True)
x_ss = sp.symbols('x_ss', real = True)
eq = sp.Eq(0, (N - x_ss) * mu - theta * x_ss / (k + x_ss**2))

sol_sym = sp.solve(eq, x_ss)

sol_sym是以符号形式包含三个根的列表。用数值代替它评估的参数，如下所示

[s.subs({N: 50, theta: 0.6, mu: 0.6/412.5, k:1}).n() for s in sol_sym]

> [0.12273605326634 - 2.29e-16*I, 10.2908641151715 + 3.09e-16*I,
> 39.5863998315621 - 7.96e-17*I]

忽略本质上为零的虚部，稳态值确实可以取数值积分得到的值，以及和。 $0.1227$ $39.586$ $10.29$

显然，在初始条件下，系统达到第一个稳态值。您需要考虑不同的初始值才能达到。例如，使用左右达到此值 $x_1(0)=0$ $40$ $x_1(0) = 35$ $t=10000$

其它你可能感兴趣的问题

上一篇如何处理中间计算中的大数字？下一篇矢量场的旋度分量和无散度分量有什么区别？