Решение системы ODE с использованием питона и графического кода (схема Chua)

Question

Я пытаюсь моделировать Chau's Circuit в Python с использованием matplotlib и scipy, что предполагает решение системы обыкновенных дифференциальных уравнений.

Это было сделано in matlab, и я просто хотел попробовать проблему в python. Связанный с Matlab код немного запутан; код слева не имеет большого значения для решения system of ode's, которые описывают схему Chua's (стр. 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() 

Что дает мне это: