Я пытаюсь визуализировать систему осцилляторов в Python. Пока в одной из переменных нет накопления (с использованием функции Дирака-дельта), система работает. Однако, когда я пытаюсь включить его, в системе не появляется никаких изменений. Ниже вы можете найти процесс, который я прошел, вместе с некоторыми вопросами, которые у меня есть.
Я написал некоторый код на python, который может правильно визуализировать поведение обычного осциллятора. Однако при включении пошагового приращения в одну из переменных система не реагирует на это. Ниже вы можете найти версию с включенным приращением. Код выполняется, но выдает результат, как будто приращения там не было. По коду вы можете найти мои проблемы с правильным решением.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
a4 = 5.68
aAPD = 9.09
deltaZ = 0.41
gamma = 1.33
negV = 1
n=500
def f(V):
if V<=-60:
return V/20 + 4
elif V>=20:
return V/20 - 1
else:
return -V/80 + 1/4
def alpha(V):
if V<0:
return a4
else:
return -aAPD
def g(q):
return q/(q+1)
def step(V):
global negV
if V>0 and negV == 1:
negV = 0
return deltaZ
elif V<=0 and negV == 0:
negV = 1
return 0
else:
return 0
# function that returns dz/dt
def model(z,t):
V = z[0]
y = z[1]
q = z[2]
dVdt = 25000*(y - f(V))
dydt = alpha(V) - a4*g(q)
dqdt = -gamma*g(q) + step(V)
dzdt = [dVdt,dydt,dqdt]
return dzdt
# initial condition
z0 = [0,0,0]
# time points
t = np.linspace(0,0.5,n)
# solve ODE
z = odeint(model,z0,t)
# plot results
plt.plot(t,z[:,0],label='Voltage')
plt.plot(t,z[:,1],label='y')
plt.plot(t,z[:,2],label='Z')
plt.ylabel('response')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()
Проблема заключается в функции step(V). Осциллятор пересекает в определенный момент значение V = 0. В этот момент значение q должно быть увеличено на deltaZ.
В одесольвере python есть несколько моментов, которые я не понимаю или не знаю, как получить доступ к правильным значениям. Поэтому вот краткий список с проблемами / вопросами:
- Чтобы найти точку, где V пересекает 0, я использую глобальное значение negV. Тем не менее, было бы лучше, если бы я мог включить заявление в соответствии с
if V[i]>0 and V[i-1]<0: return deltaZ. Эта информация должна быть доступна где-то, так как в конце вы строите значения V, но я не понимаю, как получить доступ к предыдущим значениям V.
- приращение значения q должно происходить мгновенно. Поэтому одним из решений может быть сначала решить dqdt, а затем увеличить q на значение deltaZ. Однако я не вижу, как это можно сделать. Есть ли способ прямого доступа к значениям V, y и q? Позже я также хотел бы изменить значение V вручную, так что это было бы полезно.
- другое решение может заключаться в том, чтобы сохранить значение deltaZ в dqdt, но тогда я должен знать временной шаг, с которым решаются дифференциальные уравнения. Таким образом, я могу разделить deltaZ = 0.41 на этот временной шаг, так что приращение q равно deltaZ = 0,41. При просмотре информации odeint я увидел, что временной шаг является переменным, так что, возможно, это не так уж хорошо для решения ...
Желательно, чтобы я разрешил все три пункта, но если бы вы смогли дать ответ на один из них, мне бы тоже очень помогли.