时间戳

首页 > 解决方案 > 从微分方程系统函数中绘制变量

python-2.7 - 从微分方程系统函数中绘制变量

问题描述

我在一个函数（subsystem4）中有一个 4-4 微分方程系统，我用 odeint 函数解决了它。我设法绘制了系统的结果。我的问题是我想绘制包含在同一函数（subsystem4）中的其他一些方程（例如 x、y、vcxdot ...），但我得到 NameError: name 'vcxdot' is not defined。此外，我想使用其中一些方程（不仅是方程系统的结果）作为以下微分方程系统的输入，并在同一时间段 (t) 内绘制所有方程。我已经使用 Matlab-Simulink 完成了这项工作，但由于有 Simulink 块，它变得容易得多。如何访问和绘制函数（子系统 4）的所有方程并将它们用作以下系统的输入？我是 python 新手，我使用 Python 2.7.12。先感谢您！

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

def subsystem4(u,t):
    added_mass_x = 0.03 # kg
    added_mass_y = 0.04
    mb = 0.3 # kg
    m1 = mb-added_mass_x
    m2 = mb-added_mass_y
    l1 = 0.07 # m
    l2 = 0.05 # m
    J = 0.00050797 # kgm^2
    Sa = 0.0110 # m^2
    Cd = 2.44
    Cl = 3.41
    Kd = 0.000655 # kgm^2
    r = 1000 # kg/m^3
    f = 2 # Hz
    c1 = 0.5*r*Sa*Cd
    c2 = 0.5*r*Sa*Cl
    c3 = 0.5*mb*(l1**2)
    c4 = Kd/J
    c5 = (1/(2*J))*(l1**2)*mb*l2
    c6 = (1/(3*J))*(l1**3)*mb


    vcx = u[0]
    vcy = u[1]
    psi = u[2]
    wz = u[3]

    x = 3 + 0.3*np.cos(t)
    y = 0.5 + 0.3*np.sin(t)
    xdot = -0.3*np.sin(t)
    ydot = 0.3*np.cos(t)
    xdotdot = -0.3*np.cos(t)
    ydotdot = -0.3*np.sin(t)
    vcx = xdot*np.cos(psi)-ydot*np.sin(psi)
    vcy = ydot*np.cos(psi)+xdot*np.sin(psi)

    psidot = wz
    vcxdot = xdotdot*np.cos(psi)-xdot*np.sin(psi)*psidot-ydotdot*np.sin(psi)-ydot*np.cos(psi)*psidot
    vcydot = ydotdot*np.cos(psi)-ydot*np.sin(psi)*psidot+xdotdot*np.sin(psi)+xdot*np.cos(psi)*psidot
    g1 = -(m1/c3)*vcxdot+(m2/c3)*vcy*wz-(c1/c3)*vcx*np.sqrt((vcx**2)+(vcy**2))+(c2/c3)*vcy*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)
    g2 = (m2/c3)*vcydot+(m1/c3)*vcx*wz+(c1/c3)*vcy*np.sqrt((vcx**2)+(vcy**2))+(c2/c3)*vcx*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)

    A = 12*np.sin(2*np.pi*f*t+np.pi)
    if A>=0.1:
        wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz**2*np.sign(wz)-c5*g2-c6*np.sqrt((g1**2)+(g2**2))
    elif A<-0.1:
        wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz**2*np.sign(wz)-c5*g2+c6*np.sqrt((g1**2)+(g2**2))
    else:
        wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz**2*np.sign(wz)-c5*g2


    return [vcxdot,vcydot,psidot,wzdot]

u0 = [0,0,0,0]
t = np.linspace(0,15,1000)
u = odeint(subsystem4,u0,t)

vcx = u[:,0]
vcy = u[:,1]
psi = u[:,2]
wz = u[:,3]

plt.figure(1)
plt.subplot(211)
plt.plot(t,vcx,'r-',linewidth=2,label='vcx')
plt.plot(t,vcy,'b--',linewidth=2,label='vcy')
plt.plot(t,psi,'g:',linewidth=2,label='psi')
plt.plot(t,wz,'c',linewidth=2,label='wz')
plt.xlabel('time')
plt.legend()
plt.show()

标签： python-2.7plotdifferential-equations

解决方案

对于绘制导数的直接问题，您可以通过在解决方案上再次直接调用 ODE 函数来获得速度，

u = odeint(subsystem4,u0,t)
udot = subsystem4(u.T,t)

并通过获得单独的速度数组

vcxdot,vcydot,psidot,wzdot = udot

在这种情况下，该函数涉及分支，这对它的矢量化调用不是很友好。有一些方法可以矢量化分支，但最简单的解决方法是手动循环通过解决方案点，这比工作矢量化实现要慢。这将再次产生一个像 一样的元组列表odeint，因此必须将结果转换为列表元组，以便“轻松”分配给单个数组变量。

udot = [ subsystem4(uk, tk) for uk, tk in zip(u,t) ]; 
vcxdot,vcydot,psidot,wzdot = np.asarray(udot).T

扩展情节

这似乎会使计算量加倍，但实际上并非如此，因为解点通常是从求解器的内部步长点插值而来的。积分期间 ODE 函数的评估通常发生在与解点不同的点上。

对于其他变量，将位置和速度的计算提取到函数中，以使常数和成分仅在一个位置：

def xy_pos(t): return 3 + 0.3*np.cos(t), 0.5 + 0.3*np.sin(t)
def xy_vel(t): return -0.3*np.sin(t), 0.3*np.cos(t) 
def xy_acc(t): return -0.3*np.cos(t), -0.3*np.sin(t)

或类似的，您可以在 ODE 函数内部和准备绘图时使用。

Simulink 最有可能做的是收集所有模块的所有方程，并将其形成一个大型 ODE 系统，然后一次求解整个状态。您将需要实现类似的东西。一个大的状态向量，每个子系统都知道它的状态片段。导数向量从中获取其特定的状态变量并将导数写入。然后，导数的计算可以使用子系统之间传递的值。

您正在尝试做的事情，分别解决子系统，仅适用于resp。可能会导致 1 阶集成方法。所有高阶方法都需要能够同时在从该方法的先前阶段计算的某个方向上转移状态，并在那里评估整个系统。

推荐阅读