使用 python 和绘图的 ODE 求解系统（Chua 电路）

Question

我正在尝试使用 matplotlib 和 scipy 在 Python 中对 Chau's Circuit 建模，这涉及求解常微分方程组。

这已经完成 in matlab, and I simply wanted to attempt the problem in python. The matlab code linked is a little confusing; the code on the left doesn't appear to have much relevance to solving the system of ode's，它描述了 Chua 的电路（第 3 页，等式 (2)(3) 和 (4)），而右边的代码超出了对电路组件建模的范围组件。

我不熟悉 scipy 的 odeint 函数，所以我使用了 scipy cookbook 中的一些示例作为指导。

谁能帮我排除系统故障；为什么我得到的图表看起来像这样：

与看起来像这样的人不同？

我的代码附在下面：

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

def fV_1(V_1, G_a, G_b, V_b):
    if V_1 < -V_b:
        fV_1 = G_b*V_1+(G_b-G_a)*V_b
    elif -V_b <= V_1 and V_1 <=V_b:
        fV_1 = G_a*V_1
    elif V_1 > V_b:
        fV_1 = G_b*V_1+(G_a-G_b)*V_b
    else:
        print "Error!"
    return fV_1

def ChuaDerivatives(state,t):
    #unpack the state vector
    V_1 = state[0]
    V_2 = state[1]
    I_3 = state[2]

    #definition of constant parameters
    L = 0.018 #H, or 18 mH
    C_1 = 0.00000001 #F, or 10 nF
    C_2 = 0.0000001 #F, or 100 nF
    G_a = -0.000757576 #S, or -757.576 uS
    G_b = -0.000409091 #S, or -409.091 uS
    V_b = 1 #V (E)
    G = 0.000550 #S, or 550 uS VARIABLE

    #compute state derivatives
    dV_1dt = (G/C_1)*(V_2-V_1)-(1/C_1)*fV_1(V_1, G_a, G_b, V_b)
    dV_2dt = -(G/C_2)*(V_2-V_1)+(1/C_2)*I_3
    dI_3dt = -(1/L)*V_2

    #return state derivatives
    return dV_1dt, dV_2dt, dI_3dt

#set up time series
state0 = [0.1, 0.1, 0.0001]
t = np.arange(0.0, 53.0, 0.1)

#populate state information
state = odeint(ChuaDerivatives, state0, t)

# do some fancy 3D plotting
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(state[:,0],state[:,1],state[:,2])
ax.set_xlabel('V_1')
ax.set_ylabel('V_2')
ax.set_zlabel('I_3')
plt.show()

Answer 1

所以我经过一些摆弄后设法自己解决了这个问题；我错误地解释了 odeint 函数；更仔细地阅读文档字符串并从头开始阻止我按照一个困难的方法解决它。代码如下：

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

#define universal variables
c0 = 15.6
c1 = 1.0
c2 = 28.0
m0 = -1.143
m1 = -0.714

#just a little extra, quite unimportant
def f(x):
    f = m1*x+(m0-m1)/2.0*(abs(x+1.0)-abs(x-1.0))
    return f

#the actual function calculating
def dH_dt(H, t=0):
    return np.array([c0*(H[1]-H[0]-f(H[0])),
                  c1*(H[0]-H[1]+H[2]),
                  -c2*H[1]])



#computational time steps
t = np.linspace(0, 30, 1000)
#x, y, and z initial conditions
H0 = [0.7, 0.0, 0.0]

H, infodict = integrate.odeint(dH_dt, H0, t, full_output=True)

print infodict['message']

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(H[:,0], H[:,1], H[:,2])
plt.show()

这给了我这个：

使用 python 和绘图的 ODE 求解系统（Chua 电路）

Solving system of ODE's using python and graphing (Chua's circuit)

python

numpy

matplotlib

ode