问绘制微分方程系统函数外的变量
EN

Stack Overflow用户

提问于 2019-01-08 10:15:43

回答 1查看 109关注 0票数 1

我有一个函数中的4-4微分方程组(subsystem4)，我用odeint函数解决了它。我设法画出了系统的结果。我的问题是我想要绘制和一些其他方程(例如x，y，vcxdot...)它们都包含在相同的函数(subsystem4)中，但我得到了NameError:没有定义名称'vcxdot‘。此外，我想使用其中的一些方程(不仅是方程系统的结果)作为下面的微分方程系统的输入，并绘制同一时间段(t)内的所有方程。我已经使用Matlab-Simulink做到了这一点，但由于Simulink块的存在，它要容易得多。如何访问和绘制函数(subsystem4)的所有方程式，并将它们用作以下系统的输入？我是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.7

plot

differential-equations

回答 1

Stack Overflow用户

发布于 2019-01-08 12:45:42

对于直接绘制导数图的问题，您可以通过直接在解决方案上再次调用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

这可能会使计算看起来加倍，但并不是真的，因为解点通常是从求解器的内部步长点插值的。在积分过程中，常微分方程函数的求值通常发生在与解点不同的点上。

对于其他变量，将位置和速度的计算提取到函数中，以便仅在一个位置具有常量和组合：

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)