решение дифференциального уравнения с функцией шага - PullRequest
Новая
       100

решение дифференциального уравнения с функцией шага

0 голосов
/ 22 мая 2018

Я пытаюсь решить это дифференциальное уравнение как часть моего задания.Я не могу понять, как я могу поставить условие для вас в коде.В коде, показанном ниже, я произвольно указал 

u = 5. 


2dx(t)dt=−x(t)+u(t)

5dy(t)dt=−y(t)+x(t)

u=2S(t−5)

x(0)=0

y(0)=0

where S(t−5) is a step function that changes from zero to one at t=5. When it is multiplied by two, it changes from zero to two at that same time, t=5.

def model(x,t,u):
    dxdt = (-x+u)/2
    return dxdt

def model2(y,x,t):
    dydt = -(y+x)/5
    return dydt

x0 = 0
y0 = 0
u = 5
t = np.linspace(0,40)


x = odeint(model,x0,t,args=(u,))
y = odeint(model2,y0,t,args=(u,))
plt.plot(t,x,'r-')
plt.plot(t,y,'b*')
plt.show()

Ответы [ 2 ]

0 голосов
/ 22 мая 2018

У вас есть пара вопросов здесь, и пошаговая функция - это только малая часть.Вы можете определить пошаговую функцию с помощью простого lambda, а затем просто захватить ее из внешней области, даже не передавая ее в вашу функцию.Потому что иногда это не так, мы будем явными и передадим это.Ваша следующая проблема - порядок аргументов в интегрируемой функции.Согласно документам (у, т, ...).Т.е. сначала функция, затем вектор времени, затем остальные args аргументы.Итак, для первой части мы получаем: 

u = lambda t : 2 if t>5 else 0

def model(x,t,u):
    dxdt = (-x+u(t))/2
    return dxdt

x0 = 0
y0 = 0
t = np.linspace(0,40)


x = odeint(model,x0,t,args=(u,))

Переходя к следующей части, проблема в том, что вы не можете подать x в качестве аргумента к y, потому что это вектор значений для x(t)для определенного времени и поэтому y+x не имеет смысла в функции, как вы ее написали.Вы можете следовать своей интуиции из математического класса, если передадите функцию x вместо значений x.Для этого необходимо, чтобы вы интерполировали значения x, используя конкретные значения времени, которые вас интересуют (с которыми может справиться Сципи, без проблем): 

from scipy.interpolate import interp1d
xfunc = interp1d(t.flatten(),x.flatten(),fill_value="extrapolate") 
#flatten cuz the shape is off , extrapolate because odeint will go out of bounds
def model2(y,t,x):
    dydt = -(y+x(t))/5
    return dydt
y = odeint(model2,y0,t,args=(xfunc,))

Тогда вы получите:

enter image description here

@ Ответ Свена более идиоматичен для векторного программирования, например, scipy / numpy.Но я надеюсь, что мой ответ обеспечит более четкий путь от того, что вы уже знаете, к рабочему решению.

0 голосов
/ 22 мая 2018

Я не очень хорошо знаю библиотеку SciPy, но что касается примера в документации , я бы попробовал что-то вроде этого: 

def model(x, t, K, PT)
    """
    The model consists of the state x in R^2, the time in R and the two
    parameters K and PT regarding the input u as step function, where K
    is the infimum of u and PT is the delay of the step.
    """
    x1, x2 = x   # Split the state into two variables

    u = K if t>=PT else 0    # This is the system input

    # Here comes the differential equation in vectorized form
    dx = [(-x1 + u)/2,
          (-x2 + x1)/5]
    return dx

x0 = [0, 0]
K  = 2
PT = 5
t = np.linspace(0,40)

x = odeint(model, x0, t, args=(K, PT))
plt.plot(t, x[:, 0], 'r-')
plt.plot(t, x[:, 1], 'b*')
plt.show()
...