Как использовать оператор if в дифференциальном уравнении (SciPy)? - PullRequest
Новая
       96

Как использовать оператор if в дифференциальном уравнении (SciPy)?

0 голосов
/ 05 июня 2018

Я пытаюсь решить дифференциальное уравнение с помощью Python.В этом двухсистемном дифференциальном уравнении, если значение первой переменной (v) превышает пороговое значение (30), его следует сбросить до другого значения (-65).Ниже я положил свой код.Проблема в том, что значение первой переменной после достижения 30 остается постоянным и не сбрасывается до -65.Эти уравнения описывают динамику одного нейрона.Уравнения взяты с этого веб-сайта и этого PDF-файла . 

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
from scipy.integrate import odeint
plt.close('all')

a = 0.02
b = 0.2
c = -65
d = 8
i = 0

p = [a,b,c,d,i]

def fun(u,tspan,*p):
    du = [0,0]
    if u[0] < 30: #Checking if the threshold has been reached
        du[0] = (0.04*u[0] + 5)*u[0] + 150 - u[1] - p[4]
        du[1] = p[0]*(p[1]*u[0]-u[1])
    else:
        u[0] = p[2] #reset to -65    
        u[1] = u[1] + p[3] 

    return du


p = tuple(p)

y0 = [0,0]

tspan = np.linspace(0,100,1000)
sol = odeint(fun, y0, tspan, args=p)

 fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)         
plt.plot(tspan,sol[:,0],'k',linewidth = 5)
plt.plot(tspan,sol[:,1],'r',linewidth = 5)
myleg = plt.legend(['v','u'],\
    loc='upper right',prop = {'size':28,'weight':'bold'}, bbox_to_anchor=(1,0.9))

Решение выглядит следующим образом:

enter image description here

Вот правильное решение по Julia, здесь u1 представляет v: enter image description here

Это код Julia: 

using DifferentialEquations
using Plots

a = 0.02
b = 0.2
c = -65
d = 8
i = 0

p = [a,b,c,d,i]

function fun(du,u,p,t)
    if u[1] <30
        du[1] = (0.04*u[1] + 5)*u[1] + 150 - u[2] - p[5]
        du[2] = p[1]*(p[2]*u[1]-u[2])
    else
        u[1] = p[3]
        u[2] = u[2] + p[4]
    end
end

u0 = [0.0;0.0]
tspan = (0.0,100)
prob = ODEProblem(fun,u0,tspan,p)
tic()
sol = solve(prob,reltol = 1e-8)
toc()

plot(sol)

1 Ответ

0 голосов
/ 05 июня 2018

Рекомендуемое решение

Использует события и интегрируется отдельно после каждого разрыва. 

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

a = 0.02
b = 0.2
c = -65
d = 8
i = 0

p = [a,b,c,d,i]

# Define event function and make it a terminal event
def event(t, u):
    return u[0] - 30
event.terminal = True

# Define differential equation
def fun(t, u):
    du = [(0.04*u[0] + 5)*u[0] + 150 - u[1] - p[4],
          p[0]*(p[1]*u[0]-u[1])]
    return du

u = [0,0]

ts = []
ys = []
t = 0
tend = 100
while True:
    sol = solve_ivp(fun, (t, tend), u, events=event)
    ts.append(sol.t)
    ys.append(sol.y)
    if sol.status == 1: # Event was hit
        # New start time for integration
        t = sol.t[-1]
        # Reset initial state
        u = sol.y[:, -1].copy()
        u[0] = p[2] #reset to -65    
        u[1] = u[1] + p[3]
    else:
        break

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
# We have to stitch together the separate simulation results for plotting
ax.plot(np.concatenate(ts), np.concatenate(ys, axis=1).T)
myleg = plt.legend(['v','u'])

Минимальное изменение "решение"

Кажется, что ваш подход работает просто отличнос solve_ivp.

Предупреждение Я думаю, что и у Юлии, и у solve_ivp правильный способ обработки такого рода вещей - использовать события.Я полагаю, что приведенный ниже подход основан на деталях реализации, то есть вектор состояния, передаваемый функции, является тем же объектом, что и вектор внутреннего состояния, что позволяет нам изменять его на месте.Если бы это была копия, этот подход не сработал бы.Кроме того, в этом подходе нет гарантии, что решатель предпринимает достаточно малые шаги, чтобы наступить правильная точка, где достигнут предел.Использование событий сделает это более правильным и обобщенным для других дифференциальных уравнений, которые, возможно, имеют более низкие градиенты до разрыва. 

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FormatStrFormatter
from scipy.integrate import solve_ivp
plt.close('all')

a = 0.02
b = 0.2
c = -65
d = 8
i = 0

p = [a,b,c,d,i]

def fun(t, u):
    du = [0,0]
    if u[0] < 30: #Checking if the threshold has been reached
        du[0] = (0.04*u[0] + 5)*u[0] + 150 - u[1] - p[4]
        du[1] = p[0]*(p[1]*u[0]-u[1])
    else:
        u[0] = p[2] #reset to -65    
        u[1] = u[1] + p[3] 

    return du

y0 = [0,0]

tspan = (0,100)
sol = solve_ivp(fun, tspan, y0)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)         
plt.plot(sol.t,sol.y[0, :],'k',linewidth = 5)
plt.plot(sol.t,sol.y[1, :],'r',linewidth = 5)
myleg = plt.legend(['v','u'],loc='upper right',prop = {'size':28,'weight':'bold'}, bbox_to_anchor=(1,0.9))

Результат

enter image description here

...