求解带阶跃函数的微分方程

Question

我正在尝试解决这个微分方程，作为我作业的一部分。我无法理解如何将您的条件放入代码中。下图代码中，我随意提供了

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()

Answer 1

我不太了解 SciPy 库，但关于 example in the documentation 我会尝试这样的事情：

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()

Answer 2

这里有几个问题，step 函数只是其中的一小部分。您可以使用简单的lambda 定义一个步进函数，然后简单地从外部范围捕获它，甚至无需将其传递给您的函数。因为有时情况并非如此，我们会明确并通过它。 您的下一个问题是要集成的函数中参数的顺序。根据docs (y,t,...)。即，首先是函数，然后是时间向量，然后是其他 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 作为 arg 提供给 y，因为它是特定时间 x(t) 的值向量，因此 y+x 在你写的函数。如果您传递 x 函数而不是 x 值，则可以遵循数学课的直觉。这样做需要您使用您感兴趣的特定时间值插入 x 值（scipy 可以处理，没问题）：

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,))

然后你得到：

@Sven 的答案对于像 scipy/numpy 这样的向量编程来说更惯用。但我希望我的回答提供了一条从您已经知道的到可行的解决方案的更清晰的路径。