#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Mathmatical models governed by ordinary differential equations.
Describes initial value problems as systems of first order ordinary
differential equations.
:Date: 2020-01-08
.. module:: ode_model
:platform: *nix, Windows
:synopsis: Models of ordinary differential equations.
.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
"""
from numpy import array, cos, sin, dot, square
from numpy.linalg import inv
[docs]def disk(d, e, T, method=''):
r"""Rotation of an eccentric disk.
:param d: diameter
:type d: float
:param e: eccentricity
:type e: float
:param T: torque
:type T: float
:param method: the method to use, default = ''.
:type method:
:returns: disk function. This function is independent of the time.
* For method = '': f(x, t=0) -> (xp1, xp2, xp3, xp4, xp5, xp6).
x is (x, y, phi, x', y', phi') and the return values are (x',
y', phi', x'', y'', phi'')
* For method = 'nm': f(xn, xpn, xppn, t=0) -> (N, dN, dNp, dNpp).
xn are the values of the function (x, y, phi), xpn are first
derivative values of the function (x', y', phi') and xppn are
the second derivative values of the function (x'', y'', phi'').
The return values are (N, dN, dNp, dNpp)
* For method = 'nmmdk': f(xn, xpn, xppn, t=0) ->
(rm, rmx, rmxpp, rd, rdx, rdxp, rk, rkx, f).
xn are the values of the function (x, y, phi), xpn are first
derivative values of the function (x', y', phi') and xppn are
the second derivative values of the function (x'', y'', phi'').
The return values are (rm, rmx, rmxpp, rd, rdx, rdxp, rk, rkx,
f)
:rtype: function
Model
.. math ::
\begin{vmatrix}
\ddot{x} + \cos(\varphi)\ddot{\varphi} + 2d \,\dot{x} - \sin(\varphi) \,\dot{\varphi}^2 + 2d\cos(\varphi)\, \dot{\varphi} + x &=&
0 \\
\ddot{y} - \sin(\varphi)\ddot{\varphi} + 2d \,\dot{y} - \cos(\varphi) \,\dot{\varphi}^2 + 2d\sin(\varphi)\, \dot{\varphi} + y &=&
0 \\
\ddot{\varphi} + e\,y\sin(\varphi) - e\,x\cos(\varphi) &=& t
\end{vmatrix}
\\
\begin{vmatrix}
\ddot{x} + \cos(\varphi)\ddot{\varphi} &=&
-2d \,\dot{x} + \sin(\varphi) \,\dot{\varphi}^2 -2d\cos(\varphi)\, \dot{\varphi} - x \\
\ddot{y} - \sin(\varphi)\ddot{\varphi} &=&
-2d \,\dot{y} + \cos(\varphi) \,\dot{\varphi}^2 -2d\sin(\varphi)\, \dot{\varphi} - y \\
\ddot{\varphi} &=& t - e\,y\sin(\varphi) + e\,x\cos(\varphi)
\end{vmatrix}
.. math ::
\mathbf{M}(\mathbf{x}) \cdot \mathbf{\ddot{x}} &=
\mathbf{f}(\mathbf{x}, \mathbf{\dot{x}})
\\
\begin{bmatrix}
1 & 0 & \cos \varphi \\
0 & 1 & -\sin \varphi \\
0 & 0 & 1
\end{bmatrix} \cdot
\begin{bmatrix}
\ddot{x} \\ \ddot{y} \\ \ddot{\varphi}
\end{bmatrix} &= \begin{bmatrix}
-2d \,\dot{x} + \sin(\varphi) \,\dot{\varphi}^2 -2d\cos(\varphi)\, \dot{\varphi} - x \\
-2d \,\dot{y} + \cos(\varphi) \,\dot{\varphi}^2 -2d\sin(\varphi)\, \dot{\varphi} - y \\
t - e\,y\sin(\varphi) + e\,x\cos(\varphi)
\end{bmatrix}
returns
.. math ::
x_1 &= x &\quad x_4 &= \dot{x}_1 = \dot{x} &\quad \dot{x}_4 &= \ddot{x} \\
x_2 &= y &\quad x_5 &= \dot{x}_2 = \dot{y} &\quad \dot{x}_5 &= \ddot{y} \\
x_3 &= \varphi &\quad x_6 &= \dot{x}_3 = \dot{\varphi} &\quad \dot{x}_6 &= \ddot{\varphi} \\
.. math ::
\dot{q} &= f(x) \\
\begin{bmatrix}
\dot{x}_1 \\
\dot{x}_2 \\
\dot{x}_3 \\
\dot{x}_4 \\
\dot{x}_5 \\
\dot{x}_6
\end{bmatrix} &= \begin{bmatrix}
x_4 \\
x_5 \\
x_6 \\
\begin{bmatrix}
1 & 0 & \cos x_3 \\
0 & 1 & -\sin x_3 \\
0 & 0 & 1
\end{bmatrix}^{-1} \cdot \begin{bmatrix}
-2d \,x_4 + \sin(x_3) \,x_6^2 -2d\cos(x_3)\, x_6 - x_1 \\
-2d \,x_5 + \cos(x_3) \,x_6^2 -2d\sin(x_3)\, x_6 - x_2 \\
t - e\,x_2\sin(x_3) + e\,x_1\cos(x_3)
\end{bmatrix}
\end{bmatrix}
Three explicit differential equations of order 2 reducted to a
system of 3x2 first-order differential equations.
.. seealso::
:meth:`pylib.numerical.ode.newmark_newtonraphson` and
:meth:`pylib.numerical.ode.newmark_newtonraphson_mdk`
"""
p = (d, e, T)
if method == '':
def f(x, t=0):
xp1 = x[3]
xp2 = x[4]
xp3 = x[5]
M = array([[1, 0, cos(x[2])], [0, 1, -sin(x[2])], [0, 0, 1]])
y = array([[-2*p[0]*x[3]+sin(x[2])*x[5]**2-2*p[0]*cos(x[2])*x[5]-x[0]], \
[-2*p[0]*x[4]+cos(x[2])*x[5]**2-2*p[0]*sin(x[2])*x[5]-x[1]], \
[p[2]-p[1]*x[1]*sin(x[2])+p[1]*x[0]*cos(x[2])]])
xp46 = dot(inv(M), y)
xp4, xp5, xp6 = xp46.reshape(-1,).tolist() # 2d array to 1d array to list
return xp1, xp2, xp3, xp4, xp5, xp6
elif method == 'nm':
def f(xn, xpn, xppn, t=0):
N = array([[xppn[0]+cos(xn[2])*xppn[2]+2*p[0]*xpn[0]+2*p[0]*cos(xn[2])*xpn[2]-sin(xn[2])*square(xpn[2])+xn[0]],
[xppn[1]-sin(xn[2])*xppn[2]+2*p[0]*xpn[1]-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])+xn[1]],
[xppn[2]+p[1]*(-cos(xn[2])*xn[0]+sin(xn[2])*xn[1])-p[2]]])
dN = array([[1, 0, -sin(xn[2]*xppn[2])-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
[0, 1, -cos(xn[2]*xppn[2])-2*p[0]*cos(xn[2])*xpn[2]+sin(xn[2])*square(xpn[2])],
[-p[1]*cos(xn[2]), p[1]*cos(xn[2]), p[1]*(sin(xn[2])*xn[0]+cos(xn[2])*xn[1])]])
dNp = array([[2*p[0], 0, 2*p[0]*cos(xn[2])-2*sin(xn[2])*xpn[2]],
[0, 2*p[0], -2*p[0]*sin(xn[2])-2*cos(xn[2])*xpn[2]],
[0, 0, 0]])
dNpp = array([[1, 0, cos(xn[2])],
[0, 1, -sin(xn[2])],
[0, 0, 1]])
return N, dN, dNp, dNpp
elif method == 'nmmdk':
def f(xn, xpn, xppn, t=0):
rm = array([[xppn[0]+cos(xn[2])*xppn[2]],
[xppn[1]-sin(xn[2])*xppn[2]],
[xppn[2]]])
rmx = array([[0, 0, -sin(xn[2]*xppn[2])],
[0, 0, -cos(xn[2]*xppn[2])],
[0, 0, 0]])
rmxpp = array([[1, 0, cos(xn[2])],
[0, 1, -sin(xn[2])],
[0, 0, 1]])
rd = array([[2*p[0]*xpn[0]+2*p[0]*cos(xn[2])*xpn[2]-sin(xn[2])*square(xpn[2])],
[2*p[0]*xpn[1]-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
[0]])
rdx = array([[0, 0, -2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
[0, 0, -2*p[0]*cos(xn[2])*xpn[2]+sin(xn[2])*square(xpn[2])],
[0, 0, 0]])
rdxp = array([[2*p[0], 0, 2*p[0]*cos(xn[2])-2*sin(xn[2])*xpn[2]],
[0, 2*p[0], -2*p[0]*sin(xn[2])-2*cos(xn[2])*xpn[2]],
[0, 0, 0]])
rk = array([[xn[0]],
[xn[1]],
[p[1]*(-cos(xn[2])*xn[0]+sin(xn[2])*xn[1])]])
rkx = array([[ 1, 0, 0],
[ 0, 1, 0],
[-p[1]*cos(xn[2]), p[1]*cos(xn[2]), p[1]*(sin(xn[2])*xn[0]+cos(xn[2])*xn[1])]])
f = array([[0], [0], [p[2]]])
return rm, rmx, rmxpp, rd, rdx, rdxp, rk, rkx, f
return f