用Python语言模拟两个耦合动力系统

2024-10-04 01:27:30 发布

您现在位置:Python中文网/ 问答频道 /正文

目标是绘制两个相同的耦合动力系统。你知道吗

我们有:

X=[x0,x1,x2]

U=[u0,u1,u2]

以及

Xdot=f(X)+α*(U-X)

Udot=f(U)+α*(X-U)

因此,我希望在一组轴上绘制这个大系统的解(例如在xyz中),并最终改变耦合强度来研究其行为。你知道吗

import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D

def couple(s,t,a=0.2,beta=0.2,gamma=5.7,alpha=0.03):
    [x,u] = s
    [u0,u1,u2] = u
    [x0,x1,x2] = x
    xdot = np.zeros(3)
    xdot[0] = -x1-x2
    xdot[1] = x0+a*x1
    xdot[2] = beta + x2*(x0-gamma)
    udot = np.zeros(3)
    udot[0] = -u1-u2
    udot[1] = u0+a*u1
    udot[2] = beta + u2*(u0-gamma)
    sdot = np.zeros(2)
    sdot[0] = xdot + alpha*(u-x)
    sdot[1] = udot + alpha*(x-u)
    return sdot

s_init = [0.1,0.1]

t_init=0; t_final = 300; t_step = 0.01
tpoints = np.arange(t_init,t_final,t_step)

a=0.2; beta=0.2; gamma=5.7; alpha=0.03

y = odeint(couple, s_init, tpoints,args=(a,beta,gamma,alpha), hmax = 0.01)

我认为sïinit有问题,因为它应该是两个初始条件向量,但当我尝试得到“odeint:y0应该是一维的。”另一方面,当我尝试sïinit是一个6向量时,我得到“太多的值要解包(预期是两个)。”在当前的设置中,我得到了错误

  File "C:/Users/Python Scripts/dynsys2019work.py", line 88, in couple
    [u0,u1,u2] = u

TypeError: cannot unpack non-iterable numpy.float64 object

干杯

*编辑:请注意,这基本上是我第一次尝试这种事情,将很高兴收到进一步的文件和参考。你知道吗


Tags: importalphainitnpbetax1x2gamma
1条回答
网友
1楼 · 发布于 2024-10-04 01:27:30

ode定义接受并返回scipy odeint中的1D向量,我认为您的一些困惑是,实际上有一个包含6个变量的ODEs系统。你刚刚在心里把它分成两个独立的颂歌,它们是耦合的。你知道吗

你可以这样做:

import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np

def couple(s,t,a=0.2,beta=0.2,gamma=5.7,alpha=0.03):

    x0, x1, x2, u0, u1, u2 = s
    xdot = np.zeros(3)
    xdot[0] = -x1-x2
    xdot[1] = x0+a*x1
    xdot[2] = beta + x2*(x0-gamma)
    udot = np.zeros(3)
    udot[0] = -u1-u2
    udot[1] = u0+a*u1
    udot[2] = beta + u2*(u0-gamma)
    return np.ravel([xdot, udot])

s_init = [0.1,0.1, 0.1, 0.1, 0.1, 0.1]
t_init=0; t_final = 300; t_step = 0.01
tpoints = np.arange(t_init,t_final,t_step)
a=0.2; beta=0.2; gamma=5.7; alpha=0.03
y = odeint(couple, s_init, tpoints,args=(a,beta,gamma,alpha), hmax = 0.01)
plt.plot(tpoints,y[:,0])

相关问题 更多 >