In [28]:

    

from sympy import *
init_printing(use_latex=True)
from IPython.display import display
from sympy.printing import print_python, print_fcode


    

In [31]:

    

x0,v0,x1,v1,p0,p1, m,c,k, wn, wd, z, h, t, dp= symbols(
    'x_0,v_0,x_1,v_1,p_0,p_1, m,c,k, omega_n, omega_d, z, h, t, Delta_p')
display((x0,v0,x1,v1,p0,p1, m,c,k, wn, wd, z, h, t, dp))
A,B,C,D=symbols('A,B,C,D')


    

    






$$\begin{pmatrix}x_{0}, & v_{0}, & x_{1}, & v_{1}, & p_{0}, & p_{1}, & m, & c, & k, & \omega_{n}, & \omega_{d}, & z, & h, & t, & \Delta_{p}\end{pmatrix}$$

In [32]:

    

x = exp(-z*wn*t) * (A*cos(wd*t)+B*sin(wd*t))
v = x.diff(t)


    

In [33]:

    

p =(p0*(h-t)+p1*t)/h
p = p0 +dp*t/h
x_p = C*t + D
v_p = x_p.diff(t)
aaa = (k*x_p+c*v_p-p).expand()
CD = solve((aaa.coeff(t,1),aaa.coeff(t,0)), (C,D))
x_p = x_p.subs(CD)
v_p = v_p.subs(CD)
display(CD)


    

    






$$\begin{Bmatrix}C : \frac{\Delta_{p}}{h k}, & D : - \frac{\Delta_{p} c}{h k^{2}} + \frac{p_{0}}{k}\end{Bmatrix}$$

In [37]:

    

eqx = x.subs(t,0)+x_p.subs(t,0)-x0
eqv = v.subs(t,0)+v_p.subs(t,0)-v0
AB = solve((eqx,eqv),(A,B))


    

In [38]:

    

xh = x.subs(AB).subs(t,h)
vh = v.subs(AB).subs(t,h)


    

In [39]:

    

for i, r in enumerate((xh, vh)):
    for j, s in enumerate((x0, v0, p0, dp)):
        print_fcode(r.diff(s).simplify(), 
                    source_format='free',
                    assign_to="A[%d,%d]"%(i,j))
    print


    

    




A[0,0] = (omega_d*cos(h*omega_d) + omega_n*z*sin(h*omega_d))*exp(-h* &
      omega_n*z)/omega_d
A[0,1] = exp(-h*omega_n*z)*sin(h*omega_d)/omega_d
A[0,2] = -(omega_d*cos(h*omega_d) + omega_n*z*sin(h*omega_d))*exp(-h* &
      omega_n*z)/(k*omega_d)
A[0,3] = (c*omega_d*cos(h*omega_d) + (c*omega_n*z - k)*sin(h*omega_d))* &
      exp(-h*omega_n*z)/(h*k**2*omega_d)

A[1,0] = -(omega_d**2 + omega_n**2*z**2)*exp(-h*omega_n*z)*sin(h*omega_d &
      )/omega_d
A[1,1] = (omega_d*cos(h*omega_d) - omega_n*z*sin(h*omega_d))*exp(-h* &
      omega_n*z)/omega_d
A[1,2] = (omega_d**2 + omega_n**2*z**2)*exp(-h*omega_n*z)*sin(h*omega_d) &
      /(k*omega_d)
A[1,3] = (-c*omega_d**2*sin(h*omega_d) - c*omega_n**2*z**2*sin(h*omega_d &
      ) - k*omega_d*cos(h*omega_d) + k*omega_n*z*sin(h*omega_d))*exp(-h &
      *omega_n*z)/(h*k**2*omega_d)

In [49]:

    

xh.expand().diff(dp).simplify()


    

    
Out[49]:





$$\frac{e^{- h \omega_{n} z}}{h k^{2} \omega_{d}} \left(c \omega_{d} \cos{\left (h \omega_{d} \right )} + c \omega_{n} z \sin{\left (h \omega_{d} \right )} - k \sin{\left (h \omega_{d} \right )}\right)$$

In [ ]:

    

.sim