AI Newsletter Digest improvements: fixed QP soft line break decoding, URL extraction, and content cleaning

This commit is contained in:
Krilly
2026-03-04 13:29:22 +00:00
parent 29a98137a7
commit 57dd294675
13706 changed files with 2114953 additions and 237629 deletions
@@ -0,0 +1,12 @@
"""
A module that helps solving problems in physics.
"""
from . import units
from .matrices import mgamma, msigma, minkowski_tensor, mdft
__all__ = [
'units',
'mgamma', 'msigma', 'minkowski_tensor', 'mdft',
]
@@ -0,0 +1,53 @@
"""Biomechanics extension for SymPy.
Includes biomechanics-related constructs which allows users to extend multibody
models created using `sympy.physics.mechanics` into biomechanical or
musculoskeletal models involding musculotendons and activation dynamics.
"""
from .activation import (
ActivationBase,
FirstOrderActivationDeGroote2016,
ZerothOrderActivation,
)
from .curve import (
CharacteristicCurveCollection,
CharacteristicCurveFunction,
FiberForceLengthActiveDeGroote2016,
FiberForceLengthPassiveDeGroote2016,
FiberForceLengthPassiveInverseDeGroote2016,
FiberForceVelocityDeGroote2016,
FiberForceVelocityInverseDeGroote2016,
TendonForceLengthDeGroote2016,
TendonForceLengthInverseDeGroote2016,
)
from .musculotendon import (
MusculotendonBase,
MusculotendonDeGroote2016,
MusculotendonFormulation,
)
__all__ = [
# Musculotendon characteristic curve functions
'CharacteristicCurveCollection',
'CharacteristicCurveFunction',
'FiberForceLengthActiveDeGroote2016',
'FiberForceLengthPassiveDeGroote2016',
'FiberForceLengthPassiveInverseDeGroote2016',
'FiberForceVelocityDeGroote2016',
'FiberForceVelocityInverseDeGroote2016',
'TendonForceLengthDeGroote2016',
'TendonForceLengthInverseDeGroote2016',
# Activation dynamics classes
'ActivationBase',
'FirstOrderActivationDeGroote2016',
'ZerothOrderActivation',
# Musculotendon classes
'MusculotendonBase',
'MusculotendonDeGroote2016',
'MusculotendonFormulation',
]
@@ -0,0 +1,53 @@
"""Mixin classes for sharing functionality between unrelated classes.
This module is named with a leading underscore to signify to users that it's
"private" and only intended for internal use by the biomechanics module.
"""
__all__ = ['_NamedMixin']
class _NamedMixin:
"""Mixin class for adding `name` properties.
Valid names, as will typically be used by subclasses as a suffix when
naming automatically-instantiated symbol attributes, must be nonzero length
strings.
Attributes
==========
name : str
The name identifier associated with the instance. Must be a string of
length at least 1.
"""
@property
def name(self) -> str:
"""The name associated with the class instance."""
return self._name
@name.setter
def name(self, name: str) -> None:
if hasattr(self, '_name'):
msg = (
f'Can\'t set attribute `name` to {repr(name)} as it is '
f'immutable.'
)
raise AttributeError(msg)
if not isinstance(name, str):
msg = (
f'Name {repr(name)} passed to `name` was of type '
f'{type(name)}, must be {str}.'
)
raise TypeError(msg)
if name in {''}:
msg = (
f'Name {repr(name)} is invalid, must be a nonzero length '
f'{type(str)}.'
)
raise ValueError(msg)
self._name = name
@@ -0,0 +1,869 @@
r"""Activation dynamics for musclotendon models.
Musculotendon models are able to produce active force when they are activated,
which is when a chemical process has taken place within the muscle fibers
causing them to voluntarily contract. Biologically this chemical process (the
diffusion of :math:`\textrm{Ca}^{2+}` ions) is not the input in the system,
electrical signals from the nervous system are. These are termed excitations.
Activation dynamics, which relates the normalized excitation level to the
normalized activation level, can be modeled by the models present in this
module.
"""
from abc import ABC, abstractmethod
from functools import cached_property
from sympy.core.symbol import Symbol
from sympy.core.numbers import Float, Integer, Rational
from sympy.functions.elementary.hyperbolic import tanh
from sympy.matrices.dense import MutableDenseMatrix as Matrix, zeros
from sympy.physics.biomechanics._mixin import _NamedMixin
from sympy.physics.mechanics import dynamicsymbols
__all__ = [
'ActivationBase',
'FirstOrderActivationDeGroote2016',
'ZerothOrderActivation',
]
class ActivationBase(ABC, _NamedMixin):
"""Abstract base class for all activation dynamics classes to inherit from.
Notes
=====
Instances of this class cannot be directly instantiated by users. However,
it can be used to created custom activation dynamics types through
subclassing.
"""
def __init__(self, name):
"""Initializer for ``ActivationBase``."""
self.name = str(name)
# Symbols
self._e = dynamicsymbols(f"e_{name}")
self._a = dynamicsymbols(f"a_{name}")
@classmethod
@abstractmethod
def with_defaults(cls, name):
"""Alternate constructor that provides recommended defaults for
constants."""
pass
@property
def excitation(self):
"""Dynamic symbol representing excitation.
Explanation
===========
The alias ``e`` can also be used to access the same attribute.
"""
return self._e
@property
def e(self):
"""Dynamic symbol representing excitation.
Explanation
===========
The alias ``excitation`` can also be used to access the same attribute.
"""
return self._e
@property
def activation(self):
"""Dynamic symbol representing activation.
Explanation
===========
The alias ``a`` can also be used to access the same attribute.
"""
return self._a
@property
def a(self):
"""Dynamic symbol representing activation.
Explanation
===========
The alias ``activation`` can also be used to access the same attribute.
"""
return self._a
@property
@abstractmethod
def order(self):
"""Order of the (differential) equation governing activation."""
pass
@property
@abstractmethod
def state_vars(self):
"""Ordered column matrix of functions of time that represent the state
variables.
Explanation
===========
The alias ``x`` can also be used to access the same attribute.
"""
pass
@property
@abstractmethod
def x(self):
"""Ordered column matrix of functions of time that represent the state
variables.
Explanation
===========
The alias ``state_vars`` can also be used to access the same attribute.
"""
pass
@property
@abstractmethod
def input_vars(self):
"""Ordered column matrix of functions of time that represent the input
variables.
Explanation
===========
The alias ``r`` can also be used to access the same attribute.
"""
pass
@property
@abstractmethod
def r(self):
"""Ordered column matrix of functions of time that represent the input
variables.
Explanation
===========
The alias ``input_vars`` can also be used to access the same attribute.
"""
pass
@property
@abstractmethod
def constants(self):
"""Ordered column matrix of non-time varying symbols present in ``M``
and ``F``.
Only symbolic constants are returned. If a numeric type (e.g. ``Float``)
has been used instead of ``Symbol`` for a constant then that attribute
will not be included in the matrix returned by this property. This is
because the primary use of this property attribute is to provide an
ordered sequence of the still-free symbols that require numeric values
during code generation.
Explanation
===========
The alias ``p`` can also be used to access the same attribute.
"""
pass
@property
@abstractmethod
def p(self):
"""Ordered column matrix of non-time varying symbols present in ``M``
and ``F``.
Only symbolic constants are returned. If a numeric type (e.g. ``Float``)
has been used instead of ``Symbol`` for a constant then that attribute
will not be included in the matrix returned by this property. This is
because the primary use of this property attribute is to provide an
ordered sequence of the still-free symbols that require numeric values
during code generation.
Explanation
===========
The alias ``constants`` can also be used to access the same attribute.
"""
pass
@property
@abstractmethod
def M(self):
"""Ordered square matrix of coefficients on the LHS of ``M x' = F``.
Explanation
===========
The square matrix that forms part of the LHS of the linear system of
ordinary differential equations governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
"""
pass
@property
@abstractmethod
def F(self):
"""Ordered column matrix of equations on the RHS of ``M x' = F``.
Explanation
===========
The column matrix that forms the RHS of the linear system of ordinary
differential equations governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
"""
pass
@abstractmethod
def rhs(self):
"""
Explanation
===========
The solution to the linear system of ordinary differential equations
governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
"""
pass
def __eq__(self, other):
"""Equality check for activation dynamics."""
if type(self) != type(other):
return False
if self.name != other.name:
return False
return True
def __repr__(self):
"""Default representation of activation dynamics."""
return f'{self.__class__.__name__}({self.name!r})'
class ZerothOrderActivation(ActivationBase):
"""Simple zeroth-order activation dynamics mapping excitation to
activation.
Explanation
===========
Zeroth-order activation dynamics are useful in instances where you want to
reduce the complexity of your musculotendon dynamics as they simple map
exictation to activation. As a result, no additional state equations are
introduced to your system. They also remove a potential source of delay
between the input and dynamics of your system as no (ordinary) differential
equations are involved.
"""
def __init__(self, name):
"""Initializer for ``ZerothOrderActivation``.
Parameters
==========
name : str
The name identifier associated with the instance. Must be a string
of length at least 1.
"""
super().__init__(name)
# Zeroth-order activation dynamics has activation equal excitation so
# overwrite the symbol for activation with the excitation symbol.
self._a = self._e
@classmethod
def with_defaults(cls, name):
"""Alternate constructor that provides recommended defaults for
constants.
Explanation
===========
As this concrete class doesn't implement any constants associated with
its dynamics, this ``classmethod`` simply creates a standard instance
of ``ZerothOrderActivation``. An implementation is provided to ensure
a consistent interface between all ``ActivationBase`` concrete classes.
"""
return cls(name)
@property
def order(self):
"""Order of the (differential) equation governing activation."""
return 0
@property
def state_vars(self):
"""Ordered column matrix of functions of time that represent the state
variables.
Explanation
===========
As zeroth-order activation dynamics simply maps excitation to
activation, this class has no associated state variables and so this
property return an empty column ``Matrix`` with shape (0, 1).
The alias ``x`` can also be used to access the same attribute.
"""
return zeros(0, 1)
@property
def x(self):
"""Ordered column matrix of functions of time that represent the state
variables.
Explanation
===========
As zeroth-order activation dynamics simply maps excitation to
activation, this class has no associated state variables and so this
property return an empty column ``Matrix`` with shape (0, 1).
The alias ``state_vars`` can also be used to access the same attribute.
"""
return zeros(0, 1)
@property
def input_vars(self):
"""Ordered column matrix of functions of time that represent the input
variables.
Explanation
===========
Excitation is the only input in zeroth-order activation dynamics and so
this property returns a column ``Matrix`` with one entry, ``e``, and
shape (1, 1).
The alias ``r`` can also be used to access the same attribute.
"""
return Matrix([self._e])
@property
def r(self):
"""Ordered column matrix of functions of time that represent the input
variables.
Explanation
===========
Excitation is the only input in zeroth-order activation dynamics and so
this property returns a column ``Matrix`` with one entry, ``e``, and
shape (1, 1).
The alias ``input_vars`` can also be used to access the same attribute.
"""
return Matrix([self._e])
@property
def constants(self):
"""Ordered column matrix of non-time varying symbols present in ``M``
and ``F``.
Only symbolic constants are returned. If a numeric type (e.g. ``Float``)
has been used instead of ``Symbol`` for a constant then that attribute
will not be included in the matrix returned by this property. This is
because the primary use of this property attribute is to provide an
ordered sequence of the still-free symbols that require numeric values
during code generation.
Explanation
===========
As zeroth-order activation dynamics simply maps excitation to
activation, this class has no associated constants and so this property
return an empty column ``Matrix`` with shape (0, 1).
The alias ``p`` can also be used to access the same attribute.
"""
return zeros(0, 1)
@property
def p(self):
"""Ordered column matrix of non-time varying symbols present in ``M``
and ``F``.
Only symbolic constants are returned. If a numeric type (e.g. ``Float``)
has been used instead of ``Symbol`` for a constant then that attribute
will not be included in the matrix returned by this property. This is
because the primary use of this property attribute is to provide an
ordered sequence of the still-free symbols that require numeric values
during code generation.
Explanation
===========
As zeroth-order activation dynamics simply maps excitation to
activation, this class has no associated constants and so this property
return an empty column ``Matrix`` with shape (0, 1).
The alias ``constants`` can also be used to access the same attribute.
"""
return zeros(0, 1)
@property
def M(self):
"""Ordered square matrix of coefficients on the LHS of ``M x' = F``.
Explanation
===========
The square matrix that forms part of the LHS of the linear system of
ordinary differential equations governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
As zeroth-order activation dynamics have no state variables, this
linear system has dimension 0 and therefore ``M`` is an empty square
``Matrix`` with shape (0, 0).
"""
return Matrix([])
@property
def F(self):
"""Ordered column matrix of equations on the RHS of ``M x' = F``.
Explanation
===========
The column matrix that forms the RHS of the linear system of ordinary
differential equations governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
As zeroth-order activation dynamics have no state variables, this
linear system has dimension 0 and therefore ``F`` is an empty column
``Matrix`` with shape (0, 1).
"""
return zeros(0, 1)
def rhs(self):
"""Ordered column matrix of equations for the solution of ``M x' = F``.
Explanation
===========
The solution to the linear system of ordinary differential equations
governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
As zeroth-order activation dynamics have no state variables, this
linear has dimension 0 and therefore this method returns an empty
column ``Matrix`` with shape (0, 1).
"""
return zeros(0, 1)
class FirstOrderActivationDeGroote2016(ActivationBase):
r"""First-order activation dynamics based on De Groote et al., 2016 [1]_.
Explanation
===========
Gives the first-order activation dynamics equation for the rate of change
of activation with respect to time as a function of excitation and
activation.
The function is defined by the equation:
.. math::
\frac{da}{dt} = \left(\frac{\frac{1}{2} + a0}{\tau_a \left(\frac{1}{2}
+ \frac{3a}{2}\right)} + \frac{\left(\frac{1}{2}
+ \frac{3a}{2}\right) \left(\frac{1}{2} - a0\right)}{\tau_d}\right)
\left(e - a\right)
where
.. math::
a0 = \frac{\tanh{\left(b \left(e - a\right) \right)}}{2}
with constant values of :math:`tau_a = 0.015`, :math:`tau_d = 0.060`, and
:math:`b = 10`.
References
==========
.. [1] De Groote, F., Kinney, A. L., Rao, A. V., & Fregly, B. J., Evaluation
of direct collocation optimal control problem formulations for
solving the muscle redundancy problem, Annals of biomedical
engineering, 44(10), (2016) pp. 2922-2936
"""
def __init__(self,
name,
activation_time_constant=None,
deactivation_time_constant=None,
smoothing_rate=None,
):
"""Initializer for ``FirstOrderActivationDeGroote2016``.
Parameters
==========
activation time constant : Symbol | Number | None
The value of the activation time constant governing the delay
between excitation and activation when excitation exceeds
activation.
deactivation time constant : Symbol | Number | None
The value of the deactivation time constant governing the delay
between excitation and activation when activation exceeds
excitation.
smoothing_rate : Symbol | Number | None
The slope of the hyperbolic tangent function used to smooth between
the switching of the equations where excitation exceed activation
and where activation exceeds excitation. The recommended value to
use is ``10``, but values between ``0.1`` and ``100`` can be used.
"""
super().__init__(name)
# Symbols
self.activation_time_constant = activation_time_constant
self.deactivation_time_constant = deactivation_time_constant
self.smoothing_rate = smoothing_rate
@classmethod
def with_defaults(cls, name):
r"""Alternate constructor that will use the published constants.
Explanation
===========
Returns an instance of ``FirstOrderActivationDeGroote2016`` using the
three constant values specified in the original publication.
These have the values:
:math:`tau_a = 0.015`
:math:`tau_d = 0.060`
:math:`b = 10`
"""
tau_a = Float('0.015')
tau_d = Float('0.060')
b = Float('10.0')
return cls(name, tau_a, tau_d, b)
@property
def activation_time_constant(self):
"""Delay constant for activation.
Explanation
===========
The alias ```tau_a`` can also be used to access the same attribute.
"""
return self._tau_a
@activation_time_constant.setter
def activation_time_constant(self, tau_a):
if hasattr(self, '_tau_a'):
msg = (
f'Can\'t set attribute `activation_time_constant` to '
f'{repr(tau_a)} as it is immutable and already has value '
f'{self._tau_a}.'
)
raise AttributeError(msg)
self._tau_a = Symbol(f'tau_a_{self.name}') if tau_a is None else tau_a
@property
def tau_a(self):
"""Delay constant for activation.
Explanation
===========
The alias ``activation_time_constant`` can also be used to access the
same attribute.
"""
return self._tau_a
@property
def deactivation_time_constant(self):
"""Delay constant for deactivation.
Explanation
===========
The alias ``tau_d`` can also be used to access the same attribute.
"""
return self._tau_d
@deactivation_time_constant.setter
def deactivation_time_constant(self, tau_d):
if hasattr(self, '_tau_d'):
msg = (
f'Can\'t set attribute `deactivation_time_constant` to '
f'{repr(tau_d)} as it is immutable and already has value '
f'{self._tau_d}.'
)
raise AttributeError(msg)
self._tau_d = Symbol(f'tau_d_{self.name}') if tau_d is None else tau_d
@property
def tau_d(self):
"""Delay constant for deactivation.
Explanation
===========
The alias ``deactivation_time_constant`` can also be used to access the
same attribute.
"""
return self._tau_d
@property
def smoothing_rate(self):
"""Smoothing constant for the hyperbolic tangent term.
Explanation
===========
The alias ``b`` can also be used to access the same attribute.
"""
return self._b
@smoothing_rate.setter
def smoothing_rate(self, b):
if hasattr(self, '_b'):
msg = (
f'Can\'t set attribute `smoothing_rate` to {b!r} as it is '
f'immutable and already has value {self._b!r}.'
)
raise AttributeError(msg)
self._b = Symbol(f'b_{self.name}') if b is None else b
@property
def b(self):
"""Smoothing constant for the hyperbolic tangent term.
Explanation
===========
The alias ``smoothing_rate`` can also be used to access the same
attribute.
"""
return self._b
@property
def order(self):
"""Order of the (differential) equation governing activation."""
return 1
@property
def state_vars(self):
"""Ordered column matrix of functions of time that represent the state
variables.
Explanation
===========
The alias ``x`` can also be used to access the same attribute.
"""
return Matrix([self._a])
@property
def x(self):
"""Ordered column matrix of functions of time that represent the state
variables.
Explanation
===========
The alias ``state_vars`` can also be used to access the same attribute.
"""
return Matrix([self._a])
@property
def input_vars(self):
"""Ordered column matrix of functions of time that represent the input
variables.
Explanation
===========
The alias ``r`` can also be used to access the same attribute.
"""
return Matrix([self._e])
@property
def r(self):
"""Ordered column matrix of functions of time that represent the input
variables.
Explanation
===========
The alias ``input_vars`` can also be used to access the same attribute.
"""
return Matrix([self._e])
@property
def constants(self):
"""Ordered column matrix of non-time varying symbols present in ``M``
and ``F``.
Only symbolic constants are returned. If a numeric type (e.g. ``Float``)
has been used instead of ``Symbol`` for a constant then that attribute
will not be included in the matrix returned by this property. This is
because the primary use of this property attribute is to provide an
ordered sequence of the still-free symbols that require numeric values
during code generation.
Explanation
===========
The alias ``p`` can also be used to access the same attribute.
"""
constants = [self._tau_a, self._tau_d, self._b]
symbolic_constants = [c for c in constants if not c.is_number]
return Matrix(symbolic_constants) if symbolic_constants else zeros(0, 1)
@property
def p(self):
"""Ordered column matrix of non-time varying symbols present in ``M``
and ``F``.
Explanation
===========
Only symbolic constants are returned. If a numeric type (e.g. ``Float``)
has been used instead of ``Symbol`` for a constant then that attribute
will not be included in the matrix returned by this property. This is
because the primary use of this property attribute is to provide an
ordered sequence of the still-free symbols that require numeric values
during code generation.
The alias ``constants`` can also be used to access the same attribute.
"""
constants = [self._tau_a, self._tau_d, self._b]
symbolic_constants = [c for c in constants if not c.is_number]
return Matrix(symbolic_constants) if symbolic_constants else zeros(0, 1)
@property
def M(self):
"""Ordered square matrix of coefficients on the LHS of ``M x' = F``.
Explanation
===========
The square matrix that forms part of the LHS of the linear system of
ordinary differential equations governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
"""
return Matrix([Integer(1)])
@property
def F(self):
"""Ordered column matrix of equations on the RHS of ``M x' = F``.
Explanation
===========
The column matrix that forms the RHS of the linear system of ordinary
differential equations governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
"""
return Matrix([self._da_eqn])
def rhs(self):
"""Ordered column matrix of equations for the solution of ``M x' = F``.
Explanation
===========
The solution to the linear system of ordinary differential equations
governing the activation dynamics:
``M(x, r, t, p) x' = F(x, r, t, p)``.
"""
return Matrix([self._da_eqn])
@cached_property
def _da_eqn(self):
HALF = Rational(1, 2)
a0 = HALF * tanh(self._b * (self._e - self._a))
a1 = (HALF + Rational(3, 2) * self._a)
a2 = (HALF + a0) / (self._tau_a * a1)
a3 = a1 * (HALF - a0) / self._tau_d
activation_dynamics_equation = (a2 + a3) * (self._e - self._a)
return activation_dynamics_equation
def __eq__(self, other):
"""Equality check for ``FirstOrderActivationDeGroote2016``."""
if type(self) != type(other):
return False
self_attrs = (self.name, self.tau_a, self.tau_d, self.b)
other_attrs = (other.name, other.tau_a, other.tau_d, other.b)
if self_attrs == other_attrs:
return True
return False
def __repr__(self):
"""Representation of ``FirstOrderActivationDeGroote2016``."""
return (
f'{self.__class__.__name__}({self.name!r}, '
f'activation_time_constant={self.tau_a!r}, '
f'deactivation_time_constant={self.tau_d!r}, '
f'smoothing_rate={self.b!r})'
)
@@ -0,0 +1,348 @@
"""Tests for the ``sympy.physics.biomechanics.activation.py`` module."""
import pytest
from sympy import Symbol
from sympy.core.numbers import Float, Integer, Rational
from sympy.functions.elementary.hyperbolic import tanh
from sympy.matrices import Matrix
from sympy.matrices.dense import zeros
from sympy.physics.mechanics import dynamicsymbols
from sympy.physics.biomechanics import (
ActivationBase,
FirstOrderActivationDeGroote2016,
ZerothOrderActivation,
)
from sympy.physics.biomechanics._mixin import _NamedMixin
from sympy.simplify.simplify import simplify
class TestZerothOrderActivation:
@staticmethod
def test_class():
assert issubclass(ZerothOrderActivation, ActivationBase)
assert issubclass(ZerothOrderActivation, _NamedMixin)
assert ZerothOrderActivation.__name__ == 'ZerothOrderActivation'
@pytest.fixture(autouse=True)
def _zeroth_order_activation_fixture(self):
self.name = 'name'
self.e = dynamicsymbols('e_name')
self.instance = ZerothOrderActivation(self.name)
def test_instance(self):
instance = ZerothOrderActivation(self.name)
assert isinstance(instance, ZerothOrderActivation)
def test_with_defaults(self):
instance = ZerothOrderActivation.with_defaults(self.name)
assert isinstance(instance, ZerothOrderActivation)
assert instance == ZerothOrderActivation(self.name)
def test_name(self):
assert hasattr(self.instance, 'name')
assert self.instance.name == self.name
def test_order(self):
assert hasattr(self.instance, 'order')
assert self.instance.order == 0
def test_excitation_attribute(self):
assert hasattr(self.instance, 'e')
assert hasattr(self.instance, 'excitation')
e_expected = dynamicsymbols('e_name')
assert self.instance.e == e_expected
assert self.instance.excitation == e_expected
assert self.instance.e is self.instance.excitation
def test_activation_attribute(self):
assert hasattr(self.instance, 'a')
assert hasattr(self.instance, 'activation')
a_expected = dynamicsymbols('e_name')
assert self.instance.a == a_expected
assert self.instance.activation == a_expected
assert self.instance.a is self.instance.activation is self.instance.e
def test_state_vars_attribute(self):
assert hasattr(self.instance, 'x')
assert hasattr(self.instance, 'state_vars')
assert self.instance.x == self.instance.state_vars
x_expected = zeros(0, 1)
assert self.instance.x == x_expected
assert self.instance.state_vars == x_expected
assert isinstance(self.instance.x, Matrix)
assert isinstance(self.instance.state_vars, Matrix)
assert self.instance.x.shape == (0, 1)
assert self.instance.state_vars.shape == (0, 1)
def test_input_vars_attribute(self):
assert hasattr(self.instance, 'r')
assert hasattr(self.instance, 'input_vars')
assert self.instance.r == self.instance.input_vars
r_expected = Matrix([self.e])
assert self.instance.r == r_expected
assert self.instance.input_vars == r_expected
assert isinstance(self.instance.r, Matrix)
assert isinstance(self.instance.input_vars, Matrix)
assert self.instance.r.shape == (1, 1)
assert self.instance.input_vars.shape == (1, 1)
def test_constants_attribute(self):
assert hasattr(self.instance, 'p')
assert hasattr(self.instance, 'constants')
assert self.instance.p == self.instance.constants
p_expected = zeros(0, 1)
assert self.instance.p == p_expected
assert self.instance.constants == p_expected
assert isinstance(self.instance.p, Matrix)
assert isinstance(self.instance.constants, Matrix)
assert self.instance.p.shape == (0, 1)
assert self.instance.constants.shape == (0, 1)
def test_M_attribute(self):
assert hasattr(self.instance, 'M')
M_expected = Matrix([])
assert self.instance.M == M_expected
assert isinstance(self.instance.M, Matrix)
assert self.instance.M.shape == (0, 0)
def test_F(self):
assert hasattr(self.instance, 'F')
F_expected = zeros(0, 1)
assert self.instance.F == F_expected
assert isinstance(self.instance.F, Matrix)
assert self.instance.F.shape == (0, 1)
def test_rhs(self):
assert hasattr(self.instance, 'rhs')
rhs_expected = zeros(0, 1)
rhs = self.instance.rhs()
assert rhs == rhs_expected
assert isinstance(rhs, Matrix)
assert rhs.shape == (0, 1)
def test_repr(self):
expected = 'ZerothOrderActivation(\'name\')'
assert repr(self.instance) == expected
class TestFirstOrderActivationDeGroote2016:
@staticmethod
def test_class():
assert issubclass(FirstOrderActivationDeGroote2016, ActivationBase)
assert issubclass(FirstOrderActivationDeGroote2016, _NamedMixin)
assert FirstOrderActivationDeGroote2016.__name__ == 'FirstOrderActivationDeGroote2016'
@pytest.fixture(autouse=True)
def _first_order_activation_de_groote_2016_fixture(self):
self.name = 'name'
self.e = dynamicsymbols('e_name')
self.a = dynamicsymbols('a_name')
self.tau_a = Symbol('tau_a')
self.tau_d = Symbol('tau_d')
self.b = Symbol('b')
self.instance = FirstOrderActivationDeGroote2016(
self.name,
self.tau_a,
self.tau_d,
self.b,
)
def test_instance(self):
instance = FirstOrderActivationDeGroote2016(self.name)
assert isinstance(instance, FirstOrderActivationDeGroote2016)
def test_with_defaults(self):
instance = FirstOrderActivationDeGroote2016.with_defaults(self.name)
assert isinstance(instance, FirstOrderActivationDeGroote2016)
assert instance.tau_a == Float('0.015')
assert instance.activation_time_constant == Float('0.015')
assert instance.tau_d == Float('0.060')
assert instance.deactivation_time_constant == Float('0.060')
assert instance.b == Float('10.0')
assert instance.smoothing_rate == Float('10.0')
def test_name(self):
assert hasattr(self.instance, 'name')
assert self.instance.name == self.name
def test_order(self):
assert hasattr(self.instance, 'order')
assert self.instance.order == 1
def test_excitation(self):
assert hasattr(self.instance, 'e')
assert hasattr(self.instance, 'excitation')
e_expected = dynamicsymbols('e_name')
assert self.instance.e == e_expected
assert self.instance.excitation == e_expected
assert self.instance.e is self.instance.excitation
def test_excitation_is_immutable(self):
with pytest.raises(AttributeError):
self.instance.e = None
with pytest.raises(AttributeError):
self.instance.excitation = None
def test_activation(self):
assert hasattr(self.instance, 'a')
assert hasattr(self.instance, 'activation')
a_expected = dynamicsymbols('a_name')
assert self.instance.a == a_expected
assert self.instance.activation == a_expected
def test_activation_is_immutable(self):
with pytest.raises(AttributeError):
self.instance.a = None
with pytest.raises(AttributeError):
self.instance.activation = None
@pytest.mark.parametrize(
'tau_a, expected',
[
(None, Symbol('tau_a_name')),
(Symbol('tau_a'), Symbol('tau_a')),
(Float('0.015'), Float('0.015')),
]
)
def test_activation_time_constant(self, tau_a, expected):
instance = FirstOrderActivationDeGroote2016(
'name', activation_time_constant=tau_a,
)
assert instance.tau_a == expected
assert instance.activation_time_constant == expected
assert instance.tau_a is instance.activation_time_constant
def test_activation_time_constant_is_immutable(self):
with pytest.raises(AttributeError):
self.instance.tau_a = None
with pytest.raises(AttributeError):
self.instance.activation_time_constant = None
@pytest.mark.parametrize(
'tau_d, expected',
[
(None, Symbol('tau_d_name')),
(Symbol('tau_d'), Symbol('tau_d')),
(Float('0.060'), Float('0.060')),
]
)
def test_deactivation_time_constant(self, tau_d, expected):
instance = FirstOrderActivationDeGroote2016(
'name', deactivation_time_constant=tau_d,
)
assert instance.tau_d == expected
assert instance.deactivation_time_constant == expected
assert instance.tau_d is instance.deactivation_time_constant
def test_deactivation_time_constant_is_immutable(self):
with pytest.raises(AttributeError):
self.instance.tau_d = None
with pytest.raises(AttributeError):
self.instance.deactivation_time_constant = None
@pytest.mark.parametrize(
'b, expected',
[
(None, Symbol('b_name')),
(Symbol('b'), Symbol('b')),
(Integer('10'), Integer('10')),
]
)
def test_smoothing_rate(self, b, expected):
instance = FirstOrderActivationDeGroote2016(
'name', smoothing_rate=b,
)
assert instance.b == expected
assert instance.smoothing_rate == expected
assert instance.b is instance.smoothing_rate
def test_smoothing_rate_is_immutable(self):
with pytest.raises(AttributeError):
self.instance.b = None
with pytest.raises(AttributeError):
self.instance.smoothing_rate = None
def test_state_vars(self):
assert hasattr(self.instance, 'x')
assert hasattr(self.instance, 'state_vars')
assert self.instance.x == self.instance.state_vars
x_expected = Matrix([self.a])
assert self.instance.x == x_expected
assert self.instance.state_vars == x_expected
assert isinstance(self.instance.x, Matrix)
assert isinstance(self.instance.state_vars, Matrix)
assert self.instance.x.shape == (1, 1)
assert self.instance.state_vars.shape == (1, 1)
def test_input_vars(self):
assert hasattr(self.instance, 'r')
assert hasattr(self.instance, 'input_vars')
assert self.instance.r == self.instance.input_vars
r_expected = Matrix([self.e])
assert self.instance.r == r_expected
assert self.instance.input_vars == r_expected
assert isinstance(self.instance.r, Matrix)
assert isinstance(self.instance.input_vars, Matrix)
assert self.instance.r.shape == (1, 1)
assert self.instance.input_vars.shape == (1, 1)
def test_constants(self):
assert hasattr(self.instance, 'p')
assert hasattr(self.instance, 'constants')
assert self.instance.p == self.instance.constants
p_expected = Matrix([self.tau_a, self.tau_d, self.b])
assert self.instance.p == p_expected
assert self.instance.constants == p_expected
assert isinstance(self.instance.p, Matrix)
assert isinstance(self.instance.constants, Matrix)
assert self.instance.p.shape == (3, 1)
assert self.instance.constants.shape == (3, 1)
def test_M(self):
assert hasattr(self.instance, 'M')
M_expected = Matrix([1])
assert self.instance.M == M_expected
assert isinstance(self.instance.M, Matrix)
assert self.instance.M.shape == (1, 1)
def test_F(self):
assert hasattr(self.instance, 'F')
da_expr = (
((1/(self.tau_a*(Rational(1, 2) + Rational(3, 2)*self.a)))
*(Rational(1, 2) + Rational(1, 2)*tanh(self.b*(self.e - self.a)))
+ ((Rational(1, 2) + Rational(3, 2)*self.a)/self.tau_d)
*(Rational(1, 2) - Rational(1, 2)*tanh(self.b*(self.e - self.a))))
*(self.e - self.a)
)
F_expected = Matrix([da_expr])
assert self.instance.F == F_expected
assert isinstance(self.instance.F, Matrix)
assert self.instance.F.shape == (1, 1)
def test_rhs(self):
assert hasattr(self.instance, 'rhs')
da_expr = (
((1/(self.tau_a*(Rational(1, 2) + Rational(3, 2)*self.a)))
*(Rational(1, 2) + Rational(1, 2)*tanh(self.b*(self.e - self.a)))
+ ((Rational(1, 2) + Rational(3, 2)*self.a)/self.tau_d)
*(Rational(1, 2) - Rational(1, 2)*tanh(self.b*(self.e - self.a))))
*(self.e - self.a)
)
rhs_expected = Matrix([da_expr])
rhs = self.instance.rhs()
assert rhs == rhs_expected
assert isinstance(rhs, Matrix)
assert rhs.shape == (1, 1)
assert simplify(self.instance.M.solve(self.instance.F) - rhs) == zeros(1)
def test_repr(self):
expected = (
'FirstOrderActivationDeGroote2016(\'name\', '
'activation_time_constant=tau_a, '
'deactivation_time_constant=tau_d, '
'smoothing_rate=b)'
)
assert repr(self.instance) == expected
@@ -0,0 +1,48 @@
"""Tests for the ``sympy.physics.biomechanics._mixin.py`` module."""
import pytest
from sympy.physics.biomechanics._mixin import _NamedMixin
class TestNamedMixin:
@staticmethod
def test_subclass():
class Subclass(_NamedMixin):
def __init__(self, name):
self.name = name
instance = Subclass('name')
assert instance.name == 'name'
@pytest.fixture(autouse=True)
def _named_mixin_fixture(self):
class Subclass(_NamedMixin):
def __init__(self, name):
self.name = name
self.Subclass = Subclass
@pytest.mark.parametrize('name', ['a', 'name', 'long_name'])
def test_valid_name_argument(self, name):
instance = self.Subclass(name)
assert instance.name == name
@pytest.mark.parametrize('invalid_name', [0, 0.0, None, False])
def test_invalid_name_argument_not_str(self, invalid_name):
with pytest.raises(TypeError):
_ = self.Subclass(invalid_name)
def test_invalid_name_argument_zero_length_str(self):
with pytest.raises(ValueError):
_ = self.Subclass('')
def test_name_attribute_is_immutable(self):
instance = self.Subclass('name')
with pytest.raises(AttributeError):
instance.name = 'new_name'
@@ -0,0 +1,837 @@
"""Tests for the ``sympy.physics.biomechanics.musculotendon.py`` module."""
import abc
import pytest
from sympy.core.expr import UnevaluatedExpr
from sympy.core.numbers import Float, Integer, Rational
from sympy.core.symbol import Symbol
from sympy.functions.elementary.exponential import exp
from sympy.functions.elementary.hyperbolic import tanh
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.elementary.trigonometric import sin
from sympy.matrices.dense import MutableDenseMatrix as Matrix, eye, zeros
from sympy.physics.biomechanics.activation import (
FirstOrderActivationDeGroote2016
)
from sympy.physics.biomechanics.curve import (
CharacteristicCurveCollection,
FiberForceLengthActiveDeGroote2016,
FiberForceLengthPassiveDeGroote2016,
FiberForceLengthPassiveInverseDeGroote2016,
FiberForceVelocityDeGroote2016,
FiberForceVelocityInverseDeGroote2016,
TendonForceLengthDeGroote2016,
TendonForceLengthInverseDeGroote2016,
)
from sympy.physics.biomechanics.musculotendon import (
MusculotendonBase,
MusculotendonDeGroote2016,
MusculotendonFormulation,
)
from sympy.physics.biomechanics._mixin import _NamedMixin
from sympy.physics.mechanics.actuator import ForceActuator
from sympy.physics.mechanics.pathway import LinearPathway
from sympy.physics.vector.frame import ReferenceFrame
from sympy.physics.vector.functions import dynamicsymbols
from sympy.physics.vector.point import Point
from sympy.simplify.simplify import simplify
class TestMusculotendonFormulation:
@staticmethod
def test_rigid_tendon_member():
assert MusculotendonFormulation(0) == 0
assert MusculotendonFormulation.RIGID_TENDON == 0
@staticmethod
def test_fiber_length_explicit_member():
assert MusculotendonFormulation(1) == 1
assert MusculotendonFormulation.FIBER_LENGTH_EXPLICIT == 1
@staticmethod
def test_tendon_force_explicit_member():
assert MusculotendonFormulation(2) == 2
assert MusculotendonFormulation.TENDON_FORCE_EXPLICIT == 2
@staticmethod
def test_fiber_length_implicit_member():
assert MusculotendonFormulation(3) == 3
assert MusculotendonFormulation.FIBER_LENGTH_IMPLICIT == 3
@staticmethod
def test_tendon_force_implicit_member():
assert MusculotendonFormulation(4) == 4
assert MusculotendonFormulation.TENDON_FORCE_IMPLICIT == 4
class TestMusculotendonBase:
@staticmethod
def test_is_abstract_base_class():
assert issubclass(MusculotendonBase, abc.ABC)
@staticmethod
def test_class():
assert issubclass(MusculotendonBase, ForceActuator)
assert issubclass(MusculotendonBase, _NamedMixin)
assert MusculotendonBase.__name__ == 'MusculotendonBase'
@staticmethod
def test_cannot_instantiate_directly():
with pytest.raises(TypeError):
_ = MusculotendonBase()
@pytest.mark.parametrize('musculotendon_concrete', [MusculotendonDeGroote2016])
class TestMusculotendonRigidTendon:
@pytest.fixture(autouse=True)
def _musculotendon_rigid_tendon_fixture(self, musculotendon_concrete):
self.name = 'name'
self.N = ReferenceFrame('N')
self.q = dynamicsymbols('q')
self.origin = Point('pO')
self.insertion = Point('pI')
self.insertion.set_pos(self.origin, self.q*self.N.x)
self.pathway = LinearPathway(self.origin, self.insertion)
self.activation = FirstOrderActivationDeGroote2016(self.name)
self.e = self.activation.excitation
self.a = self.activation.activation
self.tau_a = self.activation.activation_time_constant
self.tau_d = self.activation.deactivation_time_constant
self.b = self.activation.smoothing_rate
self.formulation = MusculotendonFormulation.RIGID_TENDON
self.l_T_slack = Symbol('l_T_slack')
self.F_M_max = Symbol('F_M_max')
self.l_M_opt = Symbol('l_M_opt')
self.v_M_max = Symbol('v_M_max')
self.alpha_opt = Symbol('alpha_opt')
self.beta = Symbol('beta')
self.instance = musculotendon_concrete(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=self.formulation,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
)
self.da_expr = (
(1/(self.tau_a*(Rational(1, 2) + Rational(3, 2)*self.a)))
*(Rational(1, 2) + Rational(1, 2)*tanh(self.b*(self.e - self.a)))
+ ((Rational(1, 2) + Rational(3, 2)*self.a)/self.tau_d)
*(Rational(1, 2) - Rational(1, 2)*tanh(self.b*(self.e - self.a)))
)*(self.e - self.a)
def test_state_vars(self):
assert hasattr(self.instance, 'x')
assert hasattr(self.instance, 'state_vars')
assert self.instance.x == self.instance.state_vars
x_expected = Matrix([self.a])
assert self.instance.x == x_expected
assert self.instance.state_vars == x_expected
assert isinstance(self.instance.x, Matrix)
assert isinstance(self.instance.state_vars, Matrix)
assert self.instance.x.shape == (1, 1)
assert self.instance.state_vars.shape == (1, 1)
def test_input_vars(self):
assert hasattr(self.instance, 'r')
assert hasattr(self.instance, 'input_vars')
assert self.instance.r == self.instance.input_vars
r_expected = Matrix([self.e])
assert self.instance.r == r_expected
assert self.instance.input_vars == r_expected
assert isinstance(self.instance.r, Matrix)
assert isinstance(self.instance.input_vars, Matrix)
assert self.instance.r.shape == (1, 1)
assert self.instance.input_vars.shape == (1, 1)
def test_constants(self):
assert hasattr(self.instance, 'p')
assert hasattr(self.instance, 'constants')
assert self.instance.p == self.instance.constants
p_expected = Matrix(
[
self.l_T_slack,
self.F_M_max,
self.l_M_opt,
self.v_M_max,
self.alpha_opt,
self.beta,
self.tau_a,
self.tau_d,
self.b,
Symbol('c_0_fl_T_name'),
Symbol('c_1_fl_T_name'),
Symbol('c_2_fl_T_name'),
Symbol('c_3_fl_T_name'),
Symbol('c_0_fl_M_pas_name'),
Symbol('c_1_fl_M_pas_name'),
Symbol('c_0_fl_M_act_name'),
Symbol('c_1_fl_M_act_name'),
Symbol('c_2_fl_M_act_name'),
Symbol('c_3_fl_M_act_name'),
Symbol('c_4_fl_M_act_name'),
Symbol('c_5_fl_M_act_name'),
Symbol('c_6_fl_M_act_name'),
Symbol('c_7_fl_M_act_name'),
Symbol('c_8_fl_M_act_name'),
Symbol('c_9_fl_M_act_name'),
Symbol('c_10_fl_M_act_name'),
Symbol('c_11_fl_M_act_name'),
Symbol('c_0_fv_M_name'),
Symbol('c_1_fv_M_name'),
Symbol('c_2_fv_M_name'),
Symbol('c_3_fv_M_name'),
]
)
assert self.instance.p == p_expected
assert self.instance.constants == p_expected
assert isinstance(self.instance.p, Matrix)
assert isinstance(self.instance.constants, Matrix)
assert self.instance.p.shape == (31, 1)
assert self.instance.constants.shape == (31, 1)
def test_M(self):
assert hasattr(self.instance, 'M')
M_expected = Matrix([1])
assert self.instance.M == M_expected
assert isinstance(self.instance.M, Matrix)
assert self.instance.M.shape == (1, 1)
def test_F(self):
assert hasattr(self.instance, 'F')
F_expected = Matrix([self.da_expr])
assert self.instance.F == F_expected
assert isinstance(self.instance.F, Matrix)
assert self.instance.F.shape == (1, 1)
def test_rhs(self):
assert hasattr(self.instance, 'rhs')
rhs_expected = Matrix([self.da_expr])
rhs = self.instance.rhs()
assert isinstance(rhs, Matrix)
assert rhs.shape == (1, 1)
assert simplify(rhs - rhs_expected) == zeros(1)
@pytest.mark.parametrize(
'musculotendon_concrete, curve',
[
(
MusculotendonDeGroote2016,
CharacteristicCurveCollection(
tendon_force_length=TendonForceLengthDeGroote2016,
tendon_force_length_inverse=TendonForceLengthInverseDeGroote2016,
fiber_force_length_passive=FiberForceLengthPassiveDeGroote2016,
fiber_force_length_passive_inverse=FiberForceLengthPassiveInverseDeGroote2016,
fiber_force_length_active=FiberForceLengthActiveDeGroote2016,
fiber_force_velocity=FiberForceVelocityDeGroote2016,
fiber_force_velocity_inverse=FiberForceVelocityInverseDeGroote2016,
),
)
],
)
class TestFiberLengthExplicit:
@pytest.fixture(autouse=True)
def _musculotendon_fiber_length_explicit_fixture(
self,
musculotendon_concrete,
curve,
):
self.name = 'name'
self.N = ReferenceFrame('N')
self.q = dynamicsymbols('q')
self.origin = Point('pO')
self.insertion = Point('pI')
self.insertion.set_pos(self.origin, self.q*self.N.x)
self.pathway = LinearPathway(self.origin, self.insertion)
self.activation = FirstOrderActivationDeGroote2016(self.name)
self.e = self.activation.excitation
self.a = self.activation.activation
self.tau_a = self.activation.activation_time_constant
self.tau_d = self.activation.deactivation_time_constant
self.b = self.activation.smoothing_rate
self.formulation = MusculotendonFormulation.FIBER_LENGTH_EXPLICIT
self.l_T_slack = Symbol('l_T_slack')
self.F_M_max = Symbol('F_M_max')
self.l_M_opt = Symbol('l_M_opt')
self.v_M_max = Symbol('v_M_max')
self.alpha_opt = Symbol('alpha_opt')
self.beta = Symbol('beta')
self.instance = musculotendon_concrete(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=self.formulation,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
with_defaults=True,
)
self.l_M_tilde = dynamicsymbols('l_M_tilde_name')
l_MT = self.pathway.length
l_M = self.l_M_tilde*self.l_M_opt
l_T = l_MT - sqrt(l_M**2 - (self.l_M_opt*sin(self.alpha_opt))**2)
fl_T = curve.tendon_force_length.with_defaults(l_T/self.l_T_slack)
fl_M_pas = curve.fiber_force_length_passive.with_defaults(self.l_M_tilde)
fl_M_act = curve.fiber_force_length_active.with_defaults(self.l_M_tilde)
v_M_tilde = curve.fiber_force_velocity_inverse.with_defaults(
((((fl_T*self.F_M_max)/((l_MT - l_T)/l_M))/self.F_M_max) - fl_M_pas)
/(self.a*fl_M_act)
)
self.dl_M_tilde_expr = (self.v_M_max/self.l_M_opt)*v_M_tilde
self.da_expr = (
(1/(self.tau_a*(Rational(1, 2) + Rational(3, 2)*self.a)))
*(Rational(1, 2) + Rational(1, 2)*tanh(self.b*(self.e - self.a)))
+ ((Rational(1, 2) + Rational(3, 2)*self.a)/self.tau_d)
*(Rational(1, 2) - Rational(1, 2)*tanh(self.b*(self.e - self.a)))
)*(self.e - self.a)
def test_state_vars(self):
assert hasattr(self.instance, 'x')
assert hasattr(self.instance, 'state_vars')
assert self.instance.x == self.instance.state_vars
x_expected = Matrix([self.l_M_tilde, self.a])
assert self.instance.x == x_expected
assert self.instance.state_vars == x_expected
assert isinstance(self.instance.x, Matrix)
assert isinstance(self.instance.state_vars, Matrix)
assert self.instance.x.shape == (2, 1)
assert self.instance.state_vars.shape == (2, 1)
def test_input_vars(self):
assert hasattr(self.instance, 'r')
assert hasattr(self.instance, 'input_vars')
assert self.instance.r == self.instance.input_vars
r_expected = Matrix([self.e])
assert self.instance.r == r_expected
assert self.instance.input_vars == r_expected
assert isinstance(self.instance.r, Matrix)
assert isinstance(self.instance.input_vars, Matrix)
assert self.instance.r.shape == (1, 1)
assert self.instance.input_vars.shape == (1, 1)
def test_constants(self):
assert hasattr(self.instance, 'p')
assert hasattr(self.instance, 'constants')
assert self.instance.p == self.instance.constants
p_expected = Matrix(
[
self.l_T_slack,
self.F_M_max,
self.l_M_opt,
self.v_M_max,
self.alpha_opt,
self.beta,
self.tau_a,
self.tau_d,
self.b,
]
)
assert self.instance.p == p_expected
assert self.instance.constants == p_expected
assert isinstance(self.instance.p, Matrix)
assert isinstance(self.instance.constants, Matrix)
assert self.instance.p.shape == (9, 1)
assert self.instance.constants.shape == (9, 1)
def test_M(self):
assert hasattr(self.instance, 'M')
M_expected = eye(2)
assert self.instance.M == M_expected
assert isinstance(self.instance.M, Matrix)
assert self.instance.M.shape == (2, 2)
def test_F(self):
assert hasattr(self.instance, 'F')
F_expected = Matrix([self.dl_M_tilde_expr, self.da_expr])
assert self.instance.F == F_expected
assert isinstance(self.instance.F, Matrix)
assert self.instance.F.shape == (2, 1)
def test_rhs(self):
assert hasattr(self.instance, 'rhs')
rhs_expected = Matrix([self.dl_M_tilde_expr, self.da_expr])
rhs = self.instance.rhs()
assert isinstance(rhs, Matrix)
assert rhs.shape == (2, 1)
assert simplify(rhs - rhs_expected) == zeros(2, 1)
@pytest.mark.parametrize(
'musculotendon_concrete, curve',
[
(
MusculotendonDeGroote2016,
CharacteristicCurveCollection(
tendon_force_length=TendonForceLengthDeGroote2016,
tendon_force_length_inverse=TendonForceLengthInverseDeGroote2016,
fiber_force_length_passive=FiberForceLengthPassiveDeGroote2016,
fiber_force_length_passive_inverse=FiberForceLengthPassiveInverseDeGroote2016,
fiber_force_length_active=FiberForceLengthActiveDeGroote2016,
fiber_force_velocity=FiberForceVelocityDeGroote2016,
fiber_force_velocity_inverse=FiberForceVelocityInverseDeGroote2016,
),
)
],
)
class TestTendonForceExplicit:
@pytest.fixture(autouse=True)
def _musculotendon_tendon_force_explicit_fixture(
self,
musculotendon_concrete,
curve,
):
self.name = 'name'
self.N = ReferenceFrame('N')
self.q = dynamicsymbols('q')
self.origin = Point('pO')
self.insertion = Point('pI')
self.insertion.set_pos(self.origin, self.q*self.N.x)
self.pathway = LinearPathway(self.origin, self.insertion)
self.activation = FirstOrderActivationDeGroote2016(self.name)
self.e = self.activation.excitation
self.a = self.activation.activation
self.tau_a = self.activation.activation_time_constant
self.tau_d = self.activation.deactivation_time_constant
self.b = self.activation.smoothing_rate
self.formulation = MusculotendonFormulation.TENDON_FORCE_EXPLICIT
self.l_T_slack = Symbol('l_T_slack')
self.F_M_max = Symbol('F_M_max')
self.l_M_opt = Symbol('l_M_opt')
self.v_M_max = Symbol('v_M_max')
self.alpha_opt = Symbol('alpha_opt')
self.beta = Symbol('beta')
self.instance = musculotendon_concrete(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=self.formulation,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
with_defaults=True,
)
self.F_T_tilde = dynamicsymbols('F_T_tilde_name')
l_T_tilde = curve.tendon_force_length_inverse.with_defaults(self.F_T_tilde)
l_MT = self.pathway.length
v_MT = self.pathway.extension_velocity
l_T = l_T_tilde*self.l_T_slack
l_M = sqrt((l_MT - l_T)**2 + (self.l_M_opt*sin(self.alpha_opt))**2)
l_M_tilde = l_M/self.l_M_opt
cos_alpha = (l_MT - l_T)/l_M
F_T = self.F_T_tilde*self.F_M_max
F_M = F_T/cos_alpha
F_M_tilde = F_M/self.F_M_max
fl_M_pas = curve.fiber_force_length_passive.with_defaults(l_M_tilde)
fl_M_act = curve.fiber_force_length_active.with_defaults(l_M_tilde)
fv_M = (F_M_tilde - fl_M_pas)/(self.a*fl_M_act)
v_M_tilde = curve.fiber_force_velocity_inverse.with_defaults(fv_M)
v_M = v_M_tilde*self.v_M_max
v_T = v_MT - v_M/cos_alpha
v_T_tilde = v_T/self.l_T_slack
self.dF_T_tilde_expr = (
Float('0.2')*Float('33.93669377311689')*exp(
Float('33.93669377311689')*UnevaluatedExpr(l_T_tilde - Float('0.995'))
)*v_T_tilde
)
self.da_expr = (
(1/(self.tau_a*(Rational(1, 2) + Rational(3, 2)*self.a)))
*(Rational(1, 2) + Rational(1, 2)*tanh(self.b*(self.e - self.a)))
+ ((Rational(1, 2) + Rational(3, 2)*self.a)/self.tau_d)
*(Rational(1, 2) - Rational(1, 2)*tanh(self.b*(self.e - self.a)))
)*(self.e - self.a)
def test_state_vars(self):
assert hasattr(self.instance, 'x')
assert hasattr(self.instance, 'state_vars')
assert self.instance.x == self.instance.state_vars
x_expected = Matrix([self.F_T_tilde, self.a])
assert self.instance.x == x_expected
assert self.instance.state_vars == x_expected
assert isinstance(self.instance.x, Matrix)
assert isinstance(self.instance.state_vars, Matrix)
assert self.instance.x.shape == (2, 1)
assert self.instance.state_vars.shape == (2, 1)
def test_input_vars(self):
assert hasattr(self.instance, 'r')
assert hasattr(self.instance, 'input_vars')
assert self.instance.r == self.instance.input_vars
r_expected = Matrix([self.e])
assert self.instance.r == r_expected
assert self.instance.input_vars == r_expected
assert isinstance(self.instance.r, Matrix)
assert isinstance(self.instance.input_vars, Matrix)
assert self.instance.r.shape == (1, 1)
assert self.instance.input_vars.shape == (1, 1)
def test_constants(self):
assert hasattr(self.instance, 'p')
assert hasattr(self.instance, 'constants')
assert self.instance.p == self.instance.constants
p_expected = Matrix(
[
self.l_T_slack,
self.F_M_max,
self.l_M_opt,
self.v_M_max,
self.alpha_opt,
self.beta,
self.tau_a,
self.tau_d,
self.b,
]
)
assert self.instance.p == p_expected
assert self.instance.constants == p_expected
assert isinstance(self.instance.p, Matrix)
assert isinstance(self.instance.constants, Matrix)
assert self.instance.p.shape == (9, 1)
assert self.instance.constants.shape == (9, 1)
def test_M(self):
assert hasattr(self.instance, 'M')
M_expected = eye(2)
assert self.instance.M == M_expected
assert isinstance(self.instance.M, Matrix)
assert self.instance.M.shape == (2, 2)
def test_F(self):
assert hasattr(self.instance, 'F')
F_expected = Matrix([self.dF_T_tilde_expr, self.da_expr])
assert self.instance.F == F_expected
assert isinstance(self.instance.F, Matrix)
assert self.instance.F.shape == (2, 1)
def test_rhs(self):
assert hasattr(self.instance, 'rhs')
rhs_expected = Matrix([self.dF_T_tilde_expr, self.da_expr])
rhs = self.instance.rhs()
assert isinstance(rhs, Matrix)
assert rhs.shape == (2, 1)
assert simplify(rhs - rhs_expected) == zeros(2, 1)
class TestMusculotendonDeGroote2016:
@staticmethod
def test_class():
assert issubclass(MusculotendonDeGroote2016, ForceActuator)
assert issubclass(MusculotendonDeGroote2016, _NamedMixin)
assert MusculotendonDeGroote2016.__name__ == 'MusculotendonDeGroote2016'
@staticmethod
def test_instance():
origin = Point('pO')
insertion = Point('pI')
insertion.set_pos(origin, dynamicsymbols('q')*ReferenceFrame('N').x)
pathway = LinearPathway(origin, insertion)
activation = FirstOrderActivationDeGroote2016('name')
l_T_slack = Symbol('l_T_slack')
F_M_max = Symbol('F_M_max')
l_M_opt = Symbol('l_M_opt')
v_M_max = Symbol('v_M_max')
alpha_opt = Symbol('alpha_opt')
beta = Symbol('beta')
instance = MusculotendonDeGroote2016(
'name',
pathway,
activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=l_T_slack,
peak_isometric_force=F_M_max,
optimal_fiber_length=l_M_opt,
maximal_fiber_velocity=v_M_max,
optimal_pennation_angle=alpha_opt,
fiber_damping_coefficient=beta,
)
assert isinstance(instance, MusculotendonDeGroote2016)
@pytest.fixture(autouse=True)
def _musculotendon_fixture(self):
self.name = 'name'
self.N = ReferenceFrame('N')
self.q = dynamicsymbols('q')
self.origin = Point('pO')
self.insertion = Point('pI')
self.insertion.set_pos(self.origin, self.q*self.N.x)
self.pathway = LinearPathway(self.origin, self.insertion)
self.activation = FirstOrderActivationDeGroote2016(self.name)
self.l_T_slack = Symbol('l_T_slack')
self.F_M_max = Symbol('F_M_max')
self.l_M_opt = Symbol('l_M_opt')
self.v_M_max = Symbol('v_M_max')
self.alpha_opt = Symbol('alpha_opt')
self.beta = Symbol('beta')
def test_with_defaults(self):
origin = Point('pO')
insertion = Point('pI')
insertion.set_pos(origin, dynamicsymbols('q')*ReferenceFrame('N').x)
pathway = LinearPathway(origin, insertion)
activation = FirstOrderActivationDeGroote2016('name')
l_T_slack = Symbol('l_T_slack')
F_M_max = Symbol('F_M_max')
l_M_opt = Symbol('l_M_opt')
v_M_max = Float('10.0')
alpha_opt = Float('0.0')
beta = Float('0.1')
instance = MusculotendonDeGroote2016.with_defaults(
'name',
pathway,
activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=l_T_slack,
peak_isometric_force=F_M_max,
optimal_fiber_length=l_M_opt,
)
assert instance.tendon_slack_length == l_T_slack
assert instance.peak_isometric_force == F_M_max
assert instance.optimal_fiber_length == l_M_opt
assert instance.maximal_fiber_velocity == v_M_max
assert instance.optimal_pennation_angle == alpha_opt
assert instance.fiber_damping_coefficient == beta
@pytest.mark.parametrize(
'l_T_slack, expected',
[
(None, Symbol('l_T_slack_name')),
(Symbol('l_T_slack'), Symbol('l_T_slack')),
(Rational(1, 2), Rational(1, 2)),
(Float('0.5'), Float('0.5')),
],
)
def test_tendon_slack_length(self, l_T_slack, expected):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
)
assert instance.l_T_slack == expected
assert instance.tendon_slack_length == expected
@pytest.mark.parametrize(
'F_M_max, expected',
[
(None, Symbol('F_M_max_name')),
(Symbol('F_M_max'), Symbol('F_M_max')),
(Integer(1000), Integer(1000)),
(Float('1000.0'), Float('1000.0')),
],
)
def test_peak_isometric_force(self, F_M_max, expected):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
)
assert instance.F_M_max == expected
assert instance.peak_isometric_force == expected
@pytest.mark.parametrize(
'l_M_opt, expected',
[
(None, Symbol('l_M_opt_name')),
(Symbol('l_M_opt'), Symbol('l_M_opt')),
(Rational(1, 2), Rational(1, 2)),
(Float('0.5'), Float('0.5')),
],
)
def test_optimal_fiber_length(self, l_M_opt, expected):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
)
assert instance.l_M_opt == expected
assert instance.optimal_fiber_length == expected
@pytest.mark.parametrize(
'v_M_max, expected',
[
(None, Symbol('v_M_max_name')),
(Symbol('v_M_max'), Symbol('v_M_max')),
(Integer(10), Integer(10)),
(Float('10.0'), Float('10.0')),
],
)
def test_maximal_fiber_velocity(self, v_M_max, expected):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
)
assert instance.v_M_max == expected
assert instance.maximal_fiber_velocity == expected
@pytest.mark.parametrize(
'alpha_opt, expected',
[
(None, Symbol('alpha_opt_name')),
(Symbol('alpha_opt'), Symbol('alpha_opt')),
(Integer(0), Integer(0)),
(Float('0.1'), Float('0.1')),
],
)
def test_optimal_pennation_angle(self, alpha_opt, expected):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=alpha_opt,
fiber_damping_coefficient=self.beta,
)
assert instance.alpha_opt == expected
assert instance.optimal_pennation_angle == expected
@pytest.mark.parametrize(
'beta, expected',
[
(None, Symbol('beta_name')),
(Symbol('beta'), Symbol('beta')),
(Integer(0), Integer(0)),
(Rational(1, 10), Rational(1, 10)),
(Float('0.1'), Float('0.1')),
],
)
def test_fiber_damping_coefficient(self, beta, expected):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=beta,
)
assert instance.beta == expected
assert instance.fiber_damping_coefficient == expected
def test_excitation(self):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
)
assert hasattr(instance, 'e')
assert hasattr(instance, 'excitation')
e_expected = dynamicsymbols('e_name')
assert instance.e == e_expected
assert instance.excitation == e_expected
assert instance.e is instance.excitation
def test_excitation_is_immutable(self):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
)
with pytest.raises(AttributeError):
instance.e = None
with pytest.raises(AttributeError):
instance.excitation = None
def test_activation(self):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
)
assert hasattr(instance, 'a')
assert hasattr(instance, 'activation')
a_expected = dynamicsymbols('a_name')
assert instance.a == a_expected
assert instance.activation == a_expected
def test_activation_is_immutable(self):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
)
with pytest.raises(AttributeError):
instance.a = None
with pytest.raises(AttributeError):
instance.activation = None
def test_repr(self):
instance = MusculotendonDeGroote2016(
self.name,
self.pathway,
self.activation,
musculotendon_dynamics=MusculotendonFormulation.RIGID_TENDON,
tendon_slack_length=self.l_T_slack,
peak_isometric_force=self.F_M_max,
optimal_fiber_length=self.l_M_opt,
maximal_fiber_velocity=self.v_M_max,
optimal_pennation_angle=self.alpha_opt,
fiber_damping_coefficient=self.beta,
)
expected = (
'MusculotendonDeGroote2016(\'name\', '
'pathway=LinearPathway(pO, pI), '
'activation_dynamics=FirstOrderActivationDeGroote2016(\'name\', '
'activation_time_constant=tau_a_name, '
'deactivation_time_constant=tau_d_name, '
'smoothing_rate=b_name), '
'musculotendon_dynamics=0, '
'tendon_slack_length=l_T_slack, '
'peak_isometric_force=F_M_max, '
'optimal_fiber_length=l_M_opt, '
'maximal_fiber_velocity=v_M_max, '
'optimal_pennation_angle=alpha_opt, '
'fiber_damping_coefficient=beta)'
)
assert repr(instance) == expected
@@ -0,0 +1,10 @@
__all__ = ['Beam',
'Truss',
'Cable',
'Arch'
]
from .beam import Beam
from .truss import Truss
from .cable import Cable
from .arch import Arch
@@ -0,0 +1,815 @@
"""
This module can be used to solve problems related
to 2D Cables.
"""
from sympy.core.sympify import sympify
from sympy.core.symbol import Symbol,symbols
from sympy import sin, cos, pi, atan, diff, Piecewise, solve, rad
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.solvers.solveset import linsolve
from sympy.matrices import Matrix
from sympy.plotting import plot
class Cable:
"""
Cables are structures in engineering that support
the applied transverse loads through the tensile
resistance developed in its members.
Cables are widely used in suspension bridges, tension
leg offshore platforms, transmission lines, and find
use in several other engineering applications.
Examples
========
A cable is supported at (0, 10) and (10, 10). Two point loads
acting vertically downwards act on the cable, one with magnitude 3 kN
and acting 2 meters from the left support and 3 meters below it, while
the other with magnitude 2 kN is 6 meters from the left support and
6 meters below it.
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(('A', 0, 10), ('B', 10, 10))
>>> c.apply_load(-1, ('P', 2, 7, 3, 270))
>>> c.apply_load(-1, ('Q', 6, 4, 2, 270))
>>> c.loads
{'distributed': {}, 'point_load': {'P': [3, 270], 'Q': [2, 270]}}
>>> c.loads_position
{'P': [2, 7], 'Q': [6, 4]}
"""
def __init__(self, support_1, support_2):
"""
Initializes the class.
Parameters
==========
support_1 and support_2 are tuples of the form
(label, x, y), where
label : String or symbol
The label of the support
x : Sympifyable
The x coordinate of the position of the support
y : Sympifyable
The y coordinate of the position of the support
"""
self._left_support = []
self._right_support = []
self._supports = {}
self._support_labels = []
self._loads = {"distributed": {}, "point_load": {}}
self._loads_position = {}
self._length = 0
self._reaction_loads = {}
self._tension = {}
self._lowest_x_global = sympify(0)
self._lowest_y_global = sympify(0)
self._cable_eqn = None
self._tension_func = None
if support_1[0] == support_2[0]:
raise ValueError("Supports can not have the same label")
elif support_1[1] == support_2[1]:
raise ValueError("Supports can not be at the same location")
x1 = sympify(support_1[1])
y1 = sympify(support_1[2])
self._supports[support_1[0]] = [x1, y1]
x2 = sympify(support_2[1])
y2 = sympify(support_2[2])
self._supports[support_2[0]] = [x2, y2]
if support_1[1] < support_2[1]:
self._left_support.append(x1)
self._left_support.append(y1)
self._right_support.append(x2)
self._right_support.append(y2)
self._support_labels.append(support_1[0])
self._support_labels.append(support_2[0])
else:
self._left_support.append(x2)
self._left_support.append(y2)
self._right_support.append(x1)
self._right_support.append(y1)
self._support_labels.append(support_2[0])
self._support_labels.append(support_1[0])
for i in self._support_labels:
self._reaction_loads[Symbol("R_"+ i +"_x")] = 0
self._reaction_loads[Symbol("R_"+ i +"_y")] = 0
@property
def supports(self):
"""
Returns the supports of the cable along with their
positions.
"""
return self._supports
@property
def left_support(self):
"""
Returns the position of the left support.
"""
return self._left_support
@property
def right_support(self):
"""
Returns the position of the right support.
"""
return self._right_support
@property
def loads(self):
"""
Returns the magnitude and direction of the loads
acting on the cable.
"""
return self._loads
@property
def loads_position(self):
"""
Returns the position of the point loads acting on the
cable.
"""
return self._loads_position
@property
def length(self):
"""
Returns the length of the cable.
"""
return self._length
@property
def reaction_loads(self):
"""
Returns the reaction forces at the supports, which are
initialized to 0.
"""
return self._reaction_loads
@property
def tension(self):
"""
Returns the tension developed in the cable due to the loads
applied.
"""
return self._tension
def tension_at(self, x):
"""
Returns the tension at a given value of x developed due to
distributed load.
"""
if 'distributed' not in self._tension.keys():
raise ValueError("No distributed load added or solve method not called")
if x > self._right_support[0] or x < self._left_support[0]:
raise ValueError("The value of x should be between the two supports")
A = self._tension['distributed']
X = Symbol('X')
return A.subs({X:(x-self._lowest_x_global)})
def apply_length(self, length):
"""
This method specifies the length of the cable
Parameters
==========
length : Sympifyable
The length of the cable
Examples
========
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(('A', 0, 10), ('B', 10, 10))
>>> c.apply_length(20)
>>> c.length
20
"""
dist = ((self._left_support[0] - self._right_support[0])**2
- (self._left_support[1] - self._right_support[1])**2)**(1/2)
if length < dist:
raise ValueError("length should not be less than the distance between the supports")
self._length = length
def change_support(self, label, new_support):
"""
This method changes the mentioned support with a new support.
Parameters
==========
label: String or symbol
The label of the support to be changed
new_support: Tuple of the form (new_label, x, y)
new_label: String or symbol
The label of the new support
x: Sympifyable
The x-coordinate of the position of the new support.
y: Sympifyable
The y-coordinate of the position of the new support.
Examples
========
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(('A', 0, 10), ('B', 10, 10))
>>> c.supports
{'A': [0, 10], 'B': [10, 10]}
>>> c.change_support('B', ('C', 5, 6))
>>> c.supports
{'A': [0, 10], 'C': [5, 6]}
"""
if label not in self._supports:
raise ValueError("No support exists with the given label")
i = self._support_labels.index(label)
rem_label = self._support_labels[(i+1)%2]
x1 = self._supports[rem_label][0]
y1 = self._supports[rem_label][1]
x = sympify(new_support[1])
y = sympify(new_support[2])
for l in self._loads_position:
if l[0] >= max(x, x1) or l[0] <= min(x, x1):
raise ValueError("The change in support will throw an existing load out of range")
self._supports.pop(label)
self._left_support.clear()
self._right_support.clear()
self._reaction_loads.clear()
self._support_labels.remove(label)
self._supports[new_support[0]] = [x, y]
if x1 < x:
self._left_support.append(x1)
self._left_support.append(y1)
self._right_support.append(x)
self._right_support.append(y)
self._support_labels.append(new_support[0])
else:
self._left_support.append(x)
self._left_support.append(y)
self._right_support.append(x1)
self._right_support.append(y1)
self._support_labels.insert(0, new_support[0])
for i in self._support_labels:
self._reaction_loads[Symbol("R_"+ i +"_x")] = 0
self._reaction_loads[Symbol("R_"+ i +"_y")] = 0
def apply_load(self, order, load):
"""
This method adds load to the cable.
Parameters
==========
order : Integer
The order of the applied load.
- For point loads, order = -1
- For distributed load, order = 0
load : tuple
* For point loads, load is of the form (label, x, y, magnitude, direction), where:
label : String or symbol
The label of the load
x : Sympifyable
The x coordinate of the position of the load
y : Sympifyable
The y coordinate of the position of the load
magnitude : Sympifyable
The magnitude of the load. It must always be positive
direction : Sympifyable
The angle, in degrees, that the load vector makes with the horizontal
in the counter-clockwise direction. It takes the values 0 to 360,
inclusive.
* For uniformly distributed load, load is of the form (label, magnitude)
label : String or symbol
The label of the load
magnitude : Sympifyable
The magnitude of the load. It must always be positive
Examples
========
For a point load of magnitude 12 units inclined at 30 degrees with the horizontal:
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(('A', 0, 10), ('B', 10, 10))
>>> c.apply_load(-1, ('Z', 5, 5, 12, 30))
>>> c.loads
{'distributed': {}, 'point_load': {'Z': [12, 30]}}
>>> c.loads_position
{'Z': [5, 5]}
For a uniformly distributed load of magnitude 9 units:
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(('A', 0, 10), ('B', 10, 10))
>>> c.apply_load(0, ('X', 9))
>>> c.loads
{'distributed': {'X': 9}, 'point_load': {}}
"""
if order == -1:
if len(self._loads["distributed"]) != 0:
raise ValueError("Distributed load already exists")
label = load[0]
if label in self._loads["point_load"]:
raise ValueError("Label already exists")
x = sympify(load[1])
y = sympify(load[2])
if x > self._right_support[0] or x < self._left_support[0]:
raise ValueError("The load should be positioned between the supports")
magnitude = sympify(load[3])
direction = sympify(load[4])
self._loads["point_load"][label] = [magnitude, direction]
self._loads_position[label] = [x, y]
elif order == 0:
if len(self._loads_position) != 0:
raise ValueError("Point load(s) already exist")
label = load[0]
if label in self._loads["distributed"]:
raise ValueError("Label already exists")
magnitude = sympify(load[1])
self._loads["distributed"][label] = magnitude
else:
raise ValueError("Order should be either -1 or 0")
def remove_loads(self, *args):
"""
This methods removes the specified loads.
Parameters
==========
This input takes multiple label(s) as input
label(s): String or symbol
The label(s) of the loads to be removed.
Examples
========
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(('A', 0, 10), ('B', 10, 10))
>>> c.apply_load(-1, ('Z', 5, 5, 12, 30))
>>> c.loads
{'distributed': {}, 'point_load': {'Z': [12, 30]}}
>>> c.remove_loads('Z')
>>> c.loads
{'distributed': {}, 'point_load': {}}
"""
for i in args:
if len(self._loads_position) == 0:
if i not in self._loads['distributed']:
raise ValueError("Error removing load " + i + ": no such load exists")
else:
self._loads['disrtibuted'].pop(i)
else:
if i not in self._loads['point_load']:
raise ValueError("Error removing load " + i + ": no such load exists")
else:
self._loads['point_load'].pop(i)
self._loads_position.pop(i)
def solve(self, *args):
"""
This method solves for the reaction forces at the supports, the tension developed in
the cable, and updates the length of the cable.
Parameters
==========
This method requires no input when solving for point loads
For distributed load, the x and y coordinates of the lowest point of the cable are
required as
x: Sympifyable
The x coordinate of the lowest point
y: Sympifyable
The y coordinate of the lowest point
Examples
========
For point loads,
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(("A", 0, 10), ("B", 10, 10))
>>> c.apply_load(-1, ('Z', 2, 7.26, 3, 270))
>>> c.apply_load(-1, ('X', 4, 6, 8, 270))
>>> c.solve()
>>> c.tension
{A_Z: 8.91403453669861, X_B: 19*sqrt(13)/10, Z_X: 4.79150773600774}
>>> c.reaction_loads
{R_A_x: -5.25547445255474, R_A_y: 7.2, R_B_x: 5.25547445255474, R_B_y: 3.8}
>>> c.length
5.7560958484519 + 2*sqrt(13)
For distributed load,
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c=Cable(("A", 0, 40),("B", 100, 20))
>>> c.apply_load(0, ("X", 850))
>>> c.solve(58.58)
>>> c.tension
{'distributed': 36465.0*sqrt(0.00054335718671383*X**2 + 1)}
>>> c.tension_at(0)
61717.4130533677
>>> c.reaction_loads
{R_A_x: 36465.0, R_A_y: -49793.0, R_B_x: 44399.9537590861, R_B_y: 42868.2071025955}
"""
if len(self._loads_position) != 0:
sorted_position = sorted(self._loads_position.items(), key = lambda item : item[1][0])
sorted_position.append(self._support_labels[1])
sorted_position.insert(0, self._support_labels[0])
self._tension.clear()
moment_sum_from_left_support = 0
moment_sum_from_right_support = 0
F_x = 0
F_y = 0
self._length = 0
tension_func = []
x = symbols('x')
for i in range(1, len(sorted_position)-1):
if i == 1:
self._length+=sqrt((self._left_support[0] - self._loads_position[sorted_position[i][0]][0])**2 + (self._left_support[1] - self._loads_position[sorted_position[i][0]][1])**2)
else:
self._length+=sqrt((self._loads_position[sorted_position[i-1][0]][0] - self._loads_position[sorted_position[i][0]][0])**2 + (self._loads_position[sorted_position[i-1][0]][1] - self._loads_position[sorted_position[i][0]][1])**2)
if i == len(sorted_position)-2:
self._length+=sqrt((self._right_support[0] - self._loads_position[sorted_position[i][0]][0])**2 + (self._right_support[1] - self._loads_position[sorted_position[i][0]][1])**2)
moment_sum_from_left_support += self._loads['point_load'][sorted_position[i][0]][0] * cos(pi * self._loads['point_load'][sorted_position[i][0]][1] / 180) * abs(self._left_support[1] - self._loads_position[sorted_position[i][0]][1])
moment_sum_from_left_support += self._loads['point_load'][sorted_position[i][0]][0] * sin(pi * self._loads['point_load'][sorted_position[i][0]][1] / 180) * abs(self._left_support[0] - self._loads_position[sorted_position[i][0]][0])
F_x += self._loads['point_load'][sorted_position[i][0]][0] * cos(pi * self._loads['point_load'][sorted_position[i][0]][1] / 180)
F_y += self._loads['point_load'][sorted_position[i][0]][0] * sin(pi * self._loads['point_load'][sorted_position[i][0]][1] / 180)
label = Symbol(sorted_position[i][0]+"_"+sorted_position[i+1][0])
y2 = self._loads_position[sorted_position[i][0]][1]
x2 = self._loads_position[sorted_position[i][0]][0]
y1 = 0
x1 = 0
if i == len(sorted_position)-2:
x1 = self._right_support[0]
y1 = self._right_support[1]
else:
x1 = self._loads_position[sorted_position[i+1][0]][0]
y1 = self._loads_position[sorted_position[i+1][0]][1]
angle_with_horizontal = atan((y1 - y2)/(x1 - x2))
tension = -(moment_sum_from_left_support)/(abs(self._left_support[1] - self._loads_position[sorted_position[i][0]][1])*cos(angle_with_horizontal) + abs(self._left_support[0] - self._loads_position[sorted_position[i][0]][0])*sin(angle_with_horizontal))
self._tension[label] = tension
tension_func.append((tension, x<=x1))
moment_sum_from_right_support += self._loads['point_load'][sorted_position[i][0]][0] * cos(pi * self._loads['point_load'][sorted_position[i][0]][1] / 180) * abs(self._right_support[1] - self._loads_position[sorted_position[i][0]][1])
moment_sum_from_right_support += self._loads['point_load'][sorted_position[i][0]][0] * sin(pi * self._loads['point_load'][sorted_position[i][0]][1] / 180) * abs(self._right_support[0] - self._loads_position[sorted_position[i][0]][0])
label = Symbol(sorted_position[0][0]+"_"+sorted_position[1][0])
y2 = self._loads_position[sorted_position[1][0]][1]
x2 = self._loads_position[sorted_position[1][0]][0]
x1 = self._left_support[0]
y1 = self._left_support[1]
angle_with_horizontal = -atan((y2 - y1)/(x2 - x1))
tension = -(moment_sum_from_right_support)/(abs(self._right_support[1] - self._loads_position[sorted_position[1][0]][1])*cos(angle_with_horizontal) + abs(self._right_support[0] - self._loads_position[sorted_position[1][0]][0])*sin(angle_with_horizontal))
self._tension[label] = tension
tension_func.insert(0,(tension, x<=x2))
self._tension_func = Piecewise(*tension_func)
angle_with_horizontal = pi/2 - angle_with_horizontal
label = self._support_labels[0]
self._reaction_loads[Symbol("R_"+label+"_x")] = -sin(angle_with_horizontal) * tension
F_x += -sin(angle_with_horizontal) * tension
self._reaction_loads[Symbol("R_"+label+"_y")] = cos(angle_with_horizontal) * tension
F_y += cos(angle_with_horizontal) * tension
label = self._support_labels[1]
self._reaction_loads[Symbol("R_"+label+"_x")] = -F_x
self._reaction_loads[Symbol("R_"+label+"_y")] = -F_y
elif len(self._loads['distributed']) != 0 :
if len(args) == 0:
raise ValueError("Provide the lowest point of the cable")
lowest_x = sympify(args[0])
self._lowest_x_global = lowest_x
a = Symbol('a', positive=True)
c = Symbol('c')
# augmented matrix form of linsolve
M = Matrix(
[[(self._left_support[0]-lowest_x)**2, 1, self._left_support[1]],
[(self._right_support[0]-lowest_x)**2, 1, self._right_support[1]],
])
coefficient_solution = list(linsolve(M, (a, c)))
if len(coefficient_solution) ==0 or coefficient_solution[0][0]== 0:
raise ValueError("The lowest point is inconsistent with the supports")
A = coefficient_solution[0][0]
C = coefficient_solution[0][1] + coefficient_solution[0][0]*lowest_x**2
B = -2*coefficient_solution[0][0]*lowest_x
self._lowest_y_global = coefficient_solution[0][1]
lowest_y = self._lowest_y_global
# y = A*x**2 + B*x + C
# shifting origin to lowest point
X = Symbol('X')
Y = Symbol('Y')
Y = A*(X + lowest_x)**2 + B*(X + lowest_x) + C - lowest_y
temp_list = list(self._loads['distributed'].values())
applied_force = temp_list[0]
horizontal_force_constant = (applied_force * (self._right_support[0] - lowest_x)**2) / (2 * (self._right_support[1] - lowest_y))
self._tension.clear()
tangent_slope_to_curve = diff(Y, X)
self._tension['distributed'] = horizontal_force_constant / (cos(atan(tangent_slope_to_curve)))
label = self._support_labels[0]
self._reaction_loads[Symbol("R_"+label+"_x")] = self.tension_at(self._left_support[0]) * cos(atan(tangent_slope_to_curve.subs(X, self._left_support[0] - lowest_x)))
self._reaction_loads[Symbol("R_"+label+"_y")] = self.tension_at(self._left_support[0]) * sin(atan(tangent_slope_to_curve.subs(X, self._left_support[0] - lowest_x)))
label = self._support_labels[1]
self._reaction_loads[Symbol("R_"+label+"_x")] = self.tension_at(self._left_support[0]) * cos(atan(tangent_slope_to_curve.subs(X, self._right_support[0] - lowest_x)))
self._reaction_loads[Symbol("R_"+label+"_y")] = self.tension_at(self._left_support[0]) * sin(atan(tangent_slope_to_curve.subs(X, self._right_support[0] - lowest_x)))
def draw(self):
"""
This method is used to obtain a plot for the specified cable with its supports,
shape and loads.
Examples
========
For point loads,
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(("A", 0, 10), ("B", 10, 10))
>>> c.apply_load(-1, ('Z', 2, 7.26, 3, 270))
>>> c.apply_load(-1, ('X', 4, 6, 8, 270))
>>> c.solve()
>>> p = c.draw()
>>> p # doctest: +ELLIPSIS
Plot object containing:
[0]: cartesian line: Piecewise((10 - 1.37*x, x <= 2), (8.52 - 0.63*x, x <= 4), (2*x/3 + 10/3, x <= 10)) for x over (0.0, 10.0)
...
>>> p.show()
For uniformly distributed loads,
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c=Cable(("A", 0, 40),("B", 100, 20))
>>> c.apply_load(0, ("X", 850))
>>> c.solve(58.58)
>>> p = c.draw()
>>> p # doctest: +ELLIPSIS
Plot object containing:
[0]: cartesian line: 0.0116550116550117*(x - 58.58)**2 + 0.00447086247086247 for x over (0.0, 100.0)
[1]: cartesian line: -7.49552913752915 for x over (0.0, 100.0)
...
>>> p.show()
"""
x = Symbol("x")
annotations = []
support_rectangles = self._draw_supports()
xy_min = min(self._left_support[0],self._lowest_y_global)
xy_max = max(self._right_support[0], max(self._right_support[1],self._left_support[1]))
max_diff = xy_max - xy_min
if len(self._loads_position) != 0:
self._cable_eqn = self._draw_cable(-1)
annotations += self._draw_loads(-1)
elif len(self._loads['distributed']) != 0 :
self._cable_eqn = self._draw_cable(0)
annotations += self._draw_loads(0)
if not self._cable_eqn:
raise ValueError("solve method not called and/or values provided for loads and supports not adequate")
cab_plot = plot(*self._cable_eqn,(x,self._left_support[0],self._right_support[0]),
xlim=(xy_min-0.5*max_diff,xy_max+0.5*max_diff),
ylim=(xy_min-0.5*max_diff,xy_max+0.5*max_diff),
rectangles=support_rectangles,show= False,annotations=annotations, axis=False)
return cab_plot
def _draw_supports(self):
member_rectangles = []
xy_min = min(self._left_support[0],self._lowest_y_global)
xy_max = max(self._right_support[0], max(self._right_support[1],self._left_support[1]))
max_diff = xy_max - xy_min
supp_width = 0.075*max_diff
member_rectangles.append(
{
'xy': (self._left_support[0]-supp_width,self._left_support[1]),
'width': supp_width,
'height':supp_width,
'color':'brown',
'fill': False
}
)
member_rectangles.append(
{
'xy': (self._right_support[0],self._right_support[1]),
'width': supp_width,
'height':supp_width,
'color':'brown',
'fill': False
}
)
return member_rectangles
def _draw_cable(self,order):
xy_min = min(self._left_support[0],self._lowest_y_global)
xy_max = max(self._right_support[0], max(self._right_support[1],self._left_support[1]))
max_diff = xy_max - xy_min
if order == -1 :
x,y = symbols('x y')
line_func = []
sorted_position = sorted(self._loads_position.items(), key = lambda item : item[1][0])
for i in range(len(sorted_position)):
if(i==0):
y = ((sorted_position[i][1][1] - self._left_support[1])*(x-self._left_support[0]))/(sorted_position[i][1][0]- self._left_support[0]) + self._left_support[1]
else:
y = ((sorted_position[i][1][1] - sorted_position[i-1][1][1] )*(x-sorted_position[i-1][1][0]))/(sorted_position[i][1][0]- sorted_position[i-1][1][0]) + sorted_position[i-1][1][1]
line_func.append((y,x<=sorted_position[i][1][0]))
y = ((sorted_position[len(sorted_position)-1][1][1] - self._right_support[1])*(x-self._right_support[0]))/(sorted_position[i][1][0]- self._right_support[0]) + self._right_support[1]
line_func.append((y,x<=self._right_support[0]))
return [Piecewise(*line_func)]
elif order == 0:
x0 = self._lowest_x_global
diff_force_height = max_diff*0.075
a,c,x,y = symbols('a c x y')
parabola_eqn = a*(x-x0)**2 + c - y
points = [(self._left_support[0],self._left_support[1]),(self._right_support[0],self._right_support[1])]
equations = []
for px, py in points:
equations.append(parabola_eqn.subs({x: px, y: py}))
solution = solve(equations, (a, c))
parabola_eqn = solution[a]*(x-x0)**2 + solution[c]
return [parabola_eqn, self._lowest_y_global - diff_force_height]
def _draw_loads(self,order):
xy_min = min(self._left_support[0],self._lowest_y_global)
xy_max = max(self._right_support[0], max(self._right_support[1],self._left_support[1]))
max_diff = xy_max - xy_min
if(order==-1):
arrow_length = max_diff*0.1
force_arrows = []
for key in self._loads['point_load']:
force_arrows.append(
{
'text': '',
'xy':(self._loads_position[key][0]+arrow_length*cos(rad(self._loads['point_load'][key][1])),\
self._loads_position[key][1] + arrow_length*sin(rad(self._loads['point_load'][key][1]))),
'xytext': (self._loads_position[key][0],self._loads_position[key][1]),
'arrowprops': {'width': 1, 'headlength':3, 'headwidth':3 , 'facecolor': 'black', }
}
)
mag = self._loads['point_load'][key][0]
force_arrows.append(
{
'text':f'{mag}N',
'xy': (self._loads_position[key][0]+arrow_length*1.6*cos(rad(self._loads['point_load'][key][1])),\
self._loads_position[key][1] + arrow_length*1.6*sin(rad(self._loads['point_load'][key][1]))),
}
)
return force_arrows
elif (order == 0):
x = symbols('x')
force_arrows = []
x_val = [self._left_support[0] + ((self._right_support[0]-self._left_support[0])/10)*i for i in range(1,10)]
for i in x_val:
force_arrows.append(
{
'text':'',
'xytext':(
i,
self._cable_eqn[0].subs(x,i)
),
'xy':(
i,
self._cable_eqn[1].subs(x,i)
),
'arrowprops':{'width':1, 'headlength':3.5, 'headwidth':3.5, 'facecolor':'black'}
}
)
mag = 0
for key in self._loads['distributed']:
mag += self._loads['distributed'][key]
force_arrows.append(
{
'text':f'{mag} N/m',
'xy':((self._left_support[0]+self._right_support[0])/2,self._lowest_y_global - max_diff*0.15)
}
)
return force_arrows
def plot_tension(self):
"""
Returns the diagram/plot of the tension generated in the cable at various points.
Examples
========
For point loads,
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c = Cable(("A", 0, 10), ("B", 10, 10))
>>> c.apply_load(-1, ('Z', 2, 7.26, 3, 270))
>>> c.apply_load(-1, ('X', 4, 6, 8, 270))
>>> c.solve()
>>> p = c.plot_tension()
>>> p
Plot object containing:
[0]: cartesian line: Piecewise((8.91403453669861, x <= 2), (4.79150773600774, x <= 4), (19*sqrt(13)/10, x <= 10)) for x over (0.0, 10.0)
>>> p.show()
For uniformly distributed loads,
>>> from sympy.physics.continuum_mechanics.cable import Cable
>>> c=Cable(("A", 0, 40),("B", 100, 20))
>>> c.apply_load(0, ("X", 850))
>>> c.solve(58.58)
>>> p = c.plot_tension()
>>> p
Plot object containing:
[0]: cartesian line: 36465.0*sqrt(0.00054335718671383*X**2 + 1) for X over (0.0, 100.0)
>>> p.show()
"""
if len(self._loads_position) != 0:
x = symbols('x')
tension_plot = plot(self._tension_func, (x,self._left_support[0],self._right_support[0]), show=False)
else:
X = symbols('X')
tension_plot = plot(self._tension['distributed'], (X,self._left_support[0],self._right_support[0]), show=False)
return tension_plot
@@ -0,0 +1,61 @@
from sympy.physics.continuum_mechanics.arch import Arch
from sympy import Symbol, simplify
x = Symbol('x')
t = Symbol('t')
def test_arch_init():
a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
assert a.get_loads == {'distributed': {}, 'concentrated': {}}
assert a.reaction_force == {Symbol('R_A_x'):0, Symbol('R_A_y'):0, Symbol('R_B_x'):0, Symbol('R_B_y'):0}
assert a.supports == {'left':'hinge', 'right':'hinge'}
assert a.left_support == (0,0)
assert a.right_support == (10,0)
assert a.get_shape_eqn == 5 - ((x-5)**2)/5
a = Arch((0,0),(10,1),crown_x=6)
a.change_support_type(left_support='roller')
a.add_member(0.5)
assert a.supports == {'left':'roller', 'right':'hinge'}
assert simplify(a.get_shape_eqn) == simplify(9/5 - (x - 6)**2/20)
def test_arch_support():
a = Arch((0,0),(40,0),crown_x=20,crown_y=12)
a.apply_load(-1,'C',8,150,angle=270)
a.apply_load(0,'D',start=20,end=40,mag=-4)
a.solve()
assert abs(a.reaction_force[Symbol("R_A_x")] - 83.33333333333333) < 10e-12
assert abs(a.reaction_force[Symbol("R_B_y")] - 90.00000000000000) < 10e-12
assert abs(a.reaction_force[Symbol("R_B_x")] + 83.33333333333333) < 10e-12
assert abs(a.reaction_force[Symbol("R_A_y")] - 140.00000000000000) < 10e-12
def test_arch_member():
a = Arch((0,0),(40,0),crown_x=20,crown_y=15)
a.change_support_type(right_support='roller')
a.add_member(0)
a.apply_load(-1,'D',start=12,mag=3,angle=270)
a.apply_load(-1,'E',start=6,mag=4,angle=270)
a.apply_load(-1,'C',start=30,mag=5,angle=270)
a.solve()
assert a.reaction_force[Symbol("R_A_x")] == 0
assert abs(a.reaction_force[Symbol("R_A_y")] - 6.750000000000000) < 10e-12
assert a.reaction_force[Symbol("R_B_x")] == 0
assert abs(a.reaction_force[Symbol("R_B_y")] - 5.250000000000000) < 10e-12
def test_symbol_magnitude():
a = Arch((0,0),(16,0),crown_x=8,crown_y=5)
a.apply_load(0,'C',start=3,end=5,mag=t)
a.solve()
assert a.reaction_force[Symbol("R_A_x")] == -(4*t)/5
assert a.reaction_force[Symbol("R_A_y")] == -(3*t)/2
assert a.reaction_force[Symbol("R_B_x")] == (4*t)/5
assert a.reaction_force[Symbol("R_B_y")] == -t/2
assert a.bending_moment_at(4) == -5*t/2
def test_forces():
a = Arch((0,0),(40,0),crown_x=20,crown_y=12)
a.apply_load(-1,'C',8,150,angle=270)
a.apply_load(0,'D',start=20,end=40,mag=-4)
a.solve()
assert abs(a.axial_force_at(7.999999999999999)-149.430523405935) < 1e-12
assert abs(a.shear_force_at(7.999999999999999)-64.9227473161196) < 1e-12
@@ -0,0 +1,83 @@
from sympy.physics.continuum_mechanics.cable import Cable
from sympy.core.symbol import Symbol
def test_cable():
c = Cable(('A', 0, 10), ('B', 10, 10))
assert c.supports == {'A': [0, 10], 'B': [10, 10]}
assert c.left_support == [0, 10]
assert c.right_support == [10, 10]
assert c.loads == {'distributed': {}, 'point_load': {}}
assert c.loads_position == {}
assert c.length == 0
assert c.reaction_loads == {Symbol("R_A_x"): 0, Symbol("R_A_y"): 0, Symbol("R_B_x"): 0, Symbol("R_B_y"): 0}
# tests for change_support method
c.change_support('A', ('C', 12, 3))
assert c.supports == {'B': [10, 10], 'C': [12, 3]}
assert c.left_support == [10, 10]
assert c.right_support == [12, 3]
assert c.reaction_loads == {Symbol("R_B_x"): 0, Symbol("R_B_y"): 0, Symbol("R_C_x"): 0, Symbol("R_C_y"): 0}
c.change_support('C', ('A', 0, 10))
# tests for apply_load method for point loads
c.apply_load(-1, ('X', 2, 5, 3, 30))
c.apply_load(-1, ('Y', 5, 8, 5, 60))
assert c.loads == {'distributed': {}, 'point_load': {'X': [3, 30], 'Y': [5, 60]}}
assert c.loads_position == {'X': [2, 5], 'Y': [5, 8]}
assert c.length == 0
assert c.reaction_loads == {Symbol("R_A_x"): 0, Symbol("R_A_y"): 0, Symbol("R_B_x"): 0, Symbol("R_B_y"): 0}
# tests for remove_loads method
c.remove_loads('X')
assert c.loads == {'distributed': {}, 'point_load': {'Y': [5, 60]}}
assert c.loads_position == {'Y': [5, 8]}
assert c.length == 0
assert c.reaction_loads == {Symbol("R_A_x"): 0, Symbol("R_A_y"): 0, Symbol("R_B_x"): 0, Symbol("R_B_y"): 0}
c.remove_loads('Y')
#tests for apply_load method for distributed load
c.apply_load(0, ('Z', 9))
assert c.loads == {'distributed': {'Z': 9}, 'point_load': {}}
assert c.loads_position == {}
assert c.length == 0
assert c.reaction_loads == {Symbol("R_A_x"): 0, Symbol("R_A_y"): 0, Symbol("R_B_x"): 0, Symbol("R_B_y"): 0}
# tests for apply_length method
c.apply_length(20)
assert c.length == 20
del c
# tests for solve method
# for point loads
c = Cable(("A", 0, 10), ("B", 5.5, 8))
c.apply_load(-1, ('Z', 2, 7.26, 3, 270))
c.apply_load(-1, ('X', 4, 6, 8, 270))
c.solve()
#assert c.tension == {Symbol("Z_X"): 4.79150773600774, Symbol("X_B"): 6.78571428571429, Symbol("A_Z"): 6.89488895397307}
assert abs(c.tension[Symbol("A_Z")] - 6.89488895397307) < 10e-12
assert abs(c.tension[Symbol("Z_X")] - 4.79150773600774) < 10e-12
assert abs(c.tension[Symbol("X_B")] - 6.78571428571429) < 10e-12
#assert c.reaction_loads == {Symbol("R_A_x"): -4.06504065040650, Symbol("R_A_y"): 5.56910569105691, Symbol("R_B_x"): 4.06504065040650, Symbol("R_B_y"): 5.43089430894309}
assert abs(c.reaction_loads[Symbol("R_A_x")] + 4.06504065040650) < 10e-12
assert abs(c.reaction_loads[Symbol("R_A_y")] - 5.56910569105691) < 10e-12
assert abs(c.reaction_loads[Symbol("R_B_x")] - 4.06504065040650) < 10e-12
assert abs(c.reaction_loads[Symbol("R_B_y")] - 5.43089430894309) < 10e-12
assert abs(c.length - 8.25609584845190) < 10e-12
del c
# tests for solve method
# for distributed loads
c=Cable(("A", 0, 40),("B", 100, 20))
c.apply_load(0, ("X", 850))
c.solve(58.58, 0)
# assert c.tension['distributed'] == 36456.8485*sqrt(0.000543529004799705*(X + 0.00135624381275735)**2 + 1)
assert abs(c.tension_at(0) - 61717.4130533677) < 10e-11
assert abs(c.tension_at(40) - 39738.0809048449) < 10e-11
assert abs(c.reaction_loads[Symbol("R_A_x")] - 36465.0000000000) < 10e-11
assert abs(c.reaction_loads[Symbol("R_A_y")] + 49793.0000000000) < 10e-11
assert abs(c.reaction_loads[Symbol("R_B_x")] - 44399.9537590861) < 10e-11
assert abs(c.reaction_loads[Symbol("R_B_y")] - 42868.2071025955 ) < 10e-11
@@ -0,0 +1,100 @@
from sympy.core.symbol import Symbol, symbols
from sympy.physics.continuum_mechanics.truss import Truss
from sympy import sqrt
def test_truss():
A = Symbol('A')
B = Symbol('B')
C = Symbol('C')
AB, BC, AC = symbols('AB, BC, AC')
P = Symbol('P')
t = Truss()
assert t.nodes == []
assert t.node_labels == []
assert t.node_positions == []
assert t.members == {}
assert t.loads == {}
assert t.supports == {}
assert t.reaction_loads == {}
assert t.internal_forces == {}
# testing the add_node method
t.add_node((A, 0, 0), (B, 2, 2), (C, 3, 0))
assert t.nodes == [(A, 0, 0), (B, 2, 2), (C, 3, 0)]
assert t.node_labels == [A, B, C]
assert t.node_positions == [(0, 0), (2, 2), (3, 0)]
assert t.loads == {}
assert t.supports == {}
assert t.reaction_loads == {}
# testing the remove_node method
t.remove_node(C)
assert t.nodes == [(A, 0, 0), (B, 2, 2)]
assert t.node_labels == [A, B]
assert t.node_positions == [(0, 0), (2, 2)]
assert t.loads == {}
assert t.supports == {}
t.add_node((C, 3, 0))
# testing the add_member method
t.add_member((AB, A, B), (BC, B, C), (AC, A, C))
assert t.members == {AB: [A, B], BC: [B, C], AC: [A, C]}
assert t.internal_forces == {AB: 0, BC: 0, AC: 0}
# testing the remove_member method
t.remove_member(BC)
assert t.members == {AB: [A, B], AC: [A, C]}
assert t.internal_forces == {AB: 0, AC: 0}
t.add_member((BC, B, C))
D, CD = symbols('D, CD')
# testing the change_label methods
t.change_node_label((B, D))
assert t.nodes == [(A, 0, 0), (D, 2, 2), (C, 3, 0)]
assert t.node_labels == [A, D, C]
assert t.loads == {}
assert t.supports == {}
assert t.members == {AB: [A, D], BC: [D, C], AC: [A, C]}
t.change_member_label((BC, CD))
assert t.members == {AB: [A, D], CD: [D, C], AC: [A, C]}
assert t.internal_forces == {AB: 0, CD: 0, AC: 0}
# testing the apply_load method
t.apply_load((A, P, 90), (A, P/4, 90), (A, 2*P,45), (D, P/2, 90))
assert t.loads == {A: [[P, 90], [P/4, 90], [2*P, 45]], D: [[P/2, 90]]}
assert t.loads[A] == [[P, 90], [P/4, 90], [2*P, 45]]
# testing the remove_load method
t.remove_load((A, P/4, 90))
assert t.loads == {A: [[P, 90], [2*P, 45]], D: [[P/2, 90]]}
assert t.loads[A] == [[P, 90], [2*P, 45]]
# testing the apply_support method
t.apply_support((A, "pinned"), (D, "roller"))
assert t.supports == {A: 'pinned', D: 'roller'}
assert t.reaction_loads == {}
assert t.loads == {A: [[P, 90], [2*P, 45], [Symbol('R_A_x'), 0], [Symbol('R_A_y'), 90]], D: [[P/2, 90], [Symbol('R_D_y'), 90]]}
# testing the remove_support method
t.remove_support(A)
assert t.supports == {D: 'roller'}
assert t.reaction_loads == {}
assert t.loads == {A: [[P, 90], [2*P, 45]], D: [[P/2, 90], [Symbol('R_D_y'), 90]]}
t.apply_support((A, "pinned"))
# testing the solve method
t.solve()
assert t.reaction_loads['R_A_x'] == -sqrt(2)*P
assert t.reaction_loads['R_A_y'] == -sqrt(2)*P - P
assert t.reaction_loads['R_D_y'] == -P/2
assert t.internal_forces[AB]/P == 0
assert t.internal_forces[CD] == 0
assert t.internal_forces[AC] == 0
@@ -0,0 +1,17 @@
from .lti import (TransferFunction, PIDController, Series, MIMOSeries, Parallel, MIMOParallel,
Feedback, MIMOFeedback, TransferFunctionMatrix, StateSpace, gbt, bilinear, forward_diff,
backward_diff, phase_margin, gain_margin)
from .control_plots import (pole_zero_numerical_data, pole_zero_plot, step_response_numerical_data,
step_response_plot, impulse_response_numerical_data, impulse_response_plot, ramp_response_numerical_data,
ramp_response_plot, bode_magnitude_numerical_data, bode_phase_numerical_data, bode_magnitude_plot,
bode_phase_plot, bode_plot, nyquist_plot_expr, nyquist_plot, nichols_plot_expr, nichols_plot)
__all__ = ['TransferFunction', 'PIDController', 'Series', 'MIMOSeries', 'Parallel',
'MIMOParallel', 'Feedback', 'MIMOFeedback', 'TransferFunctionMatrix', 'StateSpace',
'gbt', 'bilinear', 'forward_diff', 'backward_diff', 'phase_margin', 'gain_margin',
'pole_zero_numerical_data', 'pole_zero_plot', 'step_response_numerical_data',
'step_response_plot', 'impulse_response_numerical_data', 'impulse_response_plot',
'ramp_response_numerical_data', 'ramp_response_plot',
'bode_magnitude_numerical_data', 'bode_phase_numerical_data',
'bode_magnitude_plot', 'bode_phase_plot', 'bode_plot', 'nyquist_plot_expr', 'nyquist_plot',
'nichols_plot_expr', 'nichols_plot']
@@ -0,0 +1,332 @@
from math import isclose
from sympy.core.numbers import I, all_close
from sympy.core.symbol import Dummy
from sympy.functions.elementary.complexes import (Abs, arg)
from sympy.functions.elementary.exponential import log
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.abc import s, p, a
from sympy import pi
from sympy.external import import_module
from sympy.physics.control.control_plots import \
(pole_zero_numerical_data, pole_zero_plot, step_response_numerical_data,
step_response_plot, impulse_response_numerical_data,
impulse_response_plot, ramp_response_numerical_data,
ramp_response_plot, bode_magnitude_numerical_data,
bode_phase_numerical_data, bode_plot, nyquist_plot_expr,
nichols_plot_expr)
from sympy.physics.control.lti import (TransferFunction,
Series, Parallel, TransferFunctionMatrix)
from sympy.testing.pytest import raises, skip
matplotlib = import_module(
'matplotlib', import_kwargs={'fromlist': ['pyplot']},
catch=(RuntimeError,))
numpy = import_module('numpy')
tf1 = TransferFunction(1, p**2 + 0.5*p + 2, p)
tf2 = TransferFunction(p, 6*p**2 + 3*p + 1, p)
tf3 = TransferFunction(p, p**3 - 1, p)
tf4 = TransferFunction(10, p**3, p)
tf5 = TransferFunction(5, s**2 + 2*s + 10, s)
tf6 = TransferFunction(1, 1, s)
tf7 = TransferFunction(4*s*3 + 9*s**2 + 0.1*s + 11, 8*s**6 + 9*s**4 + 11, s)
tf8 = TransferFunction(5, s**2 + (2+I)*s + 10, s)
ser1 = Series(tf4, TransferFunction(1, p - 5, p))
ser2 = Series(tf3, TransferFunction(p, p + 2, p))
par1 = Parallel(tf1, tf2)
def _to_tuple(a, b):
return tuple(a), tuple(b)
def _trim_tuple(a, b):
a, b = _to_tuple(a, b)
return tuple(a[0: 2] + a[len(a)//2 : len(a)//2 + 1] + a[-2:]), \
tuple(b[0: 2] + b[len(b)//2 : len(b)//2 + 1] + b[-2:])
def y_coordinate_equality(plot_data_func, evalf_func, system):
"""Checks whether the y-coordinate value of the plotted
data point is equal to the value of the function at a
particular x."""
x, y = plot_data_func(system)
x, y = _trim_tuple(x, y)
y_exp = tuple(evalf_func(system, x_i) for x_i in x)
return all(Abs(y_exp_i - y_i) < 1e-8 for y_exp_i, y_i in zip(y_exp, y))
def test_errors():
if not matplotlib:
skip("Matplotlib not the default backend")
# Invalid `system` check
tfm = TransferFunctionMatrix([[tf6, tf5], [tf5, tf6]])
expr = 1/(s**2 - 1)
raises(NotImplementedError, lambda: pole_zero_plot(tfm))
raises(NotImplementedError, lambda: pole_zero_numerical_data(expr))
raises(NotImplementedError, lambda: impulse_response_plot(expr))
raises(NotImplementedError, lambda: impulse_response_numerical_data(tfm))
raises(NotImplementedError, lambda: step_response_plot(tfm))
raises(NotImplementedError, lambda: step_response_numerical_data(expr))
raises(NotImplementedError, lambda: ramp_response_plot(expr))
raises(NotImplementedError, lambda: ramp_response_numerical_data(tfm))
raises(NotImplementedError, lambda: bode_plot(tfm))
# More than 1 variables
tf_a = TransferFunction(a, s + 1, s)
raises(ValueError, lambda: pole_zero_plot(tf_a))
raises(ValueError, lambda: pole_zero_numerical_data(tf_a))
raises(ValueError, lambda: impulse_response_plot(tf_a))
raises(ValueError, lambda: impulse_response_numerical_data(tf_a))
raises(ValueError, lambda: step_response_plot(tf_a))
raises(ValueError, lambda: step_response_numerical_data(tf_a))
raises(ValueError, lambda: ramp_response_plot(tf_a))
raises(ValueError, lambda: ramp_response_numerical_data(tf_a))
raises(ValueError, lambda: bode_plot(tf_a))
# lower_limit > 0 for response plots
raises(ValueError, lambda: impulse_response_plot(tf1, lower_limit=-1))
raises(ValueError, lambda: step_response_plot(tf1, lower_limit=-0.1))
raises(ValueError, lambda: ramp_response_plot(tf1, lower_limit=-4/3))
# slope in ramp_response_plot() is negative
raises(ValueError, lambda: ramp_response_plot(tf1, slope=-0.1))
# incorrect frequency or phase unit
raises(ValueError, lambda: bode_plot(tf1,freq_unit = 'hz'))
raises(ValueError, lambda: bode_plot(tf1,phase_unit = 'degree'))
def test_pole_zero():
def pz_tester(sys, expected_value):
_z, _p = pole_zero_numerical_data(sys)
z_check = all_close(_z, expected_value[0])
p_check = all_close(_p, expected_value[1])
return p_check and z_check
exp1 = [[], [-0.24999999999999994-1.3919410907075054j, -0.24999999999999994+1.3919410907075054j]]
exp2 = [[0.0], [-0.25-0.3227486121839514j, -0.25+0.3227486121839514j]]
exp3 = [[0.0], [0.9999999999999998+0j, -0.5000000000000004-0.8660254037844395j,
-0.5000000000000004+0.8660254037844395j]]
exp4 = [[], [0.0, 0.0, 0.0, 5.0]]
exp5 = [[-5.645751311064592, -0.5000000000000008, -0.3542486889354093],
[-0.24999999999999986-0.322748612183951348j,
-0.2499999999999998+0.32274861218395134j,
-0.24999999999999986-1.3919410907075052j,
-0.2499999999999998+1.3919410907075052j]]
exp6 = [[], [-1.1641600331447917-3.545808351896439j,
-0.8358399668552097+2.5458083518964383j]]
assert pz_tester(tf1, exp1)
assert pz_tester(tf2, exp2)
assert pz_tester(tf3, exp3)
assert pz_tester(ser1, exp4)
assert pz_tester(par1, exp5)
assert pz_tester(tf8, exp6)
def test_bode():
if not numpy:
skip("NumPy is required for this test")
def bode_phase_evalf(system, point):
expr = system.to_expr()
_w = Dummy("w", real=True)
w_expr = expr.subs({system.var: I*_w})
return arg(w_expr).subs({_w: point}).evalf()
def bode_mag_evalf(system, point):
expr = system.to_expr()
_w = Dummy("w", real=True)
w_expr = expr.subs({system.var: I*_w})
return 20*log(Abs(w_expr), 10).subs({_w: point}).evalf()
def test_bode_data(sys):
return y_coordinate_equality(bode_magnitude_numerical_data, bode_mag_evalf, sys) \
and y_coordinate_equality(bode_phase_numerical_data, bode_phase_evalf, sys)
assert test_bode_data(tf1)
assert test_bode_data(tf2)
assert test_bode_data(tf3)
assert test_bode_data(tf4)
assert test_bode_data(tf5)
def check_point_accuracy(a, b):
return all(isclose(*_, rel_tol=1e-1, abs_tol=1e-6
) for _ in zip(a, b))
def test_impulse_response():
if not numpy:
skip("NumPy is required for this test")
def impulse_res_tester(sys, expected_value):
x, y = _to_tuple(*impulse_response_numerical_data(sys,
adaptive=False, n=10))
x_check = check_point_accuracy(x, expected_value[0])
y_check = check_point_accuracy(y, expected_value[1])
return x_check and y_check
exp1 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(0.0, 0.544019738507865, 0.01993849743234938, -0.31140243360893216, -0.022852779906491996, 0.1778306498155759,
0.01962941084328499, -0.1013115194573652, -0.014975541213105696, 0.0575789724730714))
exp2 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.1666666675, 0.08389223412935855,
0.02338051973475047, -0.014966807776379383, -0.034645954223054234, -0.040560075735512804,
-0.037658628907103885, -0.030149507719590022, -0.021162090730736834, -0.012721292737437523))
exp3 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (4.369893391586999e-09, 1.1750333000630964,
3.2922404058312473, 9.432290008148343, 28.37098083007151, 86.18577464367974, 261.90356653762115,
795.6538758627842, 2416.9920942096983, 7342.159505206647))
exp4 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 6.17283950617284, 24.69135802469136,
55.555555555555564, 98.76543209876544, 154.320987654321, 222.22222222222226, 302.46913580246917,
395.0617283950618, 500.0))
exp5 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, -0.10455606138085417,
0.06757671513476461, -0.03234567568833768, 0.013582514927757873, -0.005273419510705473,
0.0019364083003354075, -0.000680070134067832, 0.00022969845960406913, -7.476094359583917e-05))
exp6 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(-6.016699583000218e-09, 0.35039802056107394, 3.3728423827689884, 12.119846079276684,
25.86101014293389, 29.352480635282088, -30.49475907497664, -273.8717189554019, -863.2381702029659,
-1747.0262164682233))
exp7 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335,
4.444444444444445, 5.555555555555555, 6.666666666666667, 7.777777777777779,
8.88888888888889, 10.0), (0.0, 18.934638095560974, 5346.93244680907, 1384609.8718249386,
358161126.65801865, 92645770015.70108, 23964739753087.42, 6198974342083139.0, 1.603492601616059e+18,
4.147764422869658e+20))
assert impulse_res_tester(tf1, exp1)
assert impulse_res_tester(tf2, exp2)
assert impulse_res_tester(tf3, exp3)
assert impulse_res_tester(tf4, exp4)
assert impulse_res_tester(tf5, exp5)
assert impulse_res_tester(tf7, exp6)
assert impulse_res_tester(ser1, exp7)
def test_step_response():
if not numpy:
skip("NumPy is required for this test")
def step_res_tester(sys, expected_value):
x, y = _to_tuple(*step_response_numerical_data(sys,
adaptive=False, n=10))
x_check = check_point_accuracy(x, expected_value[0])
y_check = check_point_accuracy(y, expected_value[1])
return x_check and y_check
exp1 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(-1.9193285738516863e-08, 0.42283495488246126, 0.7840485977945262, 0.5546841805655717,
0.33903033806932087, 0.4627251747410237, 0.5909907598988051, 0.5247213989553071,
0.4486997874319281, 0.4839358435839171))
exp2 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(0.0, 0.13728409095645816, 0.19474559355325086, 0.1974909129243011, 0.16841657696573073,
0.12559777736159378, 0.08153828016664713, 0.04360471317348958, 0.015072994568868221,
-0.003636420058445484))
exp3 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(0.0, 0.6314542141914303, 2.9356520038101035, 9.37731009663807, 28.452300356688376,
86.25721933273988, 261.9236645044672, 795.6435410577224, 2416.9786984578764, 7342.154119725917))
exp4 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(0.0, 2.286236899862826, 18.28989519890261, 61.72839629629631, 146.31916159122088, 285.7796124828532,
493.8271703703705, 784.1792566529494, 1170.553292729767, 1666.6667))
exp5 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(-3.999999997894577e-09, 0.6720357068882895, 0.4429938256137113, 0.5182010838004518,
0.4944139147159695, 0.5016379853883338, 0.4995466896527733, 0.5001154784851325,
0.49997448824584123, 0.5000039745919259))
exp6 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(-1.5433688493882158e-09, 0.3428705539937336, 1.1253619102202777, 3.1849962651016517,
9.47532757182671, 28.727231099148135, 87.29426924860557, 265.2138681048606, 805.6636260007757,
2447.387582370878))
assert step_res_tester(tf1, exp1)
assert step_res_tester(tf2, exp2)
assert step_res_tester(tf3, exp3)
assert step_res_tester(tf4, exp4)
assert step_res_tester(tf5, exp5)
assert step_res_tester(ser2, exp6)
def test_ramp_response():
if not numpy:
skip("NumPy is required for this test")
def ramp_res_tester(sys, num_points, expected_value, slope=1):
x, y = _to_tuple(*ramp_response_numerical_data(sys,
slope=slope, adaptive=False, n=num_points))
x_check = check_point_accuracy(x, expected_value[0])
y_check = check_point_accuracy(y, expected_value[1])
return x_check and y_check
exp1 = ((0.0, 2.0, 4.0, 6.0, 8.0, 10.0), (0.0, 0.7324667795033895, 1.9909720978650398,
2.7956587704217783, 3.9224897567931514, 4.85022655284895))
exp2 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
(2.4360213402019326e-08, 0.10175320182493253, 0.33057612497658406, 0.5967937263298935,
0.8431511866718248, 1.0398805391471613, 1.1776043125035738, 1.2600994825747305, 1.2981042689274653,
1.304684417610106))
exp3 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (-3.9329040468771836e-08,
0.34686634635794555, 2.9998828170537903, 12.33303690737476, 40.993913948137795, 127.84145222317912,
391.41713691996, 1192.0006858708389, 3623.9808672503405, 11011.728034546572))
exp4 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 1.9051973784484078, 30.483158055174524,
154.32098765432104, 487.7305288827924, 1190.7483615302544, 2469.1358024691367, 4574.3789056546275,
7803.688462124678, 12500.0))
exp5 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 3.8844361856975635, 9.141792069209865,
14.096349157657231, 19.09783068994694, 24.10179770390321, 29.09907319114121, 34.10040420185154,
39.09983919254265, 44.10006013058409))
exp6 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 1.1111111111111112, 2.2222222222222223,
3.3333333333333335, 4.444444444444445, 5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0))
assert ramp_res_tester(tf1, 6, exp1)
assert ramp_res_tester(tf2, 10, exp2, 1.2)
assert ramp_res_tester(tf3, 10, exp3, 1.5)
assert ramp_res_tester(tf4, 10, exp4, 3)
assert ramp_res_tester(tf5, 10, exp5, 9)
assert ramp_res_tester(tf6, 10, exp6)
def test_nyquist_plot_expr():
r1, i1, w1 = nyquist_plot_expr(tf1)
r2, i2, w2 = nyquist_plot_expr(tf2)
r3, i3, w3 = nyquist_plot_expr(tf3)
r4, i4, w4 = nyquist_plot_expr(tf4)
assert r1 == (2 - w1**2)/(0.25*w1**2 + (2 - w1**2)**2)
assert i1 == -0.5*w1/(0.25*w1**2 + (2 - w1**2)**2)
assert r2 == 3*w2**2/(9*w2**2 + (1 - 6*w2**2)**2)
assert i2 == w2*(1 - 6*w2**2)/(9*w2**2 + (1 - 6*w2**2)**2)
assert r3 == -w3**4/(w3**6 + 1)
assert i3 == -w3/(w3**6 + 1)
assert r4 == 0
assert i4 == 10/w4**3
def test_nichols_expr():
m1, p1, w1 = nichols_plot_expr(tf1)
m2, p2, w2 = nichols_plot_expr(tf2)
m3, p3, w3 = nichols_plot_expr(tf3)
m4, p4, w4 = nichols_plot_expr(tf4)
assert m1 == 20*log(1/sqrt(w1**4 - 3.75*w1**2 + 4))/log(10)
assert p1 == 180*arg(1/(-w1**2 + 0.5*w1*I + 2))/pi
assert m2 == 20*log(Abs(w2)/sqrt(36*w2**4 - 3*w2**2 + 1))/log(10)
assert p2 == 180*arg(w2*I/(-6*w2**2 + 3*w2*I + 1))/pi
assert m3 == 20*log(Abs(w3)/sqrt(w3**6 + 1))/log(10)
assert p3 == 180*arg(-w3*I/(w3**3*I + 1))/pi
assert m4 == 20*log(10/(w4**2*Abs(w4)))/log(10)
assert p4 == 180*arg(I/w4**3)/pi
@@ -0,0 +1,716 @@
"""
Module to handle gamma matrices expressed as tensor objects.
Examples
========
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex
>>> from sympy.tensor.tensor import tensor_indices
>>> i = tensor_indices('i', LorentzIndex)
>>> G(i)
GammaMatrix(i)
Note that there is already an instance of GammaMatrixHead in four dimensions:
GammaMatrix, which is simply declare as
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix
>>> from sympy.tensor.tensor import tensor_indices
>>> i = tensor_indices('i', LorentzIndex)
>>> GammaMatrix(i)
GammaMatrix(i)
To access the metric tensor
>>> LorentzIndex.metric
metric(LorentzIndex,LorentzIndex)
"""
from sympy.core.mul import Mul
from sympy.core.singleton import S
from sympy.matrices.dense import eye
from sympy.matrices.expressions.trace import trace
from sympy.tensor.tensor import TensorIndexType, TensorIndex,\
TensMul, TensAdd, tensor_mul, Tensor, TensorHead, TensorSymmetry
# DiracSpinorIndex = TensorIndexType('DiracSpinorIndex', dim=4, dummy_name="S")
LorentzIndex = TensorIndexType('LorentzIndex', dim=4, dummy_name="L")
GammaMatrix = TensorHead("GammaMatrix", [LorentzIndex],
TensorSymmetry.no_symmetry(1), comm=None)
def extract_type_tens(expression, component):
"""
Extract from a ``TensExpr`` all tensors with `component`.
Returns two tensor expressions:
* the first contains all ``Tensor`` of having `component`.
* the second contains all remaining.
"""
if isinstance(expression, Tensor):
sp = [expression]
elif isinstance(expression, TensMul):
sp = expression.args
else:
raise ValueError('wrong type')
# Collect all gamma matrices of the same dimension
new_expr = S.One
residual_expr = S.One
for i in sp:
if isinstance(i, Tensor) and i.component == component:
new_expr *= i
else:
residual_expr *= i
return new_expr, residual_expr
def simplify_gamma_expression(expression):
extracted_expr, residual_expr = extract_type_tens(expression, GammaMatrix)
res_expr = _simplify_single_line(extracted_expr)
return res_expr * residual_expr
def simplify_gpgp(ex, sort=True):
"""
simplify products ``G(i)*p(-i)*G(j)*p(-j) -> p(i)*p(-i)``
Examples
========
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
LorentzIndex, simplify_gpgp
>>> from sympy.tensor.tensor import tensor_indices, tensor_heads
>>> p, q = tensor_heads('p, q', [LorentzIndex])
>>> i0,i1,i2,i3,i4,i5 = tensor_indices('i0:6', LorentzIndex)
>>> ps = p(i0)*G(-i0)
>>> qs = q(i0)*G(-i0)
>>> simplify_gpgp(ps*qs*qs)
GammaMatrix(-L_0)*p(L_0)*q(L_1)*q(-L_1)
"""
def _simplify_gpgp(ex):
components = ex.components
a = []
comp_map = []
for i, comp in enumerate(components):
comp_map.extend([i]*comp.rank)
dum = [(i[0], i[1], comp_map[i[0]], comp_map[i[1]]) for i in ex.dum]
for i in range(len(components)):
if components[i] != GammaMatrix:
continue
for dx in dum:
if dx[2] == i:
p_pos1 = dx[3]
elif dx[3] == i:
p_pos1 = dx[2]
else:
continue
comp1 = components[p_pos1]
if comp1.comm == 0 and comp1.rank == 1:
a.append((i, p_pos1))
if not a:
return ex
elim = set()
tv = []
hit = True
coeff = S.One
ta = None
while hit:
hit = False
for i, ai in enumerate(a[:-1]):
if ai[0] in elim:
continue
if ai[0] != a[i + 1][0] - 1:
continue
if components[ai[1]] != components[a[i + 1][1]]:
continue
elim.add(ai[0])
elim.add(ai[1])
elim.add(a[i + 1][0])
elim.add(a[i + 1][1])
if not ta:
ta = ex.split()
mu = TensorIndex('mu', LorentzIndex)
hit = True
if i == 0:
coeff = ex.coeff
tx = components[ai[1]](mu)*components[ai[1]](-mu)
if len(a) == 2:
tx *= 4 # eye(4)
tv.append(tx)
break
if tv:
a = [x for j, x in enumerate(ta) if j not in elim]
a.extend(tv)
t = tensor_mul(*a)*coeff
# t = t.replace(lambda x: x.is_Matrix, lambda x: 1)
return t
else:
return ex
if sort:
ex = ex.sorted_components()
# this would be better off with pattern matching
while 1:
t = _simplify_gpgp(ex)
if t != ex:
ex = t
else:
return t
def gamma_trace(t):
"""
trace of a single line of gamma matrices
Examples
========
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
gamma_trace, LorentzIndex
>>> from sympy.tensor.tensor import tensor_indices, tensor_heads
>>> p, q = tensor_heads('p, q', [LorentzIndex])
>>> i0,i1,i2,i3,i4,i5 = tensor_indices('i0:6', LorentzIndex)
>>> ps = p(i0)*G(-i0)
>>> qs = q(i0)*G(-i0)
>>> gamma_trace(G(i0)*G(i1))
4*metric(i0, i1)
>>> gamma_trace(ps*ps) - 4*p(i0)*p(-i0)
0
>>> gamma_trace(ps*qs + ps*ps) - 4*p(i0)*p(-i0) - 4*p(i0)*q(-i0)
0
"""
if isinstance(t, TensAdd):
res = TensAdd(*[gamma_trace(x) for x in t.args])
return res
t = _simplify_single_line(t)
res = _trace_single_line(t)
return res
def _simplify_single_line(expression):
"""
Simplify single-line product of gamma matrices.
Examples
========
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
LorentzIndex, _simplify_single_line
>>> from sympy.tensor.tensor import tensor_indices, TensorHead
>>> p = TensorHead('p', [LorentzIndex])
>>> i0,i1 = tensor_indices('i0:2', LorentzIndex)
>>> _simplify_single_line(G(i0)*G(i1)*p(-i1)*G(-i0)) + 2*G(i0)*p(-i0)
0
"""
t1, t2 = extract_type_tens(expression, GammaMatrix)
if t1 != 1:
t1 = kahane_simplify(t1)
res = t1*t2
return res
def _trace_single_line(t):
"""
Evaluate the trace of a single gamma matrix line inside a ``TensExpr``.
Notes
=====
If there are ``DiracSpinorIndex.auto_left`` and ``DiracSpinorIndex.auto_right``
indices trace over them; otherwise traces are not implied (explain)
Examples
========
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
LorentzIndex, _trace_single_line
>>> from sympy.tensor.tensor import tensor_indices, TensorHead
>>> p = TensorHead('p', [LorentzIndex])
>>> i0,i1,i2,i3,i4,i5 = tensor_indices('i0:6', LorentzIndex)
>>> _trace_single_line(G(i0)*G(i1))
4*metric(i0, i1)
>>> _trace_single_line(G(i0)*p(-i0)*G(i1)*p(-i1)) - 4*p(i0)*p(-i0)
0
"""
def _trace_single_line1(t):
t = t.sorted_components()
components = t.components
ncomps = len(components)
g = LorentzIndex.metric
# gamma matirices are in a[i:j]
hit = 0
for i in range(ncomps):
if components[i] == GammaMatrix:
hit = 1
break
for j in range(i + hit, ncomps):
if components[j] != GammaMatrix:
break
else:
j = ncomps
numG = j - i
if numG == 0:
tcoeff = t.coeff
return t.nocoeff if tcoeff else t
if numG % 2 == 1:
return TensMul.from_data(S.Zero, [], [], [])
elif numG > 4:
# find the open matrix indices and connect them:
a = t.split()
ind1 = a[i].get_indices()[0]
ind2 = a[i + 1].get_indices()[0]
aa = a[:i] + a[i + 2:]
t1 = tensor_mul(*aa)*g(ind1, ind2)
t1 = t1.contract_metric(g)
args = [t1]
sign = 1
for k in range(i + 2, j):
sign = -sign
ind2 = a[k].get_indices()[0]
aa = a[:i] + a[i + 1:k] + a[k + 1:]
t2 = sign*tensor_mul(*aa)*g(ind1, ind2)
t2 = t2.contract_metric(g)
t2 = simplify_gpgp(t2, False)
args.append(t2)
t3 = TensAdd(*args)
t3 = _trace_single_line(t3)
return t3
else:
a = t.split()
t1 = _gamma_trace1(*a[i:j])
a2 = a[:i] + a[j:]
t2 = tensor_mul(*a2)
t3 = t1*t2
if not t3:
return t3
t3 = t3.contract_metric(g)
return t3
t = t.expand()
if isinstance(t, TensAdd):
a = [_trace_single_line1(x)*x.coeff for x in t.args]
return TensAdd(*a)
elif isinstance(t, (Tensor, TensMul)):
r = t.coeff*_trace_single_line1(t)
return r
else:
return trace(t)
def _gamma_trace1(*a):
gctr = 4 # FIXME specific for d=4
g = LorentzIndex.metric
if not a:
return gctr
n = len(a)
if n%2 == 1:
#return TensMul.from_data(S.Zero, [], [], [])
return S.Zero
if n == 2:
ind0 = a[0].get_indices()[0]
ind1 = a[1].get_indices()[0]
return gctr*g(ind0, ind1)
if n == 4:
ind0 = a[0].get_indices()[0]
ind1 = a[1].get_indices()[0]
ind2 = a[2].get_indices()[0]
ind3 = a[3].get_indices()[0]
return gctr*(g(ind0, ind1)*g(ind2, ind3) - \
g(ind0, ind2)*g(ind1, ind3) + g(ind0, ind3)*g(ind1, ind2))
def kahane_simplify(expression):
r"""
This function cancels contracted elements in a product of four
dimensional gamma matrices, resulting in an expression equal to the given
one, without the contracted gamma matrices.
Parameters
==========
`expression` the tensor expression containing the gamma matrices to simplify.
Notes
=====
If spinor indices are given, the matrices must be given in
the order given in the product.
Algorithm
=========
The idea behind the algorithm is to use some well-known identities,
i.e., for contractions enclosing an even number of `\gamma` matrices
`\gamma^\mu \gamma_{a_1} \cdots \gamma_{a_{2N}} \gamma_\mu = 2 (\gamma_{a_{2N}} \gamma_{a_1} \cdots \gamma_{a_{2N-1}} + \gamma_{a_{2N-1}} \cdots \gamma_{a_1} \gamma_{a_{2N}} )`
for an odd number of `\gamma` matrices
`\gamma^\mu \gamma_{a_1} \cdots \gamma_{a_{2N+1}} \gamma_\mu = -2 \gamma_{a_{2N+1}} \gamma_{a_{2N}} \cdots \gamma_{a_{1}}`
Instead of repeatedly applying these identities to cancel out all contracted indices,
it is possible to recognize the links that would result from such an operation,
the problem is thus reduced to a simple rearrangement of free gamma matrices.
Examples
========
When using, always remember that the original expression coefficient
has to be handled separately
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex
>>> from sympy.physics.hep.gamma_matrices import kahane_simplify
>>> from sympy.tensor.tensor import tensor_indices
>>> i0, i1, i2 = tensor_indices('i0:3', LorentzIndex)
>>> ta = G(i0)*G(-i0)
>>> kahane_simplify(ta)
Matrix([
[4, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 4, 0],
[0, 0, 0, 4]])
>>> tb = G(i0)*G(i1)*G(-i0)
>>> kahane_simplify(tb)
-2*GammaMatrix(i1)
>>> t = G(i0)*G(-i0)
>>> kahane_simplify(t)
Matrix([
[4, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 4, 0],
[0, 0, 0, 4]])
>>> t = G(i0)*G(-i0)
>>> kahane_simplify(t)
Matrix([
[4, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 4, 0],
[0, 0, 0, 4]])
If there are no contractions, the same expression is returned
>>> tc = G(i0)*G(i1)
>>> kahane_simplify(tc)
GammaMatrix(i0)*GammaMatrix(i1)
References
==========
[1] Algorithm for Reducing Contracted Products of gamma Matrices,
Joseph Kahane, Journal of Mathematical Physics, Vol. 9, No. 10, October 1968.
"""
if isinstance(expression, Mul):
return expression
if isinstance(expression, TensAdd):
return TensAdd(*[kahane_simplify(arg) for arg in expression.args])
if isinstance(expression, Tensor):
return expression
assert isinstance(expression, TensMul)
gammas = expression.args
for gamma in gammas:
assert gamma.component == GammaMatrix
free = expression.free
# spinor_free = [_ for _ in expression.free_in_args if _[1] != 0]
# if len(spinor_free) == 2:
# spinor_free.sort(key=lambda x: x[2])
# assert spinor_free[0][1] == 1 and spinor_free[-1][1] == 2
# assert spinor_free[0][2] == 0
# elif spinor_free:
# raise ValueError('spinor indices do not match')
dum = []
for dum_pair in expression.dum:
if expression.index_types[dum_pair[0]] == LorentzIndex:
dum.append((dum_pair[0], dum_pair[1]))
dum = sorted(dum)
if len(dum) == 0: # or GammaMatrixHead:
# no contractions in `expression`, just return it.
return expression
# find the `first_dum_pos`, i.e. the position of the first contracted
# gamma matrix, Kahane's algorithm as described in his paper requires the
# gamma matrix expression to start with a contracted gamma matrix, this is
# a workaround which ignores possible initial free indices, and re-adds
# them later.
first_dum_pos = min(map(min, dum))
# for p1, p2, a1, a2 in expression.dum_in_args:
# if p1 != 0 or p2 != 0:
# # only Lorentz indices, skip Dirac indices:
# continue
# first_dum_pos = min(p1, p2)
# break
total_number = len(free) + len(dum)*2
number_of_contractions = len(dum)
free_pos = [None]*total_number
for i in free:
free_pos[i[1]] = i[0]
# `index_is_free` is a list of booleans, to identify index position
# and whether that index is free or dummy.
index_is_free = [False]*total_number
for i, indx in enumerate(free):
index_is_free[indx[1]] = True
# `links` is a dictionary containing the graph described in Kahane's paper,
# to every key correspond one or two values, representing the linked indices.
# All values in `links` are integers, negative numbers are used in the case
# where it is necessary to insert gamma matrices between free indices, in
# order to make Kahane's algorithm work (see paper).
links = {i: [] for i in range(first_dum_pos, total_number)}
# `cum_sign` is a step variable to mark the sign of every index, see paper.
cum_sign = -1
# `cum_sign_list` keeps storage for all `cum_sign` (every index).
cum_sign_list = [None]*total_number
block_free_count = 0
# multiply `resulting_coeff` by the coefficient parameter, the rest
# of the algorithm ignores a scalar coefficient.
resulting_coeff = S.One
# initialize a list of lists of indices. The outer list will contain all
# additive tensor expressions, while the inner list will contain the
# free indices (rearranged according to the algorithm).
resulting_indices = [[]]
# start to count the `connected_components`, which together with the number
# of contractions, determines a -1 or +1 factor to be multiplied.
connected_components = 1
# First loop: here we fill `cum_sign_list`, and draw the links
# among consecutive indices (they are stored in `links`). Links among
# non-consecutive indices will be drawn later.
for i, is_free in enumerate(index_is_free):
# if `expression` starts with free indices, they are ignored here;
# they are later added as they are to the beginning of all
# `resulting_indices` list of lists of indices.
if i < first_dum_pos:
continue
if is_free:
block_free_count += 1
# if previous index was free as well, draw an arch in `links`.
if block_free_count > 1:
links[i - 1].append(i)
links[i].append(i - 1)
else:
# Change the sign of the index (`cum_sign`) if the number of free
# indices preceding it is even.
cum_sign *= 1 if (block_free_count % 2) else -1
if block_free_count == 0 and i != first_dum_pos:
# check if there are two consecutive dummy indices:
# in this case create virtual indices with negative position,
# these "virtual" indices represent the insertion of two
# gamma^0 matrices to separate consecutive dummy indices, as
# Kahane's algorithm requires dummy indices to be separated by
# free indices. The product of two gamma^0 matrices is unity,
# so the new expression being examined is the same as the
# original one.
if cum_sign == -1:
links[-1-i] = [-1-i+1]
links[-1-i+1] = [-1-i]
if (i - cum_sign) in links:
if i != first_dum_pos:
links[i].append(i - cum_sign)
if block_free_count != 0:
if i - cum_sign < len(index_is_free):
if index_is_free[i - cum_sign]:
links[i - cum_sign].append(i)
block_free_count = 0
cum_sign_list[i] = cum_sign
# The previous loop has only created links between consecutive free indices,
# it is necessary to properly create links among dummy (contracted) indices,
# according to the rules described in Kahane's paper. There is only one exception
# to Kahane's rules: the negative indices, which handle the case of some
# consecutive free indices (Kahane's paper just describes dummy indices
# separated by free indices, hinting that free indices can be added without
# altering the expression result).
for i in dum:
# get the positions of the two contracted indices:
pos1 = i[0]
pos2 = i[1]
# create Kahane's upper links, i.e. the upper arcs between dummy
# (i.e. contracted) indices:
links[pos1].append(pos2)
links[pos2].append(pos1)
# create Kahane's lower links, this corresponds to the arcs below
# the line described in the paper:
# first we move `pos1` and `pos2` according to the sign of the indices:
linkpos1 = pos1 + cum_sign_list[pos1]
linkpos2 = pos2 + cum_sign_list[pos2]
# otherwise, perform some checks before creating the lower arcs:
# make sure we are not exceeding the total number of indices:
if linkpos1 >= total_number:
continue
if linkpos2 >= total_number:
continue
# make sure we are not below the first dummy index in `expression`:
if linkpos1 < first_dum_pos:
continue
if linkpos2 < first_dum_pos:
continue
# check if the previous loop created "virtual" indices between dummy
# indices, in such a case relink `linkpos1` and `linkpos2`:
if (-1-linkpos1) in links:
linkpos1 = -1-linkpos1
if (-1-linkpos2) in links:
linkpos2 = -1-linkpos2
# move only if not next to free index:
if linkpos1 >= 0 and not index_is_free[linkpos1]:
linkpos1 = pos1
if linkpos2 >=0 and not index_is_free[linkpos2]:
linkpos2 = pos2
# create the lower arcs:
if linkpos2 not in links[linkpos1]:
links[linkpos1].append(linkpos2)
if linkpos1 not in links[linkpos2]:
links[linkpos2].append(linkpos1)
# This loop starts from the `first_dum_pos` index (first dummy index)
# walks through the graph deleting the visited indices from `links`,
# it adds a gamma matrix for every free index in encounters, while it
# completely ignores dummy indices and virtual indices.
pointer = first_dum_pos
previous_pointer = 0
while True:
if pointer in links:
next_ones = links.pop(pointer)
else:
break
if previous_pointer in next_ones:
next_ones.remove(previous_pointer)
previous_pointer = pointer
if next_ones:
pointer = next_ones[0]
else:
break
if pointer == previous_pointer:
break
if pointer >=0 and free_pos[pointer] is not None:
for ri in resulting_indices:
ri.append(free_pos[pointer])
# The following loop removes the remaining connected components in `links`.
# If there are free indices inside a connected component, it gives a
# contribution to the resulting expression given by the factor
# `gamma_a gamma_b ... gamma_z + gamma_z ... gamma_b gamma_a`, in Kahanes's
# paper represented as {gamma_a, gamma_b, ... , gamma_z},
# virtual indices are ignored. The variable `connected_components` is
# increased by one for every connected component this loop encounters.
# If the connected component has virtual and dummy indices only
# (no free indices), it contributes to `resulting_indices` by a factor of two.
# The multiplication by two is a result of the
# factor {gamma^0, gamma^0} = 2 I, as it appears in Kahane's paper.
# Note: curly brackets are meant as in the paper, as a generalized
# multi-element anticommutator!
while links:
connected_components += 1
pointer = min(links.keys())
previous_pointer = pointer
# the inner loop erases the visited indices from `links`, and it adds
# all free indices to `prepend_indices` list, virtual indices are
# ignored.
prepend_indices = []
while True:
if pointer in links:
next_ones = links.pop(pointer)
else:
break
if previous_pointer in next_ones:
if len(next_ones) > 1:
next_ones.remove(previous_pointer)
previous_pointer = pointer
if next_ones:
pointer = next_ones[0]
if pointer >= first_dum_pos and free_pos[pointer] is not None:
prepend_indices.insert(0, free_pos[pointer])
# if `prepend_indices` is void, it means there are no free indices
# in the loop (and it can be shown that there must be a virtual index),
# loops of virtual indices only contribute by a factor of two:
if len(prepend_indices) == 0:
resulting_coeff *= 2
# otherwise, add the free indices in `prepend_indices` to
# the `resulting_indices`:
else:
expr1 = prepend_indices
expr2 = list(reversed(prepend_indices))
resulting_indices = [expri + ri for ri in resulting_indices for expri in (expr1, expr2)]
# sign correction, as described in Kahane's paper:
resulting_coeff *= -1 if (number_of_contractions - connected_components + 1) % 2 else 1
# power of two factor, as described in Kahane's paper:
resulting_coeff *= 2**(number_of_contractions)
# If `first_dum_pos` is not zero, it means that there are trailing free gamma
# matrices in front of `expression`, so multiply by them:
resulting_indices = [ free_pos[0:first_dum_pos] + ri for ri in resulting_indices ]
resulting_expr = S.Zero
for i in resulting_indices:
temp_expr = S.One
for j in i:
temp_expr *= GammaMatrix(j)
resulting_expr += temp_expr
t = resulting_coeff * resulting_expr
t1 = None
if isinstance(t, TensAdd):
t1 = t.args[0]
elif isinstance(t, TensMul):
t1 = t
if t1:
pass
else:
t = eye(4)*t
return t
@@ -0,0 +1,427 @@
from sympy.matrices.dense import eye, Matrix
from sympy.tensor.tensor import tensor_indices, TensorHead, tensor_heads, \
TensExpr, canon_bp
from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex, \
kahane_simplify, gamma_trace, _simplify_single_line, simplify_gamma_expression
from sympy import Symbol
def _is_tensor_eq(arg1, arg2):
arg1 = canon_bp(arg1)
arg2 = canon_bp(arg2)
if isinstance(arg1, TensExpr):
return arg1.equals(arg2)
elif isinstance(arg2, TensExpr):
return arg2.equals(arg1)
return arg1 == arg2
def execute_gamma_simplify_tests_for_function(tfunc, D):
"""
Perform tests to check if sfunc is able to simplify gamma matrix expressions.
Parameters
==========
`sfunc` a function to simplify a `TIDS`, shall return the simplified `TIDS`.
`D` the number of dimension (in most cases `D=4`).
"""
mu, nu, rho, sigma = tensor_indices("mu, nu, rho, sigma", LorentzIndex)
a1, a2, a3, a4, a5, a6 = tensor_indices("a1:7", LorentzIndex)
mu11, mu12, mu21, mu31, mu32, mu41, mu51, mu52 = tensor_indices("mu11, mu12, mu21, mu31, mu32, mu41, mu51, mu52", LorentzIndex)
mu61, mu71, mu72 = tensor_indices("mu61, mu71, mu72", LorentzIndex)
m0, m1, m2, m3, m4, m5, m6 = tensor_indices("m0:7", LorentzIndex)
def g(xx, yy):
return (G(xx)*G(yy) + G(yy)*G(xx))/2
# Some examples taken from Kahane's paper, 4 dim only:
if D == 4:
t = (G(a1)*G(mu11)*G(a2)*G(mu21)*G(-a1)*G(mu31)*G(-a2))
assert _is_tensor_eq(tfunc(t), -4*G(mu11)*G(mu31)*G(mu21) - 4*G(mu31)*G(mu11)*G(mu21))
t = (G(a1)*G(mu11)*G(mu12)*\
G(a2)*G(mu21)*\
G(a3)*G(mu31)*G(mu32)*\
G(a4)*G(mu41)*\
G(-a2)*G(mu51)*G(mu52)*\
G(-a1)*G(mu61)*\
G(-a3)*G(mu71)*G(mu72)*\
G(-a4))
assert _is_tensor_eq(tfunc(t), \
16*G(mu31)*G(mu32)*G(mu72)*G(mu71)*G(mu11)*G(mu52)*G(mu51)*G(mu12)*G(mu61)*G(mu21)*G(mu41) + 16*G(mu31)*G(mu32)*G(mu72)*G(mu71)*G(mu12)*G(mu51)*G(mu52)*G(mu11)*G(mu61)*G(mu21)*G(mu41) + 16*G(mu71)*G(mu72)*G(mu32)*G(mu31)*G(mu11)*G(mu52)*G(mu51)*G(mu12)*G(mu61)*G(mu21)*G(mu41) + 16*G(mu71)*G(mu72)*G(mu32)*G(mu31)*G(mu12)*G(mu51)*G(mu52)*G(mu11)*G(mu61)*G(mu21)*G(mu41))
# Fully Lorentz-contracted expressions, these return scalars:
def add_delta(ne):
return ne * eye(4) # DiracSpinorIndex.delta(DiracSpinorIndex.auto_left, -DiracSpinorIndex.auto_right)
t = (G(mu)*G(-mu))
ts = add_delta(D)
assert _is_tensor_eq(tfunc(t), ts)
t = (G(mu)*G(nu)*G(-mu)*G(-nu))
ts = add_delta(2*D - D**2) # -8
assert _is_tensor_eq(tfunc(t), ts)
t = (G(mu)*G(nu)*G(-nu)*G(-mu))
ts = add_delta(D**2) # 16
assert _is_tensor_eq(tfunc(t), ts)
t = (G(mu)*G(nu)*G(-rho)*G(-nu)*G(-mu)*G(rho))
ts = add_delta(4*D - 4*D**2 + D**3) # 16
assert _is_tensor_eq(tfunc(t), ts)
t = (G(mu)*G(nu)*G(rho)*G(-rho)*G(-nu)*G(-mu))
ts = add_delta(D**3) # 64
assert _is_tensor_eq(tfunc(t), ts)
t = (G(a1)*G(a2)*G(a3)*G(a4)*G(-a3)*G(-a1)*G(-a2)*G(-a4))
ts = add_delta(-8*D + 16*D**2 - 8*D**3 + D**4) # -32
assert _is_tensor_eq(tfunc(t), ts)
t = (G(-mu)*G(-nu)*G(-rho)*G(-sigma)*G(nu)*G(mu)*G(sigma)*G(rho))
ts = add_delta(-16*D + 24*D**2 - 8*D**3 + D**4) # 64
assert _is_tensor_eq(tfunc(t), ts)
t = (G(-mu)*G(nu)*G(-rho)*G(sigma)*G(rho)*G(-nu)*G(mu)*G(-sigma))
ts = add_delta(8*D - 12*D**2 + 6*D**3 - D**4) # -32
assert _is_tensor_eq(tfunc(t), ts)
t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(-a3)*G(-a2)*G(-a1)*G(-a5)*G(-a4))
ts = add_delta(64*D - 112*D**2 + 60*D**3 - 12*D**4 + D**5) # 256
assert _is_tensor_eq(tfunc(t), ts)
t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(-a3)*G(-a1)*G(-a2)*G(-a4)*G(-a5))
ts = add_delta(64*D - 120*D**2 + 72*D**3 - 16*D**4 + D**5) # -128
assert _is_tensor_eq(tfunc(t), ts)
t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(a6)*G(-a3)*G(-a2)*G(-a1)*G(-a6)*G(-a5)*G(-a4))
ts = add_delta(416*D - 816*D**2 + 528*D**3 - 144*D**4 + 18*D**5 - D**6) # -128
assert _is_tensor_eq(tfunc(t), ts)
t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(a6)*G(-a2)*G(-a3)*G(-a1)*G(-a6)*G(-a4)*G(-a5))
ts = add_delta(416*D - 848*D**2 + 584*D**3 - 172*D**4 + 22*D**5 - D**6) # -128
assert _is_tensor_eq(tfunc(t), ts)
# Expressions with free indices:
t = (G(mu)*G(nu)*G(rho)*G(sigma)*G(-mu))
assert _is_tensor_eq(tfunc(t), (-2*G(sigma)*G(rho)*G(nu) + (4-D)*G(nu)*G(rho)*G(sigma)))
t = (G(mu)*G(nu)*G(-mu))
assert _is_tensor_eq(tfunc(t), (2-D)*G(nu))
t = (G(mu)*G(nu)*G(rho)*G(-mu))
assert _is_tensor_eq(tfunc(t), 2*G(nu)*G(rho) + 2*G(rho)*G(nu) - (4-D)*G(nu)*G(rho))
t = 2*G(m2)*G(m0)*G(m1)*G(-m0)*G(-m1)
st = tfunc(t)
assert _is_tensor_eq(st, (D*(-2*D + 4))*G(m2))
t = G(m2)*G(m0)*G(m1)*G(-m0)*G(-m2)
st = tfunc(t)
assert _is_tensor_eq(st, ((-D + 2)**2)*G(m1))
t = G(m0)*G(m1)*G(m2)*G(m3)*G(-m1)
st = tfunc(t)
assert _is_tensor_eq(st, (D - 4)*G(m0)*G(m2)*G(m3) + 4*G(m0)*g(m2, m3))
t = G(m0)*G(m1)*G(m2)*G(m3)*G(-m1)*G(-m0)
st = tfunc(t)
assert _is_tensor_eq(st, ((D - 4)**2)*G(m2)*G(m3) + (8*D - 16)*g(m2, m3))
t = G(m2)*G(m0)*G(m1)*G(-m2)*G(-m0)
st = tfunc(t)
assert _is_tensor_eq(st, ((-D + 2)*(D - 4) + 4)*G(m1))
t = G(m3)*G(m1)*G(m0)*G(m2)*G(-m3)*G(-m0)*G(-m2)
st = tfunc(t)
assert _is_tensor_eq(st, (-4*D + (-D + 2)**2*(D - 4) + 8)*G(m1))
t = 2*G(m0)*G(m1)*G(m2)*G(m3)*G(-m0)
st = tfunc(t)
assert _is_tensor_eq(st, ((-2*D + 8)*G(m1)*G(m2)*G(m3) - 4*G(m3)*G(m2)*G(m1)))
t = G(m5)*G(m0)*G(m1)*G(m4)*G(m2)*G(-m4)*G(m3)*G(-m0)
st = tfunc(t)
assert _is_tensor_eq(st, (((-D + 2)*(-D + 4))*G(m5)*G(m1)*G(m2)*G(m3) + (2*D - 4)*G(m5)*G(m3)*G(m2)*G(m1)))
t = -G(m0)*G(m1)*G(m2)*G(m3)*G(-m0)*G(m4)
st = tfunc(t)
assert _is_tensor_eq(st, ((D - 4)*G(m1)*G(m2)*G(m3)*G(m4) + 2*G(m3)*G(m2)*G(m1)*G(m4)))
t = G(-m5)*G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)*G(m5)
st = tfunc(t)
result1 = ((-D + 4)**2 + 4)*G(m1)*G(m2)*G(m3)*G(m4) +\
(4*D - 16)*G(m3)*G(m2)*G(m1)*G(m4) + (4*D - 16)*G(m4)*G(m1)*G(m2)*G(m3)\
+ 4*G(m2)*G(m1)*G(m4)*G(m3) + 4*G(m3)*G(m4)*G(m1)*G(m2) +\
4*G(m4)*G(m3)*G(m2)*G(m1)
# Kahane's algorithm yields this result, which is equivalent to `result1`
# in four dimensions, but is not automatically recognized as equal:
result2 = 8*G(m1)*G(m2)*G(m3)*G(m4) + 8*G(m4)*G(m3)*G(m2)*G(m1)
if D == 4:
assert _is_tensor_eq(st, (result1)) or _is_tensor_eq(st, (result2))
else:
assert _is_tensor_eq(st, (result1))
# and a few very simple cases, with no contracted indices:
t = G(m0)
st = tfunc(t)
assert _is_tensor_eq(st, t)
t = -7*G(m0)
st = tfunc(t)
assert _is_tensor_eq(st, t)
t = 224*G(m0)*G(m1)*G(-m2)*G(m3)
st = tfunc(t)
assert _is_tensor_eq(st, t)
def test_kahane_algorithm():
# Wrap this function to convert to and from TIDS:
def tfunc(e):
return _simplify_single_line(e)
execute_gamma_simplify_tests_for_function(tfunc, D=4)
def test_kahane_simplify1():
i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15 = tensor_indices('i0:16', LorentzIndex)
mu, nu, rho, sigma = tensor_indices("mu, nu, rho, sigma", LorentzIndex)
D = 4
t = G(i0)*G(i1)
r = kahane_simplify(t)
assert r.equals(t)
t = G(i0)*G(i1)*G(-i0)
r = kahane_simplify(t)
assert r.equals(-2*G(i1))
t = G(i0)*G(i1)*G(-i0)
r = kahane_simplify(t)
assert r.equals(-2*G(i1))
t = G(i0)*G(i1)
r = kahane_simplify(t)
assert r.equals(t)
t = G(i0)*G(i1)
r = kahane_simplify(t)
assert r.equals(t)
t = G(i0)*G(-i0)
r = kahane_simplify(t)
assert r.equals(4*eye(4))
t = G(i0)*G(-i0)
r = kahane_simplify(t)
assert r.equals(4*eye(4))
t = G(i0)*G(-i0)
r = kahane_simplify(t)
assert r.equals(4*eye(4))
t = G(i0)*G(i1)*G(-i0)
r = kahane_simplify(t)
assert r.equals(-2*G(i1))
t = G(i0)*G(i1)*G(-i0)*G(-i1)
r = kahane_simplify(t)
assert r.equals((2*D - D**2)*eye(4))
t = G(i0)*G(i1)*G(-i0)*G(-i1)
r = kahane_simplify(t)
assert r.equals((2*D - D**2)*eye(4))
t = G(i0)*G(-i0)*G(i1)*G(-i1)
r = kahane_simplify(t)
assert r.equals(16*eye(4))
t = (G(mu)*G(nu)*G(-nu)*G(-mu))
r = kahane_simplify(t)
assert r.equals(D**2*eye(4))
t = (G(mu)*G(nu)*G(-nu)*G(-mu))
r = kahane_simplify(t)
assert r.equals(D**2*eye(4))
t = (G(mu)*G(nu)*G(-nu)*G(-mu))
r = kahane_simplify(t)
assert r.equals(D**2*eye(4))
t = (G(mu)*G(nu)*G(-rho)*G(-nu)*G(-mu)*G(rho))
r = kahane_simplify(t)
assert r.equals((4*D - 4*D**2 + D**3)*eye(4))
t = (G(-mu)*G(-nu)*G(-rho)*G(-sigma)*G(nu)*G(mu)*G(sigma)*G(rho))
r = kahane_simplify(t)
assert r.equals((-16*D + 24*D**2 - 8*D**3 + D**4)*eye(4))
t = (G(-mu)*G(nu)*G(-rho)*G(sigma)*G(rho)*G(-nu)*G(mu)*G(-sigma))
r = kahane_simplify(t)
assert r.equals((8*D - 12*D**2 + 6*D**3 - D**4)*eye(4))
# Expressions with free indices:
t = (G(mu)*G(nu)*G(rho)*G(sigma)*G(-mu))
r = kahane_simplify(t)
assert r.equals(-2*G(sigma)*G(rho)*G(nu))
t = (G(mu)*G(-mu)*G(rho)*G(sigma))
r = kahane_simplify(t)
assert r.equals(4*G(rho)*G(sigma))
t = (G(rho)*G(sigma)*G(mu)*G(-mu))
r = kahane_simplify(t)
assert r.equals(4*G(rho)*G(sigma))
def test_gamma_matrix_class():
i, j, k = tensor_indices('i,j,k', LorentzIndex)
# define another type of TensorHead to see if exprs are correctly handled:
A = TensorHead('A', [LorentzIndex])
t = A(k)*G(i)*G(-i)
ts = simplify_gamma_expression(t)
assert _is_tensor_eq(ts, Matrix([
[4, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 4, 0],
[0, 0, 0, 4]])*A(k))
t = G(i)*A(k)*G(j)
ts = simplify_gamma_expression(t)
assert _is_tensor_eq(ts, A(k)*G(i)*G(j))
execute_gamma_simplify_tests_for_function(simplify_gamma_expression, D=4)
def test_gamma_matrix_trace():
g = LorentzIndex.metric
m0, m1, m2, m3, m4, m5, m6 = tensor_indices('m0:7', LorentzIndex)
n0, n1, n2, n3, n4, n5 = tensor_indices('n0:6', LorentzIndex)
# working in D=4 dimensions
D = 4
# traces of odd number of gamma matrices are zero:
t = G(m0)
t1 = gamma_trace(t)
assert t1.equals(0)
t = G(m0)*G(m1)*G(m2)
t1 = gamma_trace(t)
assert t1.equals(0)
t = G(m0)*G(m1)*G(-m0)
t1 = gamma_trace(t)
assert t1.equals(0)
t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)
t1 = gamma_trace(t)
assert t1.equals(0)
# traces without internal contractions:
t = G(m0)*G(m1)
t1 = gamma_trace(t)
assert _is_tensor_eq(t1, 4*g(m0, m1))
t = G(m0)*G(m1)*G(m2)*G(m3)
t1 = gamma_trace(t)
t2 = -4*g(m0, m2)*g(m1, m3) + 4*g(m0, m1)*g(m2, m3) + 4*g(m0, m3)*g(m1, m2)
assert _is_tensor_eq(t1, t2)
t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(m5)
t1 = gamma_trace(t)
t2 = t1*g(-m0, -m5)
t2 = t2.contract_metric(g)
assert _is_tensor_eq(t2, D*gamma_trace(G(m1)*G(m2)*G(m3)*G(m4)))
# traces of expressions with internal contractions:
t = G(m0)*G(-m0)
t1 = gamma_trace(t)
assert t1.equals(4*D)
t = G(m0)*G(m1)*G(-m0)*G(-m1)
t1 = gamma_trace(t)
assert t1.equals(8*D - 4*D**2)
t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)
t1 = gamma_trace(t)
t2 = (-4*D)*g(m1, m3)*g(m2, m4) + (4*D)*g(m1, m2)*g(m3, m4) + \
(4*D)*g(m1, m4)*g(m2, m3)
assert _is_tensor_eq(t1, t2)
t = G(-m5)*G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)*G(m5)
t1 = gamma_trace(t)
t2 = (32*D + 4*(-D + 4)**2 - 64)*(g(m1, m2)*g(m3, m4) - \
g(m1, m3)*g(m2, m4) + g(m1, m4)*g(m2, m3))
assert _is_tensor_eq(t1, t2)
t = G(m0)*G(m1)*G(-m0)*G(m3)
t1 = gamma_trace(t)
assert t1.equals((-4*D + 8)*g(m1, m3))
# p, q = S1('p,q')
# ps = p(m0)*G(-m0)
# qs = q(m0)*G(-m0)
# t = ps*qs*ps*qs
# t1 = gamma_trace(t)
# assert t1 == 8*p(m0)*q(-m0)*p(m1)*q(-m1) - 4*p(m0)*p(-m0)*q(m1)*q(-m1)
t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(m5)*G(-m0)*G(-m1)*G(-m2)*G(-m3)*G(-m4)*G(-m5)
t1 = gamma_trace(t)
assert t1.equals(-4*D**6 + 120*D**5 - 1040*D**4 + 3360*D**3 - 4480*D**2 + 2048*D)
t = G(m0)*G(m1)*G(n1)*G(m2)*G(n2)*G(m3)*G(m4)*G(-n2)*G(-n1)*G(-m0)*G(-m1)*G(-m2)*G(-m3)*G(-m4)
t1 = gamma_trace(t)
tresu = -7168*D + 16768*D**2 - 14400*D**3 + 5920*D**4 - 1232*D**5 + 120*D**6 - 4*D**7
assert t1.equals(tresu)
# checked with Mathematica
# In[1]:= <<Tracer.m
# In[2]:= Spur[l];
# In[3]:= GammaTrace[l, {m0},{m1},{n1},{m2},{n2},{m3},{m4},{n3},{n4},{m0},{m1},{m2},{m3},{m4}]
t = G(m0)*G(m1)*G(n1)*G(m2)*G(n2)*G(m3)*G(m4)*G(n3)*G(n4)*G(-m0)*G(-m1)*G(-m2)*G(-m3)*G(-m4)
t1 = gamma_trace(t)
# t1 = t1.expand_coeff()
c1 = -4*D**5 + 120*D**4 - 1200*D**3 + 5280*D**2 - 10560*D + 7808
c2 = -4*D**5 + 88*D**4 - 560*D**3 + 1440*D**2 - 1600*D + 640
assert _is_tensor_eq(t1, c1*g(n1, n4)*g(n2, n3) + c2*g(n1, n2)*g(n3, n4) + \
(-c1)*g(n1, n3)*g(n2, n4))
p, q = tensor_heads('p,q', [LorentzIndex])
ps = p(m0)*G(-m0)
qs = q(m0)*G(-m0)
p2 = p(m0)*p(-m0)
q2 = q(m0)*q(-m0)
pq = p(m0)*q(-m0)
t = ps*qs*ps*qs
r = gamma_trace(t)
assert _is_tensor_eq(r, 8*pq*pq - 4*p2*q2)
t = ps*qs*ps*qs*ps*qs
r = gamma_trace(t)
assert _is_tensor_eq(r, -12*p2*pq*q2 + 16*pq*pq*pq)
t = ps*qs*ps*qs*ps*qs*ps*qs
r = gamma_trace(t)
assert _is_tensor_eq(r, -32*pq*pq*p2*q2 + 32*pq*pq*pq*pq + 4*p2*p2*q2*q2)
t = 4*p(m1)*p(m0)*p(-m0)*q(-m1)*q(m2)*q(-m2)
assert _is_tensor_eq(gamma_trace(t), t)
t = ps*ps*ps*ps*ps*ps*ps*ps
r = gamma_trace(t)
assert r.equals(4*p2*p2*p2*p2)
def test_bug_13636():
"""Test issue 13636 regarding handling traces of sums of products
of GammaMatrix mixed with other factors."""
pi, ki, pf = tensor_heads("pi, ki, pf", [LorentzIndex])
i0, i1, i2, i3, i4 = tensor_indices("i0:5", LorentzIndex)
x = Symbol("x")
pis = pi(i2) * G(-i2)
kis = ki(i3) * G(-i3)
pfs = pf(i4) * G(-i4)
a = pfs * G(i0) * kis * G(i1) * pis * G(-i1) * kis * G(-i0)
b = pfs * G(i0) * kis * G(i1) * pis * x * G(-i0) * pi(-i1)
ta = gamma_trace(a)
tb = gamma_trace(b)
t_a_plus_b = gamma_trace(a + b)
assert ta == 4 * (
-4 * ki(i0) * ki(-i0) * pf(i1) * pi(-i1)
+ 8 * ki(i0) * ki(i1) * pf(-i0) * pi(-i1)
)
assert tb == -8 * x * ki(i0) * pf(-i0) * pi(i1) * pi(-i1)
assert t_a_plus_b == ta + tb
@@ -0,0 +1,265 @@
from sympy.core.numbers import Float
from sympy.core.singleton import S
from sympy.functions.combinatorial.factorials import factorial
from sympy.functions.elementary.exponential import exp
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.special.polynomials import assoc_laguerre
from sympy.functions.special.spherical_harmonics import Ynm
def R_nl(n, l, r, Z=1):
"""
Returns the Hydrogen radial wavefunction R_{nl}.
Parameters
==========
n : integer
Principal Quantum Number which is
an integer with possible values as 1, 2, 3, 4,...
l : integer
``l`` is the Angular Momentum Quantum Number with
values ranging from 0 to ``n-1``.
r :
Radial coordinate.
Z :
Atomic number (1 for Hydrogen, 2 for Helium, ...)
Everything is in Hartree atomic units.
Examples
========
>>> from sympy.physics.hydrogen import R_nl
>>> from sympy.abc import r, Z
>>> R_nl(1, 0, r, Z)
2*sqrt(Z**3)*exp(-Z*r)
>>> R_nl(2, 0, r, Z)
sqrt(2)*(-Z*r + 2)*sqrt(Z**3)*exp(-Z*r/2)/4
>>> R_nl(2, 1, r, Z)
sqrt(6)*Z*r*sqrt(Z**3)*exp(-Z*r/2)/12
For Hydrogen atom, you can just use the default value of Z=1:
>>> R_nl(1, 0, r)
2*exp(-r)
>>> R_nl(2, 0, r)
sqrt(2)*(2 - r)*exp(-r/2)/4
>>> R_nl(3, 0, r)
2*sqrt(3)*(2*r**2/9 - 2*r + 3)*exp(-r/3)/27
For Silver atom, you would use Z=47:
>>> R_nl(1, 0, r, Z=47)
94*sqrt(47)*exp(-47*r)
>>> R_nl(2, 0, r, Z=47)
47*sqrt(94)*(2 - 47*r)*exp(-47*r/2)/4
>>> R_nl(3, 0, r, Z=47)
94*sqrt(141)*(4418*r**2/9 - 94*r + 3)*exp(-47*r/3)/27
The normalization of the radial wavefunction is:
>>> from sympy import integrate, oo
>>> integrate(R_nl(1, 0, r)**2 * r**2, (r, 0, oo))
1
>>> integrate(R_nl(2, 0, r)**2 * r**2, (r, 0, oo))
1
>>> integrate(R_nl(2, 1, r)**2 * r**2, (r, 0, oo))
1
It holds for any atomic number:
>>> integrate(R_nl(1, 0, r, Z=2)**2 * r**2, (r, 0, oo))
1
>>> integrate(R_nl(2, 0, r, Z=3)**2 * r**2, (r, 0, oo))
1
>>> integrate(R_nl(2, 1, r, Z=4)**2 * r**2, (r, 0, oo))
1
"""
# sympify arguments
n, l, r, Z = map(S, [n, l, r, Z])
# radial quantum number
n_r = n - l - 1
# rescaled "r"
a = 1/Z # Bohr radius
r0 = 2 * r / (n * a)
# normalization coefficient
C = sqrt((S(2)/(n*a))**3 * factorial(n_r) / (2*n*factorial(n + l)))
# This is an equivalent normalization coefficient, that can be found in
# some books. Both coefficients seem to be the same fast:
# C = S(2)/n**2 * sqrt(1/a**3 * factorial(n_r) / (factorial(n+l)))
return C * r0**l * assoc_laguerre(n_r, 2*l + 1, r0).expand() * exp(-r0/2)
def Psi_nlm(n, l, m, r, phi, theta, Z=1):
"""
Returns the Hydrogen wave function psi_{nlm}. It's the product of
the radial wavefunction R_{nl} and the spherical harmonic Y_{l}^{m}.
Parameters
==========
n : integer
Principal Quantum Number which is
an integer with possible values as 1, 2, 3, 4,...
l : integer
``l`` is the Angular Momentum Quantum Number with
values ranging from 0 to ``n-1``.
m : integer
``m`` is the Magnetic Quantum Number with values
ranging from ``-l`` to ``l``.
r :
radial coordinate
phi :
azimuthal angle
theta :
polar angle
Z :
atomic number (1 for Hydrogen, 2 for Helium, ...)
Everything is in Hartree atomic units.
Examples
========
>>> from sympy.physics.hydrogen import Psi_nlm
>>> from sympy import Symbol
>>> r=Symbol("r", positive=True)
>>> phi=Symbol("phi", real=True)
>>> theta=Symbol("theta", real=True)
>>> Z=Symbol("Z", positive=True, integer=True, nonzero=True)
>>> Psi_nlm(1,0,0,r,phi,theta,Z)
Z**(3/2)*exp(-Z*r)/sqrt(pi)
>>> Psi_nlm(2,1,1,r,phi,theta,Z)
-Z**(5/2)*r*exp(I*phi)*exp(-Z*r/2)*sin(theta)/(8*sqrt(pi))
Integrating the absolute square of a hydrogen wavefunction psi_{nlm}
over the whole space leads 1.
The normalization of the hydrogen wavefunctions Psi_nlm is:
>>> from sympy import integrate, conjugate, pi, oo, sin
>>> wf=Psi_nlm(2,1,1,r,phi,theta,Z)
>>> abs_sqrd=wf*conjugate(wf)
>>> jacobi=r**2*sin(theta)
>>> integrate(abs_sqrd*jacobi, (r,0,oo), (phi,0,2*pi), (theta,0,pi))
1
"""
# sympify arguments
n, l, m, r, phi, theta, Z = map(S, [n, l, m, r, phi, theta, Z])
# check if values for n,l,m make physically sense
if n.is_integer and n < 1:
raise ValueError("'n' must be positive integer")
if l.is_integer and not (n > l):
raise ValueError("'n' must be greater than 'l'")
if m.is_integer and not (abs(m) <= l):
raise ValueError("|'m'| must be less or equal 'l'")
# return the hydrogen wave function
return R_nl(n, l, r, Z)*Ynm(l, m, theta, phi).expand(func=True)
def E_nl(n, Z=1):
"""
Returns the energy of the state (n, l) in Hartree atomic units.
The energy does not depend on "l".
Parameters
==========
n : integer
Principal Quantum Number which is
an integer with possible values as 1, 2, 3, 4,...
Z :
Atomic number (1 for Hydrogen, 2 for Helium, ...)
Examples
========
>>> from sympy.physics.hydrogen import E_nl
>>> from sympy.abc import n, Z
>>> E_nl(n, Z)
-Z**2/(2*n**2)
>>> E_nl(1)
-1/2
>>> E_nl(2)
-1/8
>>> E_nl(3)
-1/18
>>> E_nl(3, 47)
-2209/18
"""
n, Z = S(n), S(Z)
if n.is_integer and (n < 1):
raise ValueError("'n' must be positive integer")
return -Z**2/(2*n**2)
def E_nl_dirac(n, l, spin_up=True, Z=1, c=Float("137.035999037")):
"""
Returns the relativistic energy of the state (n, l, spin) in Hartree atomic
units.
The energy is calculated from the Dirac equation. The rest mass energy is
*not* included.
Parameters
==========
n : integer
Principal Quantum Number which is
an integer with possible values as 1, 2, 3, 4,...
l : integer
``l`` is the Angular Momentum Quantum Number with
values ranging from 0 to ``n-1``.
spin_up :
True if the electron spin is up (default), otherwise down
Z :
Atomic number (1 for Hydrogen, 2 for Helium, ...)
c :
Speed of light in atomic units. Default value is 137.035999037,
taken from https://arxiv.org/abs/1012.3627
Examples
========
>>> from sympy.physics.hydrogen import E_nl_dirac
>>> E_nl_dirac(1, 0)
-0.500006656595360
>>> E_nl_dirac(2, 0)
-0.125002080189006
>>> E_nl_dirac(2, 1)
-0.125000416028342
>>> E_nl_dirac(2, 1, False)
-0.125002080189006
>>> E_nl_dirac(3, 0)
-0.0555562951740285
>>> E_nl_dirac(3, 1)
-0.0555558020932949
>>> E_nl_dirac(3, 1, False)
-0.0555562951740285
>>> E_nl_dirac(3, 2)
-0.0555556377366884
>>> E_nl_dirac(3, 2, False)
-0.0555558020932949
"""
n, l, Z, c = map(S, [n, l, Z, c])
if not (l >= 0):
raise ValueError("'l' must be positive or zero")
if not (n > l):
raise ValueError("'n' must be greater than 'l'")
if (l == 0 and spin_up is False):
raise ValueError("Spin must be up for l==0.")
# skappa is sign*kappa, where sign contains the correct sign
if spin_up:
skappa = -l - 1
else:
skappa = -l
beta = sqrt(skappa**2 - Z**2/c**2)
return c**2/sqrt(1 + Z**2/(n + skappa + beta)**2/c**2) - c**2
@@ -0,0 +1,176 @@
"""Known matrices related to physics"""
from sympy.core.numbers import I
from sympy.matrices.dense import MutableDenseMatrix as Matrix
from sympy.utilities.decorator import deprecated
def msigma(i):
r"""Returns a Pauli matrix `\sigma_i` with `i=1,2,3`.
References
==========
.. [1] https://en.wikipedia.org/wiki/Pauli_matrices
Examples
========
>>> from sympy.physics.matrices import msigma
>>> msigma(1)
Matrix([
[0, 1],
[1, 0]])
"""
if i == 1:
mat = (
(0, 1),
(1, 0)
)
elif i == 2:
mat = (
(0, -I),
(I, 0)
)
elif i == 3:
mat = (
(1, 0),
(0, -1)
)
else:
raise IndexError("Invalid Pauli index")
return Matrix(mat)
def pat_matrix(m, dx, dy, dz):
"""Returns the Parallel Axis Theorem matrix to translate the inertia
matrix a distance of `(dx, dy, dz)` for a body of mass m.
Examples
========
To translate a body having a mass of 2 units a distance of 1 unit along
the `x`-axis we get:
>>> from sympy.physics.matrices import pat_matrix
>>> pat_matrix(2, 1, 0, 0)
Matrix([
[0, 0, 0],
[0, 2, 0],
[0, 0, 2]])
"""
dxdy = -dx*dy
dydz = -dy*dz
dzdx = -dz*dx
dxdx = dx**2
dydy = dy**2
dzdz = dz**2
mat = ((dydy + dzdz, dxdy, dzdx),
(dxdy, dxdx + dzdz, dydz),
(dzdx, dydz, dydy + dxdx))
return m*Matrix(mat)
def mgamma(mu, lower=False):
r"""Returns a Dirac gamma matrix `\gamma^\mu` in the standard
(Dirac) representation.
Explanation
===========
If you want `\gamma_\mu`, use ``gamma(mu, True)``.
We use a convention:
`\gamma^5 = i \cdot \gamma^0 \cdot \gamma^1 \cdot \gamma^2 \cdot \gamma^3`
`\gamma_5 = i \cdot \gamma_0 \cdot \gamma_1 \cdot \gamma_2 \cdot \gamma_3 = - \gamma^5`
References
==========
.. [1] https://en.wikipedia.org/wiki/Gamma_matrices
Examples
========
>>> from sympy.physics.matrices import mgamma
>>> mgamma(1)
Matrix([
[ 0, 0, 0, 1],
[ 0, 0, 1, 0],
[ 0, -1, 0, 0],
[-1, 0, 0, 0]])
"""
if mu not in (0, 1, 2, 3, 5):
raise IndexError("Invalid Dirac index")
if mu == 0:
mat = (
(1, 0, 0, 0),
(0, 1, 0, 0),
(0, 0, -1, 0),
(0, 0, 0, -1)
)
elif mu == 1:
mat = (
(0, 0, 0, 1),
(0, 0, 1, 0),
(0, -1, 0, 0),
(-1, 0, 0, 0)
)
elif mu == 2:
mat = (
(0, 0, 0, -I),
(0, 0, I, 0),
(0, I, 0, 0),
(-I, 0, 0, 0)
)
elif mu == 3:
mat = (
(0, 0, 1, 0),
(0, 0, 0, -1),
(-1, 0, 0, 0),
(0, 1, 0, 0)
)
elif mu == 5:
mat = (
(0, 0, 1, 0),
(0, 0, 0, 1),
(1, 0, 0, 0),
(0, 1, 0, 0)
)
m = Matrix(mat)
if lower:
if mu in (1, 2, 3, 5):
m = -m
return m
#Minkowski tensor using the convention (+,-,-,-) used in the Quantum Field
#Theory
minkowski_tensor = Matrix( (
(1, 0, 0, 0),
(0, -1, 0, 0),
(0, 0, -1, 0),
(0, 0, 0, -1)
))
@deprecated(
"""
The sympy.physics.matrices.mdft method is deprecated. Use
sympy.DFT(n).as_explicit() instead.
""",
deprecated_since_version="1.9",
active_deprecations_target="deprecated-physics-mdft",
)
def mdft(n):
r"""
.. deprecated:: 1.9
Use DFT from sympy.matrices.expressions.fourier instead.
To get identical behavior to ``mdft(n)``, use ``DFT(n).as_explicit()``.
"""
from sympy.matrices.expressions.fourier import DFT
return DFT(n).as_mutable()
@@ -0,0 +1,90 @@
__all__ = [
'vector',
'CoordinateSym', 'ReferenceFrame', 'Dyadic', 'Vector', 'Point', 'cross',
'dot', 'express', 'time_derivative', 'outer', 'kinematic_equations',
'get_motion_params', 'partial_velocity', 'dynamicsymbols', 'vprint',
'vsstrrepr', 'vsprint', 'vpprint', 'vlatex', 'init_vprinting', 'curl',
'divergence', 'gradient', 'is_conservative', 'is_solenoidal',
'scalar_potential', 'scalar_potential_difference',
'KanesMethod',
'RigidBody',
'linear_momentum', 'angular_momentum', 'kinetic_energy', 'potential_energy',
'Lagrangian', 'mechanics_printing', 'mprint', 'msprint', 'mpprint',
'mlatex', 'msubs', 'find_dynamicsymbols',
'inertia', 'inertia_of_point_mass', 'Inertia',
'Force', 'Torque',
'Particle',
'LagrangesMethod',
'Linearizer',
'Body',
'SymbolicSystem', 'System',
'PinJoint', 'PrismaticJoint', 'CylindricalJoint', 'PlanarJoint',
'SphericalJoint', 'WeldJoint',
'JointsMethod',
'WrappingCylinder', 'WrappingGeometryBase', 'WrappingSphere',
'PathwayBase', 'LinearPathway', 'ObstacleSetPathway', 'WrappingPathway',
'ActuatorBase', 'ForceActuator', 'LinearDamper', 'LinearSpring',
'TorqueActuator', 'DuffingSpring', 'CoulombKineticFriction',
]
from sympy.physics import vector
from sympy.physics.vector import (CoordinateSym, ReferenceFrame, Dyadic, Vector, Point,
cross, dot, express, time_derivative, outer, kinematic_equations,
get_motion_params, partial_velocity, dynamicsymbols, vprint,
vsstrrepr, vsprint, vpprint, vlatex, init_vprinting, curl, divergence,
gradient, is_conservative, is_solenoidal, scalar_potential,
scalar_potential_difference)
from .kane import KanesMethod
from .rigidbody import RigidBody
from .functions import (linear_momentum, angular_momentum, kinetic_energy,
potential_energy, Lagrangian, mechanics_printing,
mprint, msprint, mpprint, mlatex, msubs,
find_dynamicsymbols)
from .inertia import inertia, inertia_of_point_mass, Inertia
from .loads import Force, Torque
from .particle import Particle
from .lagrange import LagrangesMethod
from .linearize import Linearizer
from .body import Body
from .system import SymbolicSystem, System
from .jointsmethod import JointsMethod
from .joint import (PinJoint, PrismaticJoint, CylindricalJoint, PlanarJoint,
SphericalJoint, WeldJoint)
from .wrapping_geometry import (WrappingCylinder, WrappingGeometryBase,
WrappingSphere)
from .pathway import (PathwayBase, LinearPathway, ObstacleSetPathway,
WrappingPathway)
from .actuator import (ActuatorBase, ForceActuator, LinearDamper, LinearSpring,
TorqueActuator, DuffingSpring, CoulombKineticFriction)
@@ -0,0 +1,710 @@
from sympy import Symbol
from sympy.physics.vector import Point, Vector, ReferenceFrame, Dyadic
from sympy.physics.mechanics import RigidBody, Particle, Inertia
from sympy.physics.mechanics.body_base import BodyBase
from sympy.utilities.exceptions import sympy_deprecation_warning
__all__ = ['Body']
# XXX: We use type:ignore because the classes RigidBody and Particle have
# inconsistent parallel axis methods that take different numbers of arguments.
class Body(RigidBody, Particle): # type: ignore
"""
Body is a common representation of either a RigidBody or a Particle SymPy
object depending on what is passed in during initialization. If a mass is
passed in and central_inertia is left as None, the Particle object is
created. Otherwise a RigidBody object will be created.
.. deprecated:: 1.13
The Body class is deprecated. Its functionality is captured by
:class:`~.RigidBody` and :class:`~.Particle`.
Explanation
===========
The attributes that Body possesses will be the same as a Particle instance
or a Rigid Body instance depending on which was created. Additional
attributes are listed below.
Attributes
==========
name : string
The body's name
masscenter : Point
The point which represents the center of mass of the rigid body
frame : ReferenceFrame
The reference frame which the body is fixed in
mass : Sympifyable
The body's mass
inertia : (Dyadic, Point)
The body's inertia around its center of mass. This attribute is specific
to the rigid body form of Body and is left undefined for the Particle
form
loads : iterable
This list contains information on the different loads acting on the
Body. Forces are listed as a (point, vector) tuple and torques are
listed as (reference frame, vector) tuples.
Parameters
==========
name : String
Defines the name of the body. It is used as the base for defining
body specific properties.
masscenter : Point, optional
A point that represents the center of mass of the body or particle.
If no point is given, a point is generated.
mass : Sympifyable, optional
A Sympifyable object which represents the mass of the body. If no
mass is passed, one is generated.
frame : ReferenceFrame, optional
The ReferenceFrame that represents the reference frame of the body.
If no frame is given, a frame is generated.
central_inertia : Dyadic, optional
Central inertia dyadic of the body. If none is passed while creating
RigidBody, a default inertia is generated.
Examples
========
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
Default behaviour. This results in the creation of a RigidBody object for
which the mass, mass center, frame and inertia attributes are given default
values. ::
>>> from sympy.physics.mechanics import Body
>>> with ignore_warnings(DeprecationWarning):
... body = Body('name_of_body')
This next example demonstrates the code required to specify all of the
values of the Body object. Note this will also create a RigidBody version of
the Body object. ::
>>> from sympy import Symbol
>>> from sympy.physics.mechanics import ReferenceFrame, Point, inertia
>>> from sympy.physics.mechanics import Body
>>> mass = Symbol('mass')
>>> masscenter = Point('masscenter')
>>> frame = ReferenceFrame('frame')
>>> ixx = Symbol('ixx')
>>> body_inertia = inertia(frame, ixx, 0, 0)
>>> with ignore_warnings(DeprecationWarning):
... body = Body('name_of_body', masscenter, mass, frame, body_inertia)
The minimal code required to create a Particle version of the Body object
involves simply passing in a name and a mass. ::
>>> from sympy import Symbol
>>> from sympy.physics.mechanics import Body
>>> mass = Symbol('mass')
>>> with ignore_warnings(DeprecationWarning):
... body = Body('name_of_body', mass=mass)
The Particle version of the Body object can also receive a masscenter point
and a reference frame, just not an inertia.
"""
def __init__(self, name, masscenter=None, mass=None, frame=None,
central_inertia=None):
sympy_deprecation_warning(
"""
Support for the Body class has been removed, as its functionality is
fully captured by RigidBody and Particle.
""",
deprecated_since_version="1.13",
active_deprecations_target="deprecated-mechanics-body-class"
)
self._loads = []
if frame is None:
frame = ReferenceFrame(name + '_frame')
if masscenter is None:
masscenter = Point(name + '_masscenter')
if central_inertia is None and mass is None:
ixx = Symbol(name + '_ixx')
iyy = Symbol(name + '_iyy')
izz = Symbol(name + '_izz')
izx = Symbol(name + '_izx')
ixy = Symbol(name + '_ixy')
iyz = Symbol(name + '_iyz')
_inertia = Inertia.from_inertia_scalars(masscenter, frame, ixx, iyy,
izz, ixy, iyz, izx)
else:
_inertia = (central_inertia, masscenter)
if mass is None:
_mass = Symbol(name + '_mass')
else:
_mass = mass
masscenter.set_vel(frame, 0)
# If user passes masscenter and mass then a particle is created
# otherwise a rigidbody. As a result a body may or may not have inertia.
# Note: BodyBase.__init__ is used to prevent problems with super() calls in
# Particle and RigidBody arising due to multiple inheritance.
if central_inertia is None and mass is not None:
BodyBase.__init__(self, name, masscenter, _mass)
self.frame = frame
self._central_inertia = Dyadic(0)
else:
BodyBase.__init__(self, name, masscenter, _mass)
self.frame = frame
self.inertia = _inertia
def __repr__(self):
if self.is_rigidbody:
return RigidBody.__repr__(self)
return Particle.__repr__(self)
@property
def loads(self):
return self._loads
@property
def x(self):
"""The basis Vector for the Body, in the x direction."""
return self.frame.x
@property
def y(self):
"""The basis Vector for the Body, in the y direction."""
return self.frame.y
@property
def z(self):
"""The basis Vector for the Body, in the z direction."""
return self.frame.z
@property
def inertia(self):
"""The body's inertia about a point; stored as (Dyadic, Point)."""
if self.is_rigidbody:
return RigidBody.inertia.fget(self)
return (self.central_inertia, self.masscenter)
@inertia.setter
def inertia(self, I):
RigidBody.inertia.fset(self, I)
@property
def is_rigidbody(self):
if hasattr(self, '_inertia'):
return True
return False
def kinetic_energy(self, frame):
"""Kinetic energy of the body.
Parameters
==========
frame : ReferenceFrame or Body
The Body's angular velocity and the velocity of it's mass
center are typically defined with respect to an inertial frame but
any relevant frame in which the velocities are known can be supplied.
Examples
========
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy.physics.mechanics import Body, ReferenceFrame, Point
>>> from sympy import symbols
>>> m, v, r, omega = symbols('m v r omega')
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> with ignore_warnings(DeprecationWarning):
... P = Body('P', masscenter=O, mass=m)
>>> P.masscenter.set_vel(N, v * N.y)
>>> P.kinetic_energy(N)
m*v**2/2
>>> N = ReferenceFrame('N')
>>> b = ReferenceFrame('b')
>>> b.set_ang_vel(N, omega * b.x)
>>> P = Point('P')
>>> P.set_vel(N, v * N.x)
>>> with ignore_warnings(DeprecationWarning):
... B = Body('B', masscenter=P, frame=b)
>>> B.kinetic_energy(N)
B_ixx*omega**2/2 + B_mass*v**2/2
See Also
========
sympy.physics.mechanics : Particle, RigidBody
"""
if isinstance(frame, Body):
frame = Body.frame
if self.is_rigidbody:
return RigidBody(self.name, self.masscenter, self.frame, self.mass,
(self.central_inertia, self.masscenter)).kinetic_energy(frame)
return Particle(self.name, self.masscenter, self.mass).kinetic_energy(frame)
def apply_force(self, force, point=None, reaction_body=None, reaction_point=None):
"""Add force to the body(s).
Explanation
===========
Applies the force on self or equal and opposite forces on
self and other body if both are given on the desired point on the bodies.
The force applied on other body is taken opposite of self, i.e, -force.
Parameters
==========
force: Vector
The force to be applied.
point: Point, optional
The point on self on which force is applied.
By default self's masscenter.
reaction_body: Body, optional
Second body on which equal and opposite force
is to be applied.
reaction_point : Point, optional
The point on other body on which equal and opposite
force is applied. By default masscenter of other body.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy import symbols
>>> from sympy.physics.mechanics import Body, Point, dynamicsymbols
>>> m, g = symbols('m g')
>>> with ignore_warnings(DeprecationWarning):
... B = Body('B')
>>> force1 = m*g*B.z
>>> B.apply_force(force1) #Applying force on B's masscenter
>>> B.loads
[(B_masscenter, g*m*B_frame.z)]
We can also remove some part of force from any point on the body by
adding the opposite force to the body on that point.
>>> f1, f2 = dynamicsymbols('f1 f2')
>>> P = Point('P') #Considering point P on body B
>>> B.apply_force(f1*B.x + f2*B.y, P)
>>> B.loads
[(B_masscenter, g*m*B_frame.z), (P, f1(t)*B_frame.x + f2(t)*B_frame.y)]
Let's remove f1 from point P on body B.
>>> B.apply_force(-f1*B.x, P)
>>> B.loads
[(B_masscenter, g*m*B_frame.z), (P, f2(t)*B_frame.y)]
To further demonstrate the use of ``apply_force`` attribute,
consider two bodies connected through a spring.
>>> from sympy.physics.mechanics import Body, dynamicsymbols
>>> with ignore_warnings(DeprecationWarning):
... N = Body('N') #Newtonion Frame
>>> x = dynamicsymbols('x')
>>> with ignore_warnings(DeprecationWarning):
... B1 = Body('B1')
... B2 = Body('B2')
>>> spring_force = x*N.x
Now let's apply equal and opposite spring force to the bodies.
>>> P1 = Point('P1')
>>> P2 = Point('P2')
>>> B1.apply_force(spring_force, point=P1, reaction_body=B2, reaction_point=P2)
We can check the loads(forces) applied to bodies now.
>>> B1.loads
[(P1, x(t)*N_frame.x)]
>>> B2.loads
[(P2, - x(t)*N_frame.x)]
Notes
=====
If a new force is applied to a body on a point which already has some
force applied on it, then the new force is added to the already applied
force on that point.
"""
if not isinstance(point, Point):
if point is None:
point = self.masscenter # masscenter
else:
raise TypeError("Force must be applied to a point on the body.")
if not isinstance(force, Vector):
raise TypeError("Force must be a vector.")
if reaction_body is not None:
reaction_body.apply_force(-force, point=reaction_point)
for load in self._loads:
if point in load:
force += load[1]
self._loads.remove(load)
break
self._loads.append((point, force))
def apply_torque(self, torque, reaction_body=None):
"""Add torque to the body(s).
Explanation
===========
Applies the torque on self or equal and opposite torques on
self and other body if both are given.
The torque applied on other body is taken opposite of self,
i.e, -torque.
Parameters
==========
torque: Vector
The torque to be applied.
reaction_body: Body, optional
Second body on which equal and opposite torque
is to be applied.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy import symbols
>>> from sympy.physics.mechanics import Body, dynamicsymbols
>>> t = symbols('t')
>>> with ignore_warnings(DeprecationWarning):
... B = Body('B')
>>> torque1 = t*B.z
>>> B.apply_torque(torque1)
>>> B.loads
[(B_frame, t*B_frame.z)]
We can also remove some part of torque from the body by
adding the opposite torque to the body.
>>> t1, t2 = dynamicsymbols('t1 t2')
>>> B.apply_torque(t1*B.x + t2*B.y)
>>> B.loads
[(B_frame, t1(t)*B_frame.x + t2(t)*B_frame.y + t*B_frame.z)]
Let's remove t1 from Body B.
>>> B.apply_torque(-t1*B.x)
>>> B.loads
[(B_frame, t2(t)*B_frame.y + t*B_frame.z)]
To further demonstrate the use, let us consider two bodies such that
a torque `T` is acting on one body, and `-T` on the other.
>>> from sympy.physics.mechanics import Body, dynamicsymbols
>>> with ignore_warnings(DeprecationWarning):
... N = Body('N') #Newtonion frame
... B1 = Body('B1')
... B2 = Body('B2')
>>> v = dynamicsymbols('v')
>>> T = v*N.y #Torque
Now let's apply equal and opposite torque to the bodies.
>>> B1.apply_torque(T, B2)
We can check the loads (torques) applied to bodies now.
>>> B1.loads
[(B1_frame, v(t)*N_frame.y)]
>>> B2.loads
[(B2_frame, - v(t)*N_frame.y)]
Notes
=====
If a new torque is applied on body which already has some torque applied on it,
then the new torque is added to the previous torque about the body's frame.
"""
if not isinstance(torque, Vector):
raise TypeError("A Vector must be supplied to add torque.")
if reaction_body is not None:
reaction_body.apply_torque(-torque)
for load in self._loads:
if self.frame in load:
torque += load[1]
self._loads.remove(load)
break
self._loads.append((self.frame, torque))
def clear_loads(self):
"""
Clears the Body's loads list.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy.physics.mechanics import Body
>>> with ignore_warnings(DeprecationWarning):
... B = Body('B')
>>> force = B.x + B.y
>>> B.apply_force(force)
>>> B.loads
[(B_masscenter, B_frame.x + B_frame.y)]
>>> B.clear_loads()
>>> B.loads
[]
"""
self._loads = []
def remove_load(self, about=None):
"""
Remove load about a point or frame.
Parameters
==========
about : Point or ReferenceFrame, optional
The point about which force is applied,
and is to be removed.
If about is None, then the torque about
self's frame is removed.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy.physics.mechanics import Body, Point
>>> with ignore_warnings(DeprecationWarning):
... B = Body('B')
>>> P = Point('P')
>>> f1 = B.x
>>> f2 = B.y
>>> B.apply_force(f1)
>>> B.apply_force(f2, P)
>>> B.loads
[(B_masscenter, B_frame.x), (P, B_frame.y)]
>>> B.remove_load(P)
>>> B.loads
[(B_masscenter, B_frame.x)]
"""
if about is not None:
if not isinstance(about, Point):
raise TypeError('Load is applied about Point or ReferenceFrame.')
else:
about = self.frame
for load in self._loads:
if about in load:
self._loads.remove(load)
break
def masscenter_vel(self, body):
"""
Returns the velocity of the mass center with respect to the provided
rigid body or reference frame.
Parameters
==========
body: Body or ReferenceFrame
The rigid body or reference frame to calculate the velocity in.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy.physics.mechanics import Body
>>> with ignore_warnings(DeprecationWarning):
... A = Body('A')
... B = Body('B')
>>> A.masscenter.set_vel(B.frame, 5*B.frame.x)
>>> A.masscenter_vel(B)
5*B_frame.x
>>> A.masscenter_vel(B.frame)
5*B_frame.x
"""
if isinstance(body, ReferenceFrame):
frame=body
elif isinstance(body, Body):
frame = body.frame
return self.masscenter.vel(frame)
def ang_vel_in(self, body):
"""
Returns this body's angular velocity with respect to the provided
rigid body or reference frame.
Parameters
==========
body: Body or ReferenceFrame
The rigid body or reference frame to calculate the angular velocity in.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy.physics.mechanics import Body, ReferenceFrame
>>> with ignore_warnings(DeprecationWarning):
... A = Body('A')
>>> N = ReferenceFrame('N')
>>> with ignore_warnings(DeprecationWarning):
... B = Body('B', frame=N)
>>> A.frame.set_ang_vel(N, 5*N.x)
>>> A.ang_vel_in(B)
5*N.x
>>> A.ang_vel_in(N)
5*N.x
"""
if isinstance(body, ReferenceFrame):
frame=body
elif isinstance(body, Body):
frame = body.frame
return self.frame.ang_vel_in(frame)
def dcm(self, body):
"""
Returns the direction cosine matrix of this body relative to the
provided rigid body or reference frame.
Parameters
==========
body: Body or ReferenceFrame
The rigid body or reference frame to calculate the dcm.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy.physics.mechanics import Body
>>> with ignore_warnings(DeprecationWarning):
... A = Body('A')
... B = Body('B')
>>> A.frame.orient_axis(B.frame, B.frame.x, 5)
>>> A.dcm(B)
Matrix([
[1, 0, 0],
[0, cos(5), sin(5)],
[0, -sin(5), cos(5)]])
>>> A.dcm(B.frame)
Matrix([
[1, 0, 0],
[0, cos(5), sin(5)],
[0, -sin(5), cos(5)]])
"""
if isinstance(body, ReferenceFrame):
frame=body
elif isinstance(body, Body):
frame = body.frame
return self.frame.dcm(frame)
def parallel_axis(self, point, frame=None):
"""Returns the inertia dyadic of the body with respect to another
point.
Parameters
==========
point : sympy.physics.vector.Point
The point to express the inertia dyadic about.
frame : sympy.physics.vector.ReferenceFrame
The reference frame used to construct the dyadic.
Returns
=======
inertia : sympy.physics.vector.Dyadic
The inertia dyadic of the rigid body expressed about the provided
point.
Example
=======
As Body has been deprecated, the following examples are for illustrative
purposes only. The functionality of Body is fully captured by
:class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
>>> from sympy.physics.mechanics import Body
>>> with ignore_warnings(DeprecationWarning):
... A = Body('A')
>>> P = A.masscenter.locatenew('point', 3 * A.x + 5 * A.y)
>>> A.parallel_axis(P).to_matrix(A.frame)
Matrix([
[A_ixx + 25*A_mass, A_ixy - 15*A_mass, A_izx],
[A_ixy - 15*A_mass, A_iyy + 9*A_mass, A_iyz],
[ A_izx, A_iyz, A_izz + 34*A_mass]])
"""
if self.is_rigidbody:
return RigidBody.parallel_axis(self, point, frame)
return Particle.parallel_axis(self, point, frame)
@@ -0,0 +1,94 @@
from abc import ABC, abstractmethod
from sympy import Symbol, sympify
from sympy.physics.vector import Point
__all__ = ['BodyBase']
class BodyBase(ABC):
"""Abstract class for body type objects."""
def __init__(self, name, masscenter=None, mass=None):
# Note: If frame=None, no auto-generated frame is created, because a
# Particle does not need to have a frame by default.
if not isinstance(name, str):
raise TypeError('Supply a valid name.')
self._name = name
if mass is None:
mass = Symbol(f'{name}_mass')
if masscenter is None:
masscenter = Point(f'{name}_masscenter')
self.mass = mass
self.masscenter = masscenter
self.potential_energy = 0
self.points = []
def __str__(self):
return self.name
def __repr__(self):
return (f'{self.__class__.__name__}({repr(self.name)}, masscenter='
f'{repr(self.masscenter)}, mass={repr(self.mass)})')
@property
def name(self):
"""The name of the body."""
return self._name
@property
def masscenter(self):
"""The body's center of mass."""
return self._masscenter
@masscenter.setter
def masscenter(self, point):
if not isinstance(point, Point):
raise TypeError("The body's center of mass must be a Point object.")
self._masscenter = point
@property
def mass(self):
"""The body's mass."""
return self._mass
@mass.setter
def mass(self, mass):
self._mass = sympify(mass)
@property
def potential_energy(self):
"""The potential energy of the body.
Examples
========
>>> from sympy.physics.mechanics import Particle, Point
>>> from sympy import symbols
>>> m, g, h = symbols('m g h')
>>> O = Point('O')
>>> P = Particle('P', O, m)
>>> P.potential_energy = m * g * h
>>> P.potential_energy
g*h*m
"""
return self._potential_energy
@potential_energy.setter
def potential_energy(self, scalar):
self._potential_energy = sympify(scalar)
@abstractmethod
def kinetic_energy(self, frame):
pass
@abstractmethod
def linear_momentum(self, frame):
pass
@abstractmethod
def angular_momentum(self, point, frame):
pass
@abstractmethod
def parallel_axis(self, point, frame):
pass
@@ -0,0 +1,735 @@
from sympy.utilities import dict_merge
from sympy.utilities.iterables import iterable
from sympy.physics.vector import (Dyadic, Vector, ReferenceFrame,
Point, dynamicsymbols)
from sympy.physics.vector.printing import (vprint, vsprint, vpprint, vlatex,
init_vprinting)
from sympy.physics.mechanics.particle import Particle
from sympy.physics.mechanics.rigidbody import RigidBody
from sympy.simplify.simplify import simplify
from sympy import Matrix, Mul, Derivative, sin, cos, tan, S
from sympy.core.function import AppliedUndef
from sympy.physics.mechanics.inertia import (inertia as _inertia,
inertia_of_point_mass as _inertia_of_point_mass)
from sympy.utilities.exceptions import sympy_deprecation_warning
__all__ = ['linear_momentum',
'angular_momentum',
'kinetic_energy',
'potential_energy',
'Lagrangian',
'mechanics_printing',
'mprint',
'msprint',
'mpprint',
'mlatex',
'msubs',
'find_dynamicsymbols']
# These are functions that we've moved and renamed during extracting the
# basic vector calculus code from the mechanics packages.
mprint = vprint
msprint = vsprint
mpprint = vpprint
mlatex = vlatex
def mechanics_printing(**kwargs):
"""
Initializes time derivative printing for all SymPy objects in
mechanics module.
"""
init_vprinting(**kwargs)
mechanics_printing.__doc__ = init_vprinting.__doc__
def inertia(frame, ixx, iyy, izz, ixy=0, iyz=0, izx=0):
sympy_deprecation_warning(
"""
The inertia function has been moved.
Import it from "sympy.physics.mechanics".
""",
deprecated_since_version="1.13",
active_deprecations_target="moved-mechanics-functions"
)
return _inertia(frame, ixx, iyy, izz, ixy, iyz, izx)
def inertia_of_point_mass(mass, pos_vec, frame):
sympy_deprecation_warning(
"""
The inertia_of_point_mass function has been moved.
Import it from "sympy.physics.mechanics".
""",
deprecated_since_version="1.13",
active_deprecations_target="moved-mechanics-functions"
)
return _inertia_of_point_mass(mass, pos_vec, frame)
def linear_momentum(frame, *body):
"""Linear momentum of the system.
Explanation
===========
This function returns the linear momentum of a system of Particle's and/or
RigidBody's. The linear momentum of a system is equal to the vector sum of
the linear momentum of its constituents. Consider a system, S, comprised of
a rigid body, A, and a particle, P. The linear momentum of the system, L,
is equal to the vector sum of the linear momentum of the particle, L1, and
the linear momentum of the rigid body, L2, i.e.
L = L1 + L2
Parameters
==========
frame : ReferenceFrame
The frame in which linear momentum is desired.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose linear momentum is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, linear_momentum
>>> N = ReferenceFrame('N')
>>> P = Point('P')
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = Point('Ac')
>>> Ac.set_vel(N, 25 * N.y)
>>> I = outer(N.x, N.x)
>>> A = RigidBody('A', Ac, N, 20, (I, Ac))
>>> linear_momentum(N, A, Pa)
10*N.x + 500*N.y
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please specify a valid ReferenceFrame')
else:
linear_momentum_sys = Vector(0)
for e in body:
if isinstance(e, (RigidBody, Particle)):
linear_momentum_sys += e.linear_momentum(frame)
else:
raise TypeError('*body must have only Particle or RigidBody')
return linear_momentum_sys
def angular_momentum(point, frame, *body):
"""Angular momentum of a system.
Explanation
===========
This function returns the angular momentum of a system of Particle's and/or
RigidBody's. The angular momentum of such a system is equal to the vector
sum of the angular momentum of its constituents. Consider a system, S,
comprised of a rigid body, A, and a particle, P. The angular momentum of
the system, H, is equal to the vector sum of the angular momentum of the
particle, H1, and the angular momentum of the rigid body, H2, i.e.
H = H1 + H2
Parameters
==========
point : Point
The point about which angular momentum of the system is desired.
frame : ReferenceFrame
The frame in which angular momentum is desired.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose angular momentum is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, angular_momentum
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> Ac.set_vel(N, 5 * N.y)
>>> a = ReferenceFrame('a')
>>> a.set_ang_vel(N, 10 * N.z)
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, 20, (I, Ac))
>>> angular_momentum(O, N, Pa, A)
10*N.z
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please enter a valid ReferenceFrame')
if not isinstance(point, Point):
raise TypeError('Please specify a valid Point')
else:
angular_momentum_sys = Vector(0)
for e in body:
if isinstance(e, (RigidBody, Particle)):
angular_momentum_sys += e.angular_momentum(point, frame)
else:
raise TypeError('*body must have only Particle or RigidBody')
return angular_momentum_sys
def kinetic_energy(frame, *body):
"""Kinetic energy of a multibody system.
Explanation
===========
This function returns the kinetic energy of a system of Particle's and/or
RigidBody's. The kinetic energy of such a system is equal to the sum of
the kinetic energies of its constituents. Consider a system, S, comprising
a rigid body, A, and a particle, P. The kinetic energy of the system, T,
is equal to the vector sum of the kinetic energy of the particle, T1, and
the kinetic energy of the rigid body, T2, i.e.
T = T1 + T2
Kinetic energy is a scalar.
Parameters
==========
frame : ReferenceFrame
The frame in which the velocity or angular velocity of the body is
defined.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose kinetic energy is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, kinetic_energy
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> Ac.set_vel(N, 5 * N.y)
>>> a = ReferenceFrame('a')
>>> a.set_ang_vel(N, 10 * N.z)
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, 20, (I, Ac))
>>> kinetic_energy(N, Pa, A)
350
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please enter a valid ReferenceFrame')
ke_sys = S.Zero
for e in body:
if isinstance(e, (RigidBody, Particle)):
ke_sys += e.kinetic_energy(frame)
else:
raise TypeError('*body must have only Particle or RigidBody')
return ke_sys
def potential_energy(*body):
"""Potential energy of a multibody system.
Explanation
===========
This function returns the potential energy of a system of Particle's and/or
RigidBody's. The potential energy of such a system is equal to the sum of
the potential energy of its constituents. Consider a system, S, comprising
a rigid body, A, and a particle, P. The potential energy of the system, V,
is equal to the vector sum of the potential energy of the particle, V1, and
the potential energy of the rigid body, V2, i.e.
V = V1 + V2
Potential energy is a scalar.
Parameters
==========
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose potential energy is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, potential_energy
>>> from sympy import symbols
>>> M, m, g, h = symbols('M m g h')
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> Pa = Particle('Pa', P, m)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> a = ReferenceFrame('a')
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, M, (I, Ac))
>>> Pa.potential_energy = m * g * h
>>> A.potential_energy = M * g * h
>>> potential_energy(Pa, A)
M*g*h + g*h*m
"""
pe_sys = S.Zero
for e in body:
if isinstance(e, (RigidBody, Particle)):
pe_sys += e.potential_energy
else:
raise TypeError('*body must have only Particle or RigidBody')
return pe_sys
def gravity(acceleration, *bodies):
from sympy.physics.mechanics.loads import gravity as _gravity
sympy_deprecation_warning(
"""
The gravity function has been moved.
Import it from "sympy.physics.mechanics.loads".
""",
deprecated_since_version="1.13",
active_deprecations_target="moved-mechanics-functions"
)
return _gravity(acceleration, *bodies)
def center_of_mass(point, *bodies):
"""
Returns the position vector from the given point to the center of mass
of the given bodies(particles or rigidbodies).
Example
=======
>>> from sympy import symbols, S
>>> from sympy.physics.vector import Point
>>> from sympy.physics.mechanics import Particle, ReferenceFrame, RigidBody, outer
>>> from sympy.physics.mechanics.functions import center_of_mass
>>> a = ReferenceFrame('a')
>>> m = symbols('m', real=True)
>>> p1 = Particle('p1', Point('p1_pt'), S(1))
>>> p2 = Particle('p2', Point('p2_pt'), S(2))
>>> p3 = Particle('p3', Point('p3_pt'), S(3))
>>> p4 = Particle('p4', Point('p4_pt'), m)
>>> b_f = ReferenceFrame('b_f')
>>> b_cm = Point('b_cm')
>>> mb = symbols('mb')
>>> b = RigidBody('b', b_cm, b_f, mb, (outer(b_f.x, b_f.x), b_cm))
>>> p2.point.set_pos(p1.point, a.x)
>>> p3.point.set_pos(p1.point, a.x + a.y)
>>> p4.point.set_pos(p1.point, a.y)
>>> b.masscenter.set_pos(p1.point, a.y + a.z)
>>> point_o=Point('o')
>>> point_o.set_pos(p1.point, center_of_mass(p1.point, p1, p2, p3, p4, b))
>>> expr = 5/(m + mb + 6)*a.x + (m + mb + 3)/(m + mb + 6)*a.y + mb/(m + mb + 6)*a.z
>>> point_o.pos_from(p1.point)
5/(m + mb + 6)*a.x + (m + mb + 3)/(m + mb + 6)*a.y + mb/(m + mb + 6)*a.z
"""
if not bodies:
raise TypeError("No bodies(instances of Particle or Rigidbody) were passed.")
total_mass = 0
vec = Vector(0)
for i in bodies:
total_mass += i.mass
masscenter = getattr(i, 'masscenter', None)
if masscenter is None:
masscenter = i.point
vec += i.mass*masscenter.pos_from(point)
return vec/total_mass
def Lagrangian(frame, *body):
"""Lagrangian of a multibody system.
Explanation
===========
This function returns the Lagrangian of a system of Particle's and/or
RigidBody's. The Lagrangian of such a system is equal to the difference
between the kinetic energies and potential energies of its constituents. If
T and V are the kinetic and potential energies of a system then it's
Lagrangian, L, is defined as
L = T - V
The Lagrangian is a scalar.
Parameters
==========
frame : ReferenceFrame
The frame in which the velocity or angular velocity of the body is
defined to determine the kinetic energy.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose Lagrangian is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, Lagrangian
>>> from sympy import symbols
>>> M, m, g, h = symbols('M m g h')
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> Ac.set_vel(N, 5 * N.y)
>>> a = ReferenceFrame('a')
>>> a.set_ang_vel(N, 10 * N.z)
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, 20, (I, Ac))
>>> Pa.potential_energy = m * g * h
>>> A.potential_energy = M * g * h
>>> Lagrangian(N, Pa, A)
-M*g*h - g*h*m + 350
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please supply a valid ReferenceFrame')
for e in body:
if not isinstance(e, (RigidBody, Particle)):
raise TypeError('*body must have only Particle or RigidBody')
return kinetic_energy(frame, *body) - potential_energy(*body)
def find_dynamicsymbols(expression, exclude=None, reference_frame=None):
"""Find all dynamicsymbols in expression.
Explanation
===========
If the optional ``exclude`` kwarg is used, only dynamicsymbols
not in the iterable ``exclude`` are returned.
If we intend to apply this function on a vector, the optional
``reference_frame`` is also used to inform about the corresponding frame
with respect to which the dynamic symbols of the given vector is to be
determined.
Parameters
==========
expression : SymPy expression
exclude : iterable of dynamicsymbols, optional
reference_frame : ReferenceFrame, optional
The frame with respect to which the dynamic symbols of the
given vector is to be determined.
Examples
========
>>> from sympy.physics.mechanics import dynamicsymbols, find_dynamicsymbols
>>> from sympy.physics.mechanics import ReferenceFrame
>>> x, y = dynamicsymbols('x, y')
>>> expr = x + x.diff()*y
>>> find_dynamicsymbols(expr)
{x(t), y(t), Derivative(x(t), t)}
>>> find_dynamicsymbols(expr, exclude=[x, y])
{Derivative(x(t), t)}
>>> a, b, c = dynamicsymbols('a, b, c')
>>> A = ReferenceFrame('A')
>>> v = a * A.x + b * A.y + c * A.z
>>> find_dynamicsymbols(v, reference_frame=A)
{a(t), b(t), c(t)}
"""
t_set = {dynamicsymbols._t}
if exclude:
if iterable(exclude):
exclude_set = set(exclude)
else:
raise TypeError("exclude kwarg must be iterable")
else:
exclude_set = set()
if isinstance(expression, Vector):
if reference_frame is None:
raise ValueError("You must provide reference_frame when passing a "
"vector expression, got %s." % reference_frame)
else:
expression = expression.to_matrix(reference_frame)
return {i for i in expression.atoms(AppliedUndef, Derivative) if
i.free_symbols == t_set} - exclude_set
def msubs(expr, *sub_dicts, smart=False, **kwargs):
"""A custom subs for use on expressions derived in physics.mechanics.
Traverses the expression tree once, performing the subs found in sub_dicts.
Terms inside ``Derivative`` expressions are ignored:
Examples
========
>>> from sympy.physics.mechanics import dynamicsymbols, msubs
>>> x = dynamicsymbols('x')
>>> msubs(x.diff() + x, {x: 1})
Derivative(x(t), t) + 1
Note that sub_dicts can be a single dictionary, or several dictionaries:
>>> x, y, z = dynamicsymbols('x, y, z')
>>> sub1 = {x: 1, y: 2}
>>> sub2 = {z: 3, x.diff(): 4}
>>> msubs(x.diff() + x + y + z, sub1, sub2)
10
If smart=True (default False), also checks for conditions that may result
in ``nan``, but if simplified would yield a valid expression. For example:
>>> from sympy import sin, tan
>>> (sin(x)/tan(x)).subs(x, 0)
nan
>>> msubs(sin(x)/tan(x), {x: 0}, smart=True)
1
It does this by first replacing all ``tan`` with ``sin/cos``. Then each
node is traversed. If the node is a fraction, subs is first evaluated on
the denominator. If this results in 0, simplification of the entire
fraction is attempted. Using this selective simplification, only
subexpressions that result in 1/0 are targeted, resulting in faster
performance.
"""
sub_dict = dict_merge(*sub_dicts)
if smart:
func = _smart_subs
elif hasattr(expr, 'msubs'):
return expr.msubs(sub_dict)
else:
func = lambda expr, sub_dict: _crawl(expr, _sub_func, sub_dict)
if isinstance(expr, (Matrix, Vector, Dyadic)):
return expr.applyfunc(lambda x: func(x, sub_dict))
else:
return func(expr, sub_dict)
def _crawl(expr, func, *args, **kwargs):
"""Crawl the expression tree, and apply func to every node."""
val = func(expr, *args, **kwargs)
if val is not None:
return val
new_args = (_crawl(arg, func, *args, **kwargs) for arg in expr.args)
return expr.func(*new_args)
def _sub_func(expr, sub_dict):
"""Perform direct matching substitution, ignoring derivatives."""
if expr in sub_dict:
return sub_dict[expr]
elif not expr.args or expr.is_Derivative:
return expr
def _tan_repl_func(expr):
"""Replace tan with sin/cos."""
if isinstance(expr, tan):
return sin(*expr.args) / cos(*expr.args)
elif not expr.args or expr.is_Derivative:
return expr
def _smart_subs(expr, sub_dict):
"""Performs subs, checking for conditions that may result in `nan` or
`oo`, and attempts to simplify them out.
The expression tree is traversed twice, and the following steps are
performed on each expression node:
- First traverse:
Replace all `tan` with `sin/cos`.
- Second traverse:
If node is a fraction, check if the denominator evaluates to 0.
If so, attempt to simplify it out. Then if node is in sub_dict,
sub in the corresponding value.
"""
expr = _crawl(expr, _tan_repl_func)
def _recurser(expr, sub_dict):
# Decompose the expression into num, den
num, den = _fraction_decomp(expr)
if den != 1:
# If there is a non trivial denominator, we need to handle it
denom_subbed = _recurser(den, sub_dict)
if denom_subbed.evalf() == 0:
# If denom is 0 after this, attempt to simplify the bad expr
expr = simplify(expr)
else:
# Expression won't result in nan, find numerator
num_subbed = _recurser(num, sub_dict)
return num_subbed / denom_subbed
# We have to crawl the tree manually, because `expr` may have been
# modified in the simplify step. First, perform subs as normal:
val = _sub_func(expr, sub_dict)
if val is not None:
return val
new_args = (_recurser(arg, sub_dict) for arg in expr.args)
return expr.func(*new_args)
return _recurser(expr, sub_dict)
def _fraction_decomp(expr):
"""Return num, den such that expr = num/den."""
if not isinstance(expr, Mul):
return expr, 1
num = []
den = []
for a in expr.args:
if a.is_Pow and a.args[1] < 0:
den.append(1 / a)
else:
num.append(a)
if not den:
return expr, 1
num = Mul(*num)
den = Mul(*den)
return num, den
def _f_list_parser(fl, ref_frame):
"""Parses the provided forcelist composed of items
of the form (obj, force).
Returns a tuple containing:
vel_list: The velocity (ang_vel for Frames, vel for Points) in
the provided reference frame.
f_list: The forces.
Used internally in the KanesMethod and LagrangesMethod classes.
"""
def flist_iter():
for pair in fl:
obj, force = pair
if isinstance(obj, ReferenceFrame):
yield obj.ang_vel_in(ref_frame), force
elif isinstance(obj, Point):
yield obj.vel(ref_frame), force
else:
raise TypeError('First entry in each forcelist pair must '
'be a point or frame.')
if not fl:
vel_list, f_list = (), ()
else:
unzip = lambda l: list(zip(*l)) if l[0] else [(), ()]
vel_list, f_list = unzip(list(flist_iter()))
return vel_list, f_list
def _validate_coordinates(coordinates=None, speeds=None, check_duplicates=True,
is_dynamicsymbols=True, u_auxiliary=None):
"""Validate the generalized coordinates and generalized speeds.
Parameters
==========
coordinates : iterable, optional
Generalized coordinates to be validated.
speeds : iterable, optional
Generalized speeds to be validated.
check_duplicates : bool, optional
Checks if there are duplicates in the generalized coordinates and
generalized speeds. If so it will raise a ValueError. The default is
True.
is_dynamicsymbols : iterable, optional
Checks if all the generalized coordinates and generalized speeds are
dynamicsymbols. If any is not a dynamicsymbol, a ValueError will be
raised. The default is True.
u_auxiliary : iterable, optional
Auxiliary generalized speeds to be validated.
"""
t_set = {dynamicsymbols._t}
# Convert input to iterables
if coordinates is None:
coordinates = []
elif not iterable(coordinates):
coordinates = [coordinates]
if speeds is None:
speeds = []
elif not iterable(speeds):
speeds = [speeds]
if u_auxiliary is None:
u_auxiliary = []
elif not iterable(u_auxiliary):
u_auxiliary = [u_auxiliary]
msgs = []
if check_duplicates: # Check for duplicates
seen = set()
coord_duplicates = {x for x in coordinates if x in seen or seen.add(x)}
seen = set()
speed_duplicates = {x for x in speeds if x in seen or seen.add(x)}
seen = set()
aux_duplicates = {x for x in u_auxiliary if x in seen or seen.add(x)}
overlap_coords = set(coordinates).intersection(speeds)
overlap_aux = set(coordinates).union(speeds).intersection(u_auxiliary)
if coord_duplicates:
msgs.append(f'The generalized coordinates {coord_duplicates} are '
f'duplicated, all generalized coordinates should be '
f'unique.')
if speed_duplicates:
msgs.append(f'The generalized speeds {speed_duplicates} are '
f'duplicated, all generalized speeds should be unique.')
if aux_duplicates:
msgs.append(f'The auxiliary speeds {aux_duplicates} are duplicated,'
f' all auxiliary speeds should be unique.')
if overlap_coords:
msgs.append(f'{overlap_coords} are defined as both generalized '
f'coordinates and generalized speeds.')
if overlap_aux:
msgs.append(f'The auxiliary speeds {overlap_aux} are also defined '
f'as generalized coordinates or generalized speeds.')
if is_dynamicsymbols: # Check whether all coordinates are dynamicsymbols
for coordinate in coordinates:
if not (isinstance(coordinate, (AppliedUndef, Derivative)) and
coordinate.free_symbols == t_set):
msgs.append(f'Generalized coordinate "{coordinate}" is not a '
f'dynamicsymbol.')
for speed in speeds:
if not (isinstance(speed, (AppliedUndef, Derivative)) and
speed.free_symbols == t_set):
msgs.append(
f'Generalized speed "{speed}" is not a dynamicsymbol.')
for aux in u_auxiliary:
if not (isinstance(aux, (AppliedUndef, Derivative)) and
aux.free_symbols == t_set):
msgs.append(
f'Auxiliary speed "{aux}" is not a dynamicsymbol.')
if msgs:
raise ValueError('\n'.join(msgs))
def _parse_linear_solver(linear_solver):
"""Helper function to retrieve a specified linear solver."""
if callable(linear_solver):
return linear_solver
return lambda A, b: Matrix.solve(A, b, method=linear_solver)
@@ -0,0 +1,199 @@
from sympy import sympify
from sympy.physics.vector import Point, Dyadic, ReferenceFrame, outer
from collections import namedtuple
__all__ = ['inertia', 'inertia_of_point_mass', 'Inertia']
def inertia(frame, ixx, iyy, izz, ixy=0, iyz=0, izx=0):
"""Simple way to create inertia Dyadic object.
Explanation
===========
Creates an inertia Dyadic based on the given tensor values and a body-fixed
reference frame.
Parameters
==========
frame : ReferenceFrame
The frame the inertia is defined in.
ixx : Sympifyable
The xx element in the inertia dyadic.
iyy : Sympifyable
The yy element in the inertia dyadic.
izz : Sympifyable
The zz element in the inertia dyadic.
ixy : Sympifyable
The xy element in the inertia dyadic.
iyz : Sympifyable
The yz element in the inertia dyadic.
izx : Sympifyable
The zx element in the inertia dyadic.
Examples
========
>>> from sympy.physics.mechanics import ReferenceFrame, inertia
>>> N = ReferenceFrame('N')
>>> inertia(N, 1, 2, 3)
(N.x|N.x) + 2*(N.y|N.y) + 3*(N.z|N.z)
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Need to define the inertia in a frame')
ixx, iyy, izz = sympify(ixx), sympify(iyy), sympify(izz)
ixy, iyz, izx = sympify(ixy), sympify(iyz), sympify(izx)
return (ixx*outer(frame.x, frame.x) + ixy*outer(frame.x, frame.y) +
izx*outer(frame.x, frame.z) + ixy*outer(frame.y, frame.x) +
iyy*outer(frame.y, frame.y) + iyz*outer(frame.y, frame.z) +
izx*outer(frame.z, frame.x) + iyz*outer(frame.z, frame.y) +
izz*outer(frame.z, frame.z))
def inertia_of_point_mass(mass, pos_vec, frame):
"""Inertia dyadic of a point mass relative to point O.
Parameters
==========
mass : Sympifyable
Mass of the point mass
pos_vec : Vector
Position from point O to point mass
frame : ReferenceFrame
Reference frame to express the dyadic in
Examples
========
>>> from sympy import symbols
>>> from sympy.physics.mechanics import ReferenceFrame, inertia_of_point_mass
>>> N = ReferenceFrame('N')
>>> r, m = symbols('r m')
>>> px = r * N.x
>>> inertia_of_point_mass(m, px, N)
m*r**2*(N.y|N.y) + m*r**2*(N.z|N.z)
"""
return mass*(
(outer(frame.x, frame.x) +
outer(frame.y, frame.y) +
outer(frame.z, frame.z)) *
(pos_vec.dot(pos_vec)) - outer(pos_vec, pos_vec))
class Inertia(namedtuple('Inertia', ['dyadic', 'point'])):
"""Inertia object consisting of a Dyadic and a Point of reference.
Explanation
===========
This is a simple class to store the Point and Dyadic, belonging to an
inertia.
Attributes
==========
dyadic : Dyadic
The dyadic of the inertia.
point : Point
The reference point of the inertia.
Examples
========
>>> from sympy.physics.mechanics import ReferenceFrame, Point, Inertia
>>> N = ReferenceFrame('N')
>>> Po = Point('Po')
>>> Inertia(N.x.outer(N.x) + N.y.outer(N.y) + N.z.outer(N.z), Po)
((N.x|N.x) + (N.y|N.y) + (N.z|N.z), Po)
In the example above the Dyadic was created manually, one can however also
use the ``inertia`` function for this or the class method ``from_tensor`` as
shown below.
>>> Inertia.from_inertia_scalars(Po, N, 1, 1, 1)
((N.x|N.x) + (N.y|N.y) + (N.z|N.z), Po)
"""
__slots__ = ()
def __new__(cls, dyadic, point):
# Switch order if given in the wrong order
if isinstance(dyadic, Point) and isinstance(point, Dyadic):
point, dyadic = dyadic, point
if not isinstance(point, Point):
raise TypeError('Reference point should be of type Point')
if not isinstance(dyadic, Dyadic):
raise TypeError('Inertia value should be expressed as a Dyadic')
return super().__new__(cls, dyadic, point)
@classmethod
def from_inertia_scalars(cls, point, frame, ixx, iyy, izz, ixy=0, iyz=0,
izx=0):
"""Simple way to create an Inertia object based on the tensor values.
Explanation
===========
This class method uses the :func`~.inertia` to create the Dyadic based
on the tensor values.
Parameters
==========
point : Point
The reference point of the inertia.
frame : ReferenceFrame
The frame the inertia is defined in.
ixx : Sympifyable
The xx element in the inertia dyadic.
iyy : Sympifyable
The yy element in the inertia dyadic.
izz : Sympifyable
The zz element in the inertia dyadic.
ixy : Sympifyable
The xy element in the inertia dyadic.
iyz : Sympifyable
The yz element in the inertia dyadic.
izx : Sympifyable
The zx element in the inertia dyadic.
Examples
========
>>> from sympy import symbols
>>> from sympy.physics.mechanics import ReferenceFrame, Point, Inertia
>>> ixx, iyy, izz, ixy, iyz, izx = symbols('ixx iyy izz ixy iyz izx')
>>> N = ReferenceFrame('N')
>>> P = Point('P')
>>> I = Inertia.from_inertia_scalars(P, N, ixx, iyy, izz, ixy, iyz, izx)
The tensor values can easily be seen when converting the dyadic to a
matrix.
>>> I.dyadic.to_matrix(N)
Matrix([
[ixx, ixy, izx],
[ixy, iyy, iyz],
[izx, iyz, izz]])
"""
return cls(inertia(frame, ixx, iyy, izz, ixy, iyz, izx), point)
def __add__(self, other):
raise TypeError(f"unsupported operand type(s) for +: "
f"'{self.__class__.__name__}' and "
f"'{other.__class__.__name__}'")
def __mul__(self, other):
raise TypeError(f"unsupported operand type(s) for *: "
f"'{self.__class__.__name__}' and "
f"'{other.__class__.__name__}'")
__radd__ = __add__
__rmul__ = __mul__
@@ -0,0 +1,318 @@
from sympy.physics.mechanics import (Body, Lagrangian, KanesMethod, LagrangesMethod,
RigidBody, Particle)
from sympy.physics.mechanics.body_base import BodyBase
from sympy.physics.mechanics.method import _Methods
from sympy import Matrix
from sympy.utilities.exceptions import sympy_deprecation_warning
__all__ = ['JointsMethod']
class JointsMethod(_Methods):
"""Method for formulating the equations of motion using a set of interconnected bodies with joints.
.. deprecated:: 1.13
The JointsMethod class is deprecated. Its functionality has been
replaced by the new :class:`~.System` class.
Parameters
==========
newtonion : Body or ReferenceFrame
The newtonion(inertial) frame.
*joints : Joint
The joints in the system
Attributes
==========
q, u : iterable
Iterable of the generalized coordinates and speeds
bodies : iterable
Iterable of Body objects in the system.
loads : iterable
Iterable of (Point, vector) or (ReferenceFrame, vector) tuples
describing the forces on the system.
mass_matrix : Matrix, shape(n, n)
The system's mass matrix
forcing : Matrix, shape(n, 1)
The system's forcing vector
mass_matrix_full : Matrix, shape(2*n, 2*n)
The "mass matrix" for the u's and q's
forcing_full : Matrix, shape(2*n, 1)
The "forcing vector" for the u's and q's
method : KanesMethod or Lagrange's method
Method's object.
kdes : iterable
Iterable of kde in they system.
Examples
========
As Body and JointsMethod have been deprecated, the following examples are
for illustrative purposes only. The functionality of Body is fully captured
by :class:`~.RigidBody` and :class:`~.Particle` and the functionality of
JointsMethod is fully captured by :class:`~.System`. To ignore the
deprecation warning we can use the ignore_warnings context manager.
>>> from sympy.utilities.exceptions import ignore_warnings
This is a simple example for a one degree of freedom translational
spring-mass-damper.
>>> from sympy import symbols
>>> from sympy.physics.mechanics import Body, JointsMethod, PrismaticJoint
>>> from sympy.physics.vector import dynamicsymbols
>>> c, k = symbols('c k')
>>> x, v = dynamicsymbols('x v')
>>> with ignore_warnings(DeprecationWarning):
... wall = Body('W')
... body = Body('B')
>>> J = PrismaticJoint('J', wall, body, coordinates=x, speeds=v)
>>> wall.apply_force(c*v*wall.x, reaction_body=body)
>>> wall.apply_force(k*x*wall.x, reaction_body=body)
>>> with ignore_warnings(DeprecationWarning):
... method = JointsMethod(wall, J)
>>> method.form_eoms()
Matrix([[-B_mass*Derivative(v(t), t) - c*v(t) - k*x(t)]])
>>> M = method.mass_matrix_full
>>> F = method.forcing_full
>>> rhs = M.LUsolve(F)
>>> rhs
Matrix([
[ v(t)],
[(-c*v(t) - k*x(t))/B_mass]])
Notes
=====
``JointsMethod`` currently only works with systems that do not have any
configuration or motion constraints.
"""
def __init__(self, newtonion, *joints):
sympy_deprecation_warning(
"""
The JointsMethod class is deprecated.
Its functionality has been replaced by the new System class.
""",
deprecated_since_version="1.13",
active_deprecations_target="deprecated-mechanics-jointsmethod"
)
if isinstance(newtonion, BodyBase):
self.frame = newtonion.frame
else:
self.frame = newtonion
self._joints = joints
self._bodies = self._generate_bodylist()
self._loads = self._generate_loadlist()
self._q = self._generate_q()
self._u = self._generate_u()
self._kdes = self._generate_kdes()
self._method = None
@property
def bodies(self):
"""List of bodies in they system."""
return self._bodies
@property
def loads(self):
"""List of loads on the system."""
return self._loads
@property
def q(self):
"""List of the generalized coordinates."""
return self._q
@property
def u(self):
"""List of the generalized speeds."""
return self._u
@property
def kdes(self):
"""List of the generalized coordinates."""
return self._kdes
@property
def forcing_full(self):
"""The "forcing vector" for the u's and q's."""
return self.method.forcing_full
@property
def mass_matrix_full(self):
"""The "mass matrix" for the u's and q's."""
return self.method.mass_matrix_full
@property
def mass_matrix(self):
"""The system's mass matrix."""
return self.method.mass_matrix
@property
def forcing(self):
"""The system's forcing vector."""
return self.method.forcing
@property
def method(self):
"""Object of method used to form equations of systems."""
return self._method
def _generate_bodylist(self):
bodies = []
for joint in self._joints:
if joint.child not in bodies:
bodies.append(joint.child)
if joint.parent not in bodies:
bodies.append(joint.parent)
return bodies
def _generate_loadlist(self):
load_list = []
for body in self.bodies:
if isinstance(body, Body):
load_list.extend(body.loads)
return load_list
def _generate_q(self):
q_ind = []
for joint in self._joints:
for coordinate in joint.coordinates:
if coordinate in q_ind:
raise ValueError('Coordinates of joints should be unique.')
q_ind.append(coordinate)
return Matrix(q_ind)
def _generate_u(self):
u_ind = []
for joint in self._joints:
for speed in joint.speeds:
if speed in u_ind:
raise ValueError('Speeds of joints should be unique.')
u_ind.append(speed)
return Matrix(u_ind)
def _generate_kdes(self):
kd_ind = Matrix(1, 0, []).T
for joint in self._joints:
kd_ind = kd_ind.col_join(joint.kdes)
return kd_ind
def _convert_bodies(self):
# Convert `Body` to `Particle` and `RigidBody`
bodylist = []
for body in self.bodies:
if not isinstance(body, Body):
bodylist.append(body)
continue
if body.is_rigidbody:
rb = RigidBody(body.name, body.masscenter, body.frame, body.mass,
(body.central_inertia, body.masscenter))
rb.potential_energy = body.potential_energy
bodylist.append(rb)
else:
part = Particle(body.name, body.masscenter, body.mass)
part.potential_energy = body.potential_energy
bodylist.append(part)
return bodylist
def form_eoms(self, method=KanesMethod):
"""Method to form system's equation of motions.
Parameters
==========
method : Class
Class name of method.
Returns
========
Matrix
Vector of equations of motions.
Examples
========
As Body and JointsMethod have been deprecated, the following examples
are for illustrative purposes only. The functionality of Body is fully
captured by :class:`~.RigidBody` and :class:`~.Particle` and the
functionality of JointsMethod is fully captured by :class:`~.System`. To
ignore the deprecation warning we can use the ignore_warnings context
manager.
>>> from sympy.utilities.exceptions import ignore_warnings
This is a simple example for a one degree of freedom translational
spring-mass-damper.
>>> from sympy import S, symbols
>>> from sympy.physics.mechanics import LagrangesMethod, dynamicsymbols, Body
>>> from sympy.physics.mechanics import PrismaticJoint, JointsMethod
>>> q = dynamicsymbols('q')
>>> qd = dynamicsymbols('q', 1)
>>> m, k, b = symbols('m k b')
>>> with ignore_warnings(DeprecationWarning):
... wall = Body('W')
... part = Body('P', mass=m)
>>> part.potential_energy = k * q**2 / S(2)
>>> J = PrismaticJoint('J', wall, part, coordinates=q, speeds=qd)
>>> wall.apply_force(b * qd * wall.x, reaction_body=part)
>>> with ignore_warnings(DeprecationWarning):
... method = JointsMethod(wall, J)
>>> method.form_eoms(LagrangesMethod)
Matrix([[b*Derivative(q(t), t) + k*q(t) + m*Derivative(q(t), (t, 2))]])
We can also solve for the states using the 'rhs' method.
>>> method.rhs()
Matrix([
[ Derivative(q(t), t)],
[(-b*Derivative(q(t), t) - k*q(t))/m]])
"""
bodylist = self._convert_bodies()
if issubclass(method, LagrangesMethod): #LagrangesMethod or similar
L = Lagrangian(self.frame, *bodylist)
self._method = method(L, self.q, self.loads, bodylist, self.frame)
else: #KanesMethod or similar
self._method = method(self.frame, q_ind=self.q, u_ind=self.u, kd_eqs=self.kdes,
forcelist=self.loads, bodies=bodylist)
soln = self.method._form_eoms()
return soln
def rhs(self, inv_method=None):
"""Returns equations that can be solved numerically.
Parameters
==========
inv_method : str
The specific sympy inverse matrix calculation method to use. For a
list of valid methods, see
:meth:`~sympy.matrices.matrixbase.MatrixBase.inv`
Returns
========
Matrix
Numerically solvable equations.
See Also
========
sympy.physics.mechanics.kane.KanesMethod.rhs:
KanesMethod's rhs function.
sympy.physics.mechanics.lagrange.LagrangesMethod.rhs:
LagrangesMethod's rhs function.
"""
return self.method.rhs(inv_method=inv_method)
@@ -0,0 +1,859 @@
from sympy import zeros, Matrix, diff, eye, linear_eq_to_matrix
from sympy.core.sorting import default_sort_key
from sympy.physics.vector import (ReferenceFrame, dynamicsymbols,
partial_velocity)
from sympy.physics.mechanics.method import _Methods
from sympy.physics.mechanics.particle import Particle
from sympy.physics.mechanics.rigidbody import RigidBody
from sympy.physics.mechanics.functions import (msubs, find_dynamicsymbols,
_f_list_parser,
_validate_coordinates,
_parse_linear_solver)
from sympy.physics.mechanics.linearize import Linearizer
from sympy.utilities.iterables import iterable
__all__ = ['KanesMethod']
class KanesMethod(_Methods):
r"""Kane's method object.
Explanation
===========
This object is used to do the "book-keeping" as you go through and form
equations of motion in the way Kane presents in:
Kane, T., Levinson, D. Dynamics Theory and Applications. 1985 McGraw-Hill
The attributes are for equations in the form [M] udot = forcing.
Attributes
==========
q, u : Matrix
Matrices of the generalized coordinates and speeds
bodies : iterable
Iterable of Particle and RigidBody objects in the system.
loads : iterable
Iterable of (Point, vector) or (ReferenceFrame, vector) tuples
describing the forces on the system.
auxiliary_eqs : Matrix
If applicable, the set of auxiliary Kane's
equations used to solve for non-contributing
forces.
mass_matrix : Matrix
The system's dynamics mass matrix: [k_d; k_dnh]
forcing : Matrix
The system's dynamics forcing vector: -[f_d; f_dnh]
mass_matrix_kin : Matrix
The "mass matrix" for kinematic differential equations: k_kqdot
forcing_kin : Matrix
The forcing vector for kinematic differential equations: -(k_ku*u + f_k)
mass_matrix_full : Matrix
The "mass matrix" for the u's and q's with dynamics and kinematics
forcing_full : Matrix
The "forcing vector" for the u's and q's with dynamics and kinematics
Parameters
==========
frame : ReferenceFrame
The inertial reference frame for the system.
q_ind : iterable of dynamicsymbols
Independent generalized coordinates.
u_ind : iterable of dynamicsymbols
Independent generalized speeds.
kd_eqs : iterable of Expr, optional
Kinematic differential equations, which linearly relate the generalized
speeds to the time-derivatives of the generalized coordinates.
q_dependent : iterable of dynamicsymbols, optional
Dependent generalized coordinates.
configuration_constraints : iterable of Expr, optional
Constraints on the system's configuration, i.e. holonomic constraints.
u_dependent : iterable of dynamicsymbols, optional
Dependent generalized speeds.
velocity_constraints : iterable of Expr, optional
Constraints on the system's velocity, i.e. the combination of the
nonholonomic constraints and the time-derivative of the holonomic
constraints.
acceleration_constraints : iterable of Expr, optional
Constraints on the system's acceleration, by default these are the
time-derivative of the velocity constraints.
u_auxiliary : iterable of dynamicsymbols, optional
Auxiliary generalized speeds.
bodies : iterable of Particle and/or RigidBody, optional
The particles and rigid bodies in the system.
forcelist : iterable of tuple[Point | ReferenceFrame, Vector], optional
Forces and torques applied on the system.
explicit_kinematics : bool
Boolean whether the mass matrices and forcing vectors should use the
explicit form (default) or implicit form for kinematics.
See the notes for more details.
kd_eqs_solver : str, callable
Method used to solve the kinematic differential equations. If a string
is supplied, it should be a valid method that can be used with the
:meth:`sympy.matrices.matrixbase.MatrixBase.solve`. If a callable is
supplied, it should have the format ``f(A, rhs)``, where it solves the
equations and returns the solution. The default utilizes LU solve. See
the notes for more information.
constraint_solver : str, callable
Method used to solve the velocity constraints. If a string is
supplied, it should be a valid method that can be used with the
:meth:`sympy.matrices.matrixbase.MatrixBase.solve`. If a callable is
supplied, it should have the format ``f(A, rhs)``, where it solves the
equations and returns the solution. The default utilizes LU solve. See
the notes for more information.
Notes
=====
The mass matrices and forcing vectors related to kinematic equations
are given in the explicit form by default. In other words, the kinematic
mass matrix is $\mathbf{k_{k\dot{q}}} = \mathbf{I}$.
In order to get the implicit form of those matrices/vectors, you can set the
``explicit_kinematics`` attribute to ``False``. So $\mathbf{k_{k\dot{q}}}$
is not necessarily an identity matrix. This can provide more compact
equations for non-simple kinematics.
Two linear solvers can be supplied to ``KanesMethod``: one for solving the
kinematic differential equations and one to solve the velocity constraints.
Both of these sets of equations can be expressed as a linear system ``Ax = rhs``,
which have to be solved in order to obtain the equations of motion.
The default solver ``'LU'``, which stands for LU solve, results relatively low
number of operations. The weakness of this method is that it can result in zero
division errors.
If zero divisions are encountered, a possible solver which may solve the problem
is ``"CRAMER"``. This method uses Cramer's rule to solve the system. This method
is slower and results in more operations than the default solver. However it only
uses a single division by default per entry of the solution.
While a valid list of solvers can be found at
:meth:`sympy.matrices.matrixbase.MatrixBase.solve`, it is also possible to supply a
`callable`. This way it is possible to use a different solver routine. If the
kinematic differential equations are not too complex it can be worth it to simplify
the solution by using ``lambda A, b: simplify(Matrix.LUsolve(A, b))``. Another
option solver one may use is :func:`sympy.solvers.solveset.linsolve`. This can be
done using `lambda A, b: tuple(linsolve((A, b)))[0]`, where we select the first
solution as our system should have only one unique solution.
Examples
========
This is a simple example for a one degree of freedom translational
spring-mass-damper.
In this example, we first need to do the kinematics.
This involves creating generalized speeds and coordinates and their
derivatives.
Then we create a point and set its velocity in a frame.
>>> from sympy import symbols
>>> from sympy.physics.mechanics import dynamicsymbols, ReferenceFrame
>>> from sympy.physics.mechanics import Point, Particle, KanesMethod
>>> q, u = dynamicsymbols('q u')
>>> qd, ud = dynamicsymbols('q u', 1)
>>> m, c, k = symbols('m c k')
>>> N = ReferenceFrame('N')
>>> P = Point('P')
>>> P.set_vel(N, u * N.x)
Next we need to arrange/store information in the way that KanesMethod
requires. The kinematic differential equations should be an iterable of
expressions. A list of forces/torques must be constructed, where each entry
in the list is a (Point, Vector) or (ReferenceFrame, Vector) tuple, where
the Vectors represent the Force or Torque.
Next a particle needs to be created, and it needs to have a point and mass
assigned to it.
Finally, a list of all bodies and particles needs to be created.
>>> kd = [qd - u]
>>> FL = [(P, (-k * q - c * u) * N.x)]
>>> pa = Particle('pa', P, m)
>>> BL = [pa]
Finally we can generate the equations of motion.
First we create the KanesMethod object and supply an inertial frame,
coordinates, generalized speeds, and the kinematic differential equations.
Additional quantities such as configuration and motion constraints,
dependent coordinates and speeds, and auxiliary speeds are also supplied
here (see the online documentation).
Next we form FR* and FR to complete: Fr + Fr* = 0.
We have the equations of motion at this point.
It makes sense to rearrange them though, so we calculate the mass matrix and
the forcing terms, for E.o.M. in the form: [MM] udot = forcing, where MM is
the mass matrix, udot is a vector of the time derivatives of the
generalized speeds, and forcing is a vector representing "forcing" terms.
>>> KM = KanesMethod(N, q_ind=[q], u_ind=[u], kd_eqs=kd)
>>> (fr, frstar) = KM.kanes_equations(BL, FL)
>>> MM = KM.mass_matrix
>>> forcing = KM.forcing
>>> rhs = MM.inv() * forcing
>>> rhs
Matrix([[(-c*u(t) - k*q(t))/m]])
>>> KM.linearize(A_and_B=True)[0]
Matrix([
[ 0, 1],
[-k/m, -c/m]])
Please look at the documentation pages for more information on how to
perform linearization and how to deal with dependent coordinates & speeds,
and how do deal with bringing non-contributing forces into evidence.
"""
def __init__(self, frame, q_ind, u_ind, kd_eqs=None, q_dependent=None,
configuration_constraints=None, u_dependent=None,
velocity_constraints=None, acceleration_constraints=None,
u_auxiliary=None, bodies=None, forcelist=None,
explicit_kinematics=True, kd_eqs_solver='LU',
constraint_solver='LU'):
"""Please read the online documentation. """
if not q_ind:
q_ind = [dynamicsymbols('dummy_q')]
kd_eqs = [dynamicsymbols('dummy_kd')]
if not isinstance(frame, ReferenceFrame):
raise TypeError('An inertial ReferenceFrame must be supplied')
self._inertial = frame
self._fr = None
self._frstar = None
self._forcelist = forcelist
self._bodylist = bodies
self.explicit_kinematics = explicit_kinematics
self._constraint_solver = constraint_solver
self._initialize_vectors(q_ind, q_dependent, u_ind, u_dependent,
u_auxiliary)
_validate_coordinates(self.q, self.u)
self._initialize_kindiffeq_matrices(kd_eqs, kd_eqs_solver)
self._initialize_constraint_matrices(
configuration_constraints, velocity_constraints,
acceleration_constraints, constraint_solver)
def _initialize_vectors(self, q_ind, q_dep, u_ind, u_dep, u_aux):
"""Initialize the coordinate and speed vectors."""
none_handler = lambda x: Matrix(x) if x else Matrix()
# Initialize generalized coordinates
q_dep = none_handler(q_dep)
if not iterable(q_ind):
raise TypeError('Generalized coordinates must be an iterable.')
if not iterable(q_dep):
raise TypeError('Dependent coordinates must be an iterable.')
q_ind = Matrix(q_ind)
self._qdep = q_dep
self._q = Matrix([q_ind, q_dep])
self._qdot = self.q.diff(dynamicsymbols._t)
# Initialize generalized speeds
u_dep = none_handler(u_dep)
if not iterable(u_ind):
raise TypeError('Generalized speeds must be an iterable.')
if not iterable(u_dep):
raise TypeError('Dependent speeds must be an iterable.')
u_ind = Matrix(u_ind)
self._udep = u_dep
self._u = Matrix([u_ind, u_dep])
self._udot = self.u.diff(dynamicsymbols._t)
self._uaux = none_handler(u_aux)
def _initialize_constraint_matrices(self, config, vel, acc, linear_solver='LU'):
"""Initializes constraint matrices."""
linear_solver = _parse_linear_solver(linear_solver)
# Define vector dimensions
o = len(self.u)
m = len(self._udep)
p = o - m
none_handler = lambda x: Matrix(x) if x else Matrix()
# Initialize configuration constraints
config = none_handler(config)
if len(self._qdep) != len(config):
raise ValueError('There must be an equal number of dependent '
'coordinates and configuration constraints.')
self._f_h = none_handler(config)
# Initialize velocity and acceleration constraints
vel = none_handler(vel)
acc = none_handler(acc)
if len(vel) != m:
raise ValueError('There must be an equal number of dependent '
'speeds and velocity constraints.')
if acc and (len(acc) != m):
raise ValueError('There must be an equal number of dependent '
'speeds and acceleration constraints.')
if vel:
# When calling kanes_equations, another class instance will be
# created if auxiliary u's are present. In this case, the
# computation of kinetic differential equation matrices will be
# skipped as this was computed during the original KanesMethod
# object, and the qd_u_map will not be available.
if self._qdot_u_map is not None:
vel = msubs(vel, self._qdot_u_map)
self._k_nh, f_nh_neg = linear_eq_to_matrix(vel, self.u[:])
self._f_nh = -f_nh_neg
# If no acceleration constraints given, calculate them.
if not acc:
_f_dnh = (self._k_nh.diff(dynamicsymbols._t) * self.u +
self._f_nh.diff(dynamicsymbols._t))
if self._qdot_u_map is not None:
_f_dnh = msubs(_f_dnh, self._qdot_u_map)
self._f_dnh = _f_dnh
self._k_dnh = self._k_nh
else:
if self._qdot_u_map is not None:
acc = msubs(acc, self._qdot_u_map)
self._k_dnh, f_dnh_neg = linear_eq_to_matrix(acc, self._udot[:])
self._f_dnh = -f_dnh_neg
# Form of non-holonomic constraints is B*u + C = 0.
# We partition B into independent and dependent columns:
# Ars is then -B_dep.inv() * B_ind, and it relates dependent speeds
# to independent speeds as: udep = Ars*uind, neglecting the C term.
B_ind = self._k_nh[:, :p]
B_dep = self._k_nh[:, p:o]
self._Ars = -linear_solver(B_dep, B_ind)
else:
self._f_nh = Matrix()
self._k_nh = Matrix()
self._f_dnh = Matrix()
self._k_dnh = Matrix()
self._Ars = Matrix()
def _initialize_kindiffeq_matrices(self, kdeqs, linear_solver='LU'):
"""Initialize the kinematic differential equation matrices.
Parameters
==========
kdeqs : sequence of sympy expressions
Kinematic differential equations in the form of f(u,q',q,t) where
f() = 0. The equations have to be linear in the time-derivatives of
the generalized coordinates and in the generalized speeds.
"""
linear_solver = _parse_linear_solver(linear_solver)
if kdeqs:
if len(self.q) != len(kdeqs):
raise ValueError('There must be an equal number of kinematic '
'differential equations and coordinates.')
u = self.u
qdot = self._qdot
kdeqs = Matrix(kdeqs)
u_zero = dict.fromkeys(u, 0)
uaux_zero = dict.fromkeys(self._uaux, 0)
qdot_zero = dict.fromkeys(qdot, 0)
# Extract the linear coefficient matrices as per the following
# equation:
#
# k_ku(q,t)*u(t) + k_kqdot(q,t)*q'(t) + f_k(q,t) = 0
#
k_ku = kdeqs.jacobian(u)
k_kqdot = kdeqs.jacobian(qdot)
f_k = kdeqs.xreplace(u_zero).xreplace(qdot_zero)
# The kinematic differential equations should be linear in both q'
# and u so check for u and q' in the components.
dy_syms = find_dynamicsymbols(k_ku.row_join(k_kqdot).row_join(f_k))
nonlin_vars = [vari for vari in u[:] + qdot[:] if vari in dy_syms]
if nonlin_vars:
msg = ('The provided kinematic differential equations are '
'nonlinear in {}. They must be linear in the '
'generalized speeds and derivatives of the generalized '
'coordinates.')
raise ValueError(msg.format(nonlin_vars))
self._f_k_implicit = f_k.xreplace(uaux_zero)
self._k_ku_implicit = k_ku.xreplace(uaux_zero)
self._k_kqdot_implicit = k_kqdot
# Solve for q'(t) such that the coefficient matrices are now in
# this form:
#
# k_kqdot^-1*k_ku*u(t) + I*q'(t) + k_kqdot^-1*f_k = 0
#
# NOTE : Solving the kinematic differential equations here is not
# necessary and prevents the equations from being provided in fully
# implicit form.
f_k_explicit = linear_solver(k_kqdot, f_k)
k_ku_explicit = linear_solver(k_kqdot, k_ku)
self._qdot_u_map = dict(zip(qdot, -(k_ku_explicit*u + f_k_explicit)))
self._f_k = f_k_explicit.xreplace(uaux_zero)
self._k_ku = k_ku_explicit.xreplace(uaux_zero)
self._k_kqdot = eye(len(qdot))
else:
self._qdot_u_map = None
self._f_k_implicit = self._f_k = Matrix()
self._k_ku_implicit = self._k_ku = Matrix()
self._k_kqdot_implicit = self._k_kqdot = Matrix()
def _form_fr(self, fl):
"""Form the generalized active force."""
if fl is not None and (len(fl) == 0 or not iterable(fl)):
raise ValueError('Force pairs must be supplied in an '
'non-empty iterable or None.')
N = self._inertial
# pull out relevant velocities for constructing partial velocities
vel_list, f_list = _f_list_parser(fl, N)
vel_list = [msubs(i, self._qdot_u_map) for i in vel_list]
f_list = [msubs(i, self._qdot_u_map) for i in f_list]
# Fill Fr with dot product of partial velocities and forces
o = len(self.u)
b = len(f_list)
FR = zeros(o, 1)
partials = partial_velocity(vel_list, self.u, N)
for i in range(o):
FR[i] = sum(partials[j][i].dot(f_list[j]) for j in range(b))
# In case there are dependent speeds
if self._udep:
p = o - len(self._udep)
FRtilde = FR[:p, 0]
FRold = FR[p:o, 0]
FRtilde += self._Ars.T * FRold
FR = FRtilde
self._forcelist = fl
self._fr = FR
return FR
def _form_frstar(self, bl):
"""Form the generalized inertia force."""
if not iterable(bl):
raise TypeError('Bodies must be supplied in an iterable.')
t = dynamicsymbols._t
N = self._inertial
# Dicts setting things to zero
udot_zero = dict.fromkeys(self._udot, 0)
uaux_zero = dict.fromkeys(self._uaux, 0)
uauxdot = [diff(i, t) for i in self._uaux]
uauxdot_zero = dict.fromkeys(uauxdot, 0)
# Dictionary of q' and q'' to u and u'
q_ddot_u_map = {k.diff(t): v.diff(t).xreplace(
self._qdot_u_map) for (k, v) in self._qdot_u_map.items()}
q_ddot_u_map.update(self._qdot_u_map)
# Fill up the list of partials: format is a list with num elements
# equal to number of entries in body list. Each of these elements is a
# list - either of length 1 for the translational components of
# particles or of length 2 for the translational and rotational
# components of rigid bodies. The inner most list is the list of
# partial velocities.
def get_partial_velocity(body):
if isinstance(body, RigidBody):
vlist = [body.masscenter.vel(N), body.frame.ang_vel_in(N)]
elif isinstance(body, Particle):
vlist = [body.point.vel(N),]
else:
raise TypeError('The body list may only contain either '
'RigidBody or Particle as list elements.')
v = [msubs(vel, self._qdot_u_map) for vel in vlist]
return partial_velocity(v, self.u, N)
partials = [get_partial_velocity(body) for body in bl]
# Compute fr_star in two components:
# fr_star = -(MM*u' + nonMM)
o = len(self.u)
MM = zeros(o, o)
nonMM = zeros(o, 1)
zero_uaux = lambda expr: msubs(expr, uaux_zero)
zero_udot_uaux = lambda expr: msubs(msubs(expr, udot_zero), uaux_zero)
for i, body in enumerate(bl):
if isinstance(body, RigidBody):
M = zero_uaux(body.mass)
I = zero_uaux(body.central_inertia)
vel = zero_uaux(body.masscenter.vel(N))
omega = zero_uaux(body.frame.ang_vel_in(N))
acc = zero_udot_uaux(body.masscenter.acc(N))
inertial_force = (M.diff(t) * vel + M * acc)
inertial_torque = zero_uaux((I.dt(body.frame).dot(omega)) +
msubs(I.dot(body.frame.ang_acc_in(N)), udot_zero) +
(omega.cross(I.dot(omega))))
for j in range(o):
tmp_vel = zero_uaux(partials[i][0][j])
tmp_ang = zero_uaux(I.dot(partials[i][1][j]))
for k in range(o):
# translational
MM[j, k] += M*tmp_vel.dot(partials[i][0][k])
# rotational
MM[j, k] += tmp_ang.dot(partials[i][1][k])
nonMM[j] += inertial_force.dot(partials[i][0][j])
nonMM[j] += inertial_torque.dot(partials[i][1][j])
else:
M = zero_uaux(body.mass)
vel = zero_uaux(body.point.vel(N))
acc = zero_udot_uaux(body.point.acc(N))
inertial_force = (M.diff(t) * vel + M * acc)
for j in range(o):
temp = zero_uaux(partials[i][0][j])
for k in range(o):
MM[j, k] += M*temp.dot(partials[i][0][k])
nonMM[j] += inertial_force.dot(partials[i][0][j])
# Compose fr_star out of MM and nonMM
MM = zero_uaux(msubs(MM, q_ddot_u_map))
nonMM = msubs(msubs(nonMM, q_ddot_u_map),
udot_zero, uauxdot_zero, uaux_zero)
fr_star = -(MM * msubs(Matrix(self._udot), uauxdot_zero) + nonMM)
# If there are dependent speeds, we need to find fr_star_tilde
if self._udep:
p = o - len(self._udep)
fr_star_ind = fr_star[:p, 0]
fr_star_dep = fr_star[p:o, 0]
fr_star = fr_star_ind + (self._Ars.T * fr_star_dep)
# Apply the same to MM
MMi = MM[:p, :]
MMd = MM[p:o, :]
MM = MMi + (self._Ars.T * MMd)
# Apply the same to nonMM
nonMM = nonMM[:p, :] + (self._Ars.T * nonMM[p:o, :])
self._bodylist = bl
self._frstar = fr_star
self._k_d = MM
self._f_d = -(self._fr - nonMM)
return fr_star
def to_linearizer(self, linear_solver='LU'):
"""Returns an instance of the Linearizer class, initiated from the
data in the KanesMethod class. This may be more desirable than using
the linearize class method, as the Linearizer object will allow more
efficient recalculation (i.e. about varying operating points).
Parameters
==========
linear_solver : str, callable
Method used to solve the several symbolic linear systems of the
form ``A*x=b`` in the linearization process. If a string is
supplied, it should be a valid method that can be used with the
:meth:`sympy.matrices.matrixbase.MatrixBase.solve`. If a callable is
supplied, it should have the format ``x = f(A, b)``, where it
solves the equations and returns the solution. The default is
``'LU'`` which corresponds to SymPy's ``A.LUsolve(b)``.
``LUsolve()`` is fast to compute but will often result in
divide-by-zero and thus ``nan`` results.
Returns
=======
Linearizer
An instantiated
:class:`sympy.physics.mechanics.linearize.Linearizer`.
"""
if (self._fr is None) or (self._frstar is None):
raise ValueError('Need to compute Fr, Fr* first.')
# Get required equation components. The Kane's method class breaks
# these into pieces. Need to reassemble
f_c = self._f_h
if self._f_nh and self._k_nh:
f_v = self._f_nh + self._k_nh*Matrix(self.u)
else:
f_v = Matrix()
if self._f_dnh and self._k_dnh:
f_a = self._f_dnh + self._k_dnh*Matrix(self._udot)
else:
f_a = Matrix()
# Dicts to sub to zero, for splitting up expressions
u_zero = dict.fromkeys(self.u, 0)
ud_zero = dict.fromkeys(self._udot, 0)
qd_zero = dict.fromkeys(self._qdot, 0)
qd_u_zero = dict.fromkeys(Matrix([self._qdot, self.u]), 0)
# Break the kinematic differential eqs apart into f_0 and f_1
f_0 = msubs(self._f_k, u_zero) + self._k_kqdot*Matrix(self._qdot)
f_1 = msubs(self._f_k, qd_zero) + self._k_ku*Matrix(self.u)
# Break the dynamic differential eqs into f_2 and f_3
f_2 = msubs(self._frstar, qd_u_zero)
f_3 = msubs(self._frstar, ud_zero) + self._fr
f_4 = zeros(len(f_2), 1)
# Get the required vector components
q = self.q
u = self.u
if self._qdep:
q_i = q[:-len(self._qdep)]
else:
q_i = q
q_d = self._qdep
if self._udep:
u_i = u[:-len(self._udep)]
else:
u_i = u
u_d = self._udep
# Form dictionary to set auxiliary speeds & their derivatives to 0.
uaux = self._uaux
uauxdot = uaux.diff(dynamicsymbols._t)
uaux_zero = dict.fromkeys(Matrix([uaux, uauxdot]), 0)
# Checking for dynamic symbols outside the dynamic differential
# equations; throws error if there is.
sym_list = set(Matrix([q, self._qdot, u, self._udot, uaux, uauxdot]))
if any(find_dynamicsymbols(i, sym_list) for i in [self._k_kqdot,
self._k_ku, self._f_k, self._k_dnh, self._f_dnh, self._k_d]):
raise ValueError('Cannot have dynamicsymbols outside dynamic \
forcing vector.')
# Find all other dynamic symbols, forming the forcing vector r.
# Sort r to make it canonical.
r = list(find_dynamicsymbols(msubs(self._f_d, uaux_zero), sym_list))
r.sort(key=default_sort_key)
# Check for any derivatives of variables in r that are also found in r.
for i in r:
if diff(i, dynamicsymbols._t) in r:
raise ValueError('Cannot have derivatives of specified \
quantities when linearizing forcing terms.')
return Linearizer(f_0, f_1, f_2, f_3, f_4, f_c, f_v, f_a, q, u, q_i,
q_d, u_i, u_d, r, linear_solver=linear_solver)
# TODO : Remove `new_method` after 1.1 has been released.
def linearize(self, *, new_method=None, linear_solver='LU', **kwargs):
""" Linearize the equations of motion about a symbolic operating point.
Parameters
==========
new_method
Deprecated, does nothing and will be removed.
linear_solver : str, callable
Method used to solve the several symbolic linear systems of the
form ``A*x=b`` in the linearization process. If a string is
supplied, it should be a valid method that can be used with the
:meth:`sympy.matrices.matrixbase.MatrixBase.solve`. If a callable is
supplied, it should have the format ``x = f(A, b)``, where it
solves the equations and returns the solution. The default is
``'LU'`` which corresponds to SymPy's ``A.LUsolve(b)``.
``LUsolve()`` is fast to compute but will often result in
divide-by-zero and thus ``nan`` results.
**kwargs
Extra keyword arguments are passed to
:meth:`sympy.physics.mechanics.linearize.Linearizer.linearize`.
Explanation
===========
If kwarg A_and_B is False (default), returns M, A, B, r for the
linearized form, M*[q', u']^T = A*[q_ind, u_ind]^T + B*r.
If kwarg A_and_B is True, returns A, B, r for the linearized form
dx = A*x + B*r, where x = [q_ind, u_ind]^T. Note that this is
computationally intensive if there are many symbolic parameters. For
this reason, it may be more desirable to use the default A_and_B=False,
returning M, A, and B. Values may then be substituted in to these
matrices, and the state space form found as
A = P.T*M.inv()*A, B = P.T*M.inv()*B, where P = Linearizer.perm_mat.
In both cases, r is found as all dynamicsymbols in the equations of
motion that are not part of q, u, q', or u'. They are sorted in
canonical form.
The operating points may be also entered using the ``op_point`` kwarg.
This takes a dictionary of {symbol: value}, or a an iterable of such
dictionaries. The values may be numeric or symbolic. The more values
you can specify beforehand, the faster this computation will run.
For more documentation, please see the ``Linearizer`` class.
"""
linearizer = self.to_linearizer(linear_solver=linear_solver)
result = linearizer.linearize(**kwargs)
return result + (linearizer.r,)
def kanes_equations(self, bodies=None, loads=None):
""" Method to form Kane's equations, Fr + Fr* = 0.
Explanation
===========
Returns (Fr, Fr*). In the case where auxiliary generalized speeds are
present (say, s auxiliary speeds, o generalized speeds, and m motion
constraints) the length of the returned vectors will be o - m + s in
length. The first o - m equations will be the constrained Kane's
equations, then the s auxiliary Kane's equations. These auxiliary
equations can be accessed with the auxiliary_eqs property.
Parameters
==========
bodies : iterable
An iterable of all RigidBody's and Particle's in the system.
A system must have at least one body.
loads : iterable
Takes in an iterable of (Particle, Vector) or (ReferenceFrame, Vector)
tuples which represent the force at a point or torque on a frame.
Must be either a non-empty iterable of tuples or None which corresponds
to a system with no constraints.
"""
if bodies is None:
bodies = self.bodies
if loads is None and self._forcelist is not None:
loads = self._forcelist
if loads == []:
loads = None
if not self._k_kqdot:
raise AttributeError('Create an instance of KanesMethod with '
'kinematic differential equations to use this method.')
fr = self._form_fr(loads)
frstar = self._form_frstar(bodies)
if self._uaux:
if not self._udep:
km = KanesMethod(self._inertial, self.q, self._uaux,
u_auxiliary=self._uaux, constraint_solver=self._constraint_solver)
else:
km = KanesMethod(self._inertial, self.q, self._uaux,
u_auxiliary=self._uaux, u_dependent=self._udep,
velocity_constraints=(self._k_nh * self.u +
self._f_nh),
acceleration_constraints=(self._k_dnh * self._udot +
self._f_dnh),
constraint_solver=self._constraint_solver
)
km._qdot_u_map = self._qdot_u_map
self._km = km
fraux = km._form_fr(loads)
frstaraux = km._form_frstar(bodies)
self._aux_eq = fraux + frstaraux
self._fr = fr.col_join(fraux)
self._frstar = frstar.col_join(frstaraux)
return (self._fr, self._frstar)
def _form_eoms(self):
fr, frstar = self.kanes_equations(self.bodylist, self.forcelist)
return fr + frstar
def rhs(self, inv_method=None):
"""Returns the system's equations of motion in first order form. The
output is the right hand side of::
x' = |q'| =: f(q, u, r, p, t)
|u'|
The right hand side is what is needed by most numerical ODE
integrators.
Parameters
==========
inv_method : str
The specific sympy inverse matrix calculation method to use. For a
list of valid methods, see
:meth:`~sympy.matrices.matrixbase.MatrixBase.inv`
"""
rhs = zeros(len(self.q) + len(self.u), 1)
kdes = self.kindiffdict()
for i, q_i in enumerate(self.q):
rhs[i] = kdes[q_i.diff()]
if inv_method is None:
rhs[len(self.q):, 0] = self.mass_matrix.LUsolve(self.forcing)
else:
rhs[len(self.q):, 0] = (self.mass_matrix.inv(inv_method,
try_block_diag=True) *
self.forcing)
return rhs
def kindiffdict(self):
"""Returns a dictionary mapping q' to u."""
if not self._qdot_u_map:
raise AttributeError('Create an instance of KanesMethod with '
'kinematic differential equations to use this method.')
return self._qdot_u_map
@property
def auxiliary_eqs(self):
"""A matrix containing the auxiliary equations."""
if not self._fr or not self._frstar:
raise ValueError('Need to compute Fr, Fr* first.')
if not self._uaux:
raise ValueError('No auxiliary speeds have been declared.')
return self._aux_eq
@property
def mass_matrix_kin(self):
r"""The kinematic "mass matrix" $\mathbf{k_{k\dot{q}}}$ of the system."""
return self._k_kqdot if self.explicit_kinematics else self._k_kqdot_implicit
@property
def forcing_kin(self):
"""The kinematic "forcing vector" of the system."""
if self.explicit_kinematics:
return -(self._k_ku * Matrix(self.u) + self._f_k)
else:
return -(self._k_ku_implicit * Matrix(self.u) + self._f_k_implicit)
@property
def mass_matrix(self):
"""The mass matrix of the system."""
if not self._fr or not self._frstar:
raise ValueError('Need to compute Fr, Fr* first.')
return Matrix([self._k_d, self._k_dnh])
@property
def forcing(self):
"""The forcing vector of the system."""
if not self._fr or not self._frstar:
raise ValueError('Need to compute Fr, Fr* first.')
return -Matrix([self._f_d, self._f_dnh])
@property
def mass_matrix_full(self):
"""The mass matrix of the system, augmented by the kinematic
differential equations in explicit or implicit form."""
if not self._fr or not self._frstar:
raise ValueError('Need to compute Fr, Fr* first.')
o, n = len(self.u), len(self.q)
return (self.mass_matrix_kin.row_join(zeros(n, o))).col_join(
zeros(o, n).row_join(self.mass_matrix))
@property
def forcing_full(self):
"""The forcing vector of the system, augmented by the kinematic
differential equations in explicit or implicit form."""
return Matrix([self.forcing_kin, self.forcing])
@property
def q(self):
return self._q
@property
def u(self):
return self._u
@property
def bodylist(self):
return self._bodylist
@property
def forcelist(self):
return self._forcelist
@property
def bodies(self):
return self._bodylist
@property
def loads(self):
return self._forcelist

Some files were not shown because too many files have changed in this diff Show More