{ "cells": [ { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ ".. _user_guide_linalg:" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# Linear Algebra" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "📖 :ref:`Linear Algebra API Reference `" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ ".. _Awkward: https://awkward-array.org/doc/main/\n", ".. _NumPy: https://numpy.org/\n", ".. _Numba: https://numba.pydata.org/" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "It is safe to say that linear algebra is the cornerstone of modern mathematical physics, providing the essential tools for understanding and solving a wide range of problems. From the study of vector spaces and transformations to eigenvalue problems and matrix operations, linear algebra forms the foundation for many physical theories. Whether in classical mechanics, quantum physics, relativity, or electrodynamics, linear algebra enables the representation of physical quantities, the formulation of physical laws, and the manipulation of complex systems in an efficient and structured way. Its concepts permeate the formulation of everything from equations of motion to the quantum states of particles, making it indispensable for both theoretical insights and practical computations." ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Opposed to the widespread disbelief, vectors are not merely 1D arrays, and tensors are not simply multi-dimensional arrays either. While it’s common to represent vectors and tensors in this form for computational convenience, their true nature is far more profound. Vectors exist in an abstract vector space and only take on specific numerical values when projected onto a basis. Similarly, tensors are generalized multi-linear mappings that relate vectors, scalars, and other tensors in ways that transcend their array-like representations. Their intrinsic properties remain invariant under transformations, making them essential tools in describing physical phenomena, from forces and fields to stress and strain in materials, in a way that doesn’t depend on arbitrary choices of coordinates." ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The linear algebra submodule in `sigmaepsilon.math` provides **implementations of Vectors and Tensors as they are fundamentally understood in physics.**" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ ".. note::\n", "\n", " Some of the cell outputs on this page don't render properly in dark mode." ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Arrays" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "**Arrays** are ordered collections of elements, typically numbers, arranged in one or more dimensions. In mathematics and computing, they represent data in a structured form, making it easier to manipulate large datasets. A **one-dimensional array** is a simple list (similar to a vector), while **multi-dimensional arrays** can represent more complex structures like matrices or tensors (but they are not tensors). Arrays are fundamental in fields like data science, machine learning, and numerical computation, where they are used to store and process large volumes of data efficiently." ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The library provides some classes for arrays, namely :class:`~sigmaepsilon.math.linalg.meta.Array`, :class:`~sigmaepsilon.math.linalg.sparse.jaggedarray.JaggedArray` and :class:`~sigmaepsilon.math.linalg.sparse.csr.csr_matrix`." ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### Array" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The :class:`~sigmaepsilon.math.linalg.meta.Array` class is basically the same as the `ndarray` class in `NumPy`_, with the minor difference that it has a safe metaclass, meaning that there is a safety mechanism that prevents you from unintentionally crashing the internal behaviour of the class upon subclassing. This practically means that you will see an error if you try to shadow a definition in any of the base classes of the class. For this reason, it is safer to subclass this class rather than to directly subclass NumPy's `ndarray` class." ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### JaggedArray" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Jagged arrays are like matrices (2d arrays) that have rows of different length. The :class:`~sigmaepsilon.math.linalg.sparse.jaggedarray.JaggedArray` class wraps either a `NumPy`_ array or an `Awkward`_ array in the background, depending on the input and generalizes some matrix operations for the jagged case." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([[1, 2, 3],\n", " [4, 5, 6],\n", " [7, 8, 9]])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sigmaepsilon.math.linalg import JaggedArray\n", "import numpy as np\n", "\n", "data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])\n", "arr = JaggedArray(data, cuts=[3, 3, 3])\n", "arr" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(arr._wrapped)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr.is_jagged()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([3, 3, 3], dtype=int64)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr.widths()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The underlying array in this case is a NumPy array and it should behave like one, you can use it as a drop-in replacement. However, the constructor of the class accepts a `force_numpy` parameter, which is `True` by default. It means that the constructor will try to use NumPy whenever possible and only uses Awkward if the input is truly jagged." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
[[1, 2, 3],\n",
       " [4, 5, 6],\n",
       " [7, 8, 9, 10]]\n",
       "---------------------\n",
       "type: 3 * var * int32
" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])\n", "arr = JaggedArray(data, cuts=[3, 3, 4])\n", "arr" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "awkward.highlevel.Array" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(arr._wrapped)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr.is_jagged()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([3, 3, 4], dtype=int64)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr.widths()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The class generalizes some NumPy functions to jagged arrays, like `unique` and `concatenate`." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 2, 3, 4, 5, 6, 7, 8, 9])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 9])\n", "arr = JaggedArray(data, cuts=[3, 3, 4])\n", "np.unique(arr)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
[[1, 2, 3],\n",
       " [4, 5, 6],\n",
       " [7, 8, 9, 10],\n",
       " [1, 2, 3],\n",
       " [4, 5],\n",
       " [6, 7, 8, 9, 10, 11]]\n",
       "----------------------\n",
       "type: 6 * var * int32
" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])\n", "arr1 = JaggedArray(data, cuts=[3, 3, 4])\n", "\n", "data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])\n", "arr2 = JaggedArray(data, cuts=[3, 2, 6])\n", "\n", "np.concatenate([arr1, arr2])" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "For other operations, use the `to_array` method of the instance, which will return either a NumPy or an Awkward array, both of which has support for all the operations you might need, and you can use the result to initialize another :class:`~sigmaepsilon.math.linalg.sparse.jaggedarray.JaggedArray` instance if necessary." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
[[3, 6, 9],\n",
       " [12, 15, 18],\n",
       " [21, 24, 27, 27]]\n",
       "---------------------\n",
       "type: 3 * var * int32
" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 9])\n", "arr = JaggedArray(data, cuts=[3, 3, 4])\n", "arr = JaggedArray(arr.to_array()*3)\n", "arr" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "You can also use the :class:`~sigmaepsilon.math.linalg.sparse.jaggedarray.JaggedArray` class to transform between different data formats. For the details of this transformation functions, consult the API Reference." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr.to_array() # returns a NumPy or an Awkward array\n", "arr.to_ak() # returns an Awkward array\n", "arr.to_numpy() # returns a dense NumPy array\n", "arr.to_list() # returns a list of lists\n", "arr.to_csr() # returns an instance of sigmaepsilon.math.linalg.sparse.csr_matrix\n", "arr.to_scipy() # returns an instance of scipy.sparse.csr_matrix" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### csr_matrix" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The :class:`~sigmaepsilon.math.linalg.sparse.csr.csr_matrix` class is basically a remake of a similar class in SciPy, with the addition that it is *Numba-jittable*. To create an instance, you can start from a :class:`~sigmaepsilon.math.linalg.sparse.jaggedarray.JaggedArray` instance, or use any other input you would use with a SciPy csr matrix." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "3x4 CSR matrix of 12 values." ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.sparse import csr_matrix as csr_scipy\n", "from sigmaepsilon.math.linalg import csr_matrix\n", "\n", "# create a JaggedArray and convert it to a scipy.sparse.csr_matrix\n", "data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])\n", "scipy_matrix = JaggedArray(data, cuts=[3, 3, 4]).to_csr()\n", "\n", "# create a csr_matrix instance from a dense NumPy array\n", "scipy_matrix = csr_scipy((3, 4), dtype=np.int8).toarray()\n", "csr_matrix(scipy_matrix)" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Instances of the :class:`~sigmaepsilon.math.linalg.sparse.csr.csr_matrix` class are Numba-jittable, with some attributes and methods accessible inside Numba-jitted functions, even in nopython-mode." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([1., 2.]), array([0, 2]))" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from numba import jit\n", "\n", "\n", "@jit(nopython=True, nogil=True)\n", "def numba_nopython(csr: csr_matrix, i: int):\n", " csr.data\n", " csr.indices\n", " csr.indptr\n", " csr.shape\n", " return csr.row(i), csr.irow(i)\n", "\n", "\n", "row = np.array([0, 0, 1, 2, 2, 2])\n", "col = np.array([0, 2, 2, 0, 1, 2])\n", "data = np.array([1, 2, 3, 4, 5, 6])\n", "matrix = csr_scipy((data, (row, col)), shape=(3, 3))\n", "\n", "csr = csr_matrix(matrix)\n", "numba_nopython(csr, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reference Frames" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Whether you are working with vectors or tensors, in order to represent them as one or multi-dimensional arrays, you need a reference frame. The :ref:`sigmaepsilon.math.linalg ` module provides you with the following classes to represent reference frames:\n", "\n", "- :class:`~sigmaepsilon.math.linalg.frame.ReferenceFrame`: The most general reference frame, allowing for arbitrary settings. Base vectors do not need to be orthogonal, nor normalized.\n", "- :class:`~sigmaepsilon.math.linalg.frame.RectangularFrame`: For the case when base vectors are not normalized, but orthogonal.\n", "- :class:`~sigmaepsilon.math.linalg.frame.CartesianFrame`: For frames with orthonormal base vectors.\n", "\n", "Of course, you can just use the most general :class:`~sigmaepsilon.math.linalg.frame.ReferenceFrame` class and accept that some calculations might involve computing the Gram matrix, which in some cases would simplify to the identity matrix, meaning that the class is a bit overpower for your needs. Or, you can go with the :class:`~sigmaepsilon.math.linalg.frame.RectangularFrame` and :class:`~sigmaepsilon.math.linalg.frame.CartesianFrame` classes that might be a better choice for your problem and they are computationally more efficient. The good thing is that you don't have to worry about crossing the boundaries between these types, the library takes care of everything for you." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "from sigmaepsilon.math.linalg import ReferenceFrame, RectangularFrame, CartesianFrame" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's say you get your hands on a cartesian frame (one with orthonormal basis vectors)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "frame = CartesianFrame(dim=3)" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Now there are some operations you can apply on a frame. Some of these operations would result in another frame, which is not cartesian. Some operations would modify the frame inplace and result in the instance losing the property of orthogonality or normality. For instance, if all basis vectors were multiplied by a scalar, the new basis vectors would be still orthogonal, but not normal." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "sigmaepsilon.math.linalg.frame.RectangularFrame" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(frame*2)" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "As you can see, in this case the instance created by the operation is the next most general one, which is a :class:`~sigmaepsilon.math.linalg.frame.RectangularFrame`. If let's say no you wanted to increment all coordinates of the basis vectors by the same scalar, there would be no guarantee of the basis vectors of the resulting frame to being either orthogonal nor normal, hence the frame created by the operation is the most general one, an instance of `ReferenceFrame`." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "sigmaepsilon.math.linalg.frame.ReferenceFrame" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(frame+2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " The good new is that you don't really have to worry about these things as the library will choose the most appropriate class for the output, depending on the nature of the operation." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "sigmaepsilon.math.linalg.frame.CartesianFrame" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(frame)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "sigmaepsilon.math.linalg.frame.RectangularFrame" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame *=2\n", "type(frame)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "sigmaepsilon.math.linalg.frame.ReferenceFrame" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame +=2\n", "type(frame)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "Array([[4., 2., 2.],\n", " [2., 4., 2.],\n", " [2., 2., 4.]])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.axes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can see that the basis vectors of the frame are no longer orthogonal or normed, but they are still independent." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.is_cartesian" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.is_rectangular" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.is_independent" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "It is also possible to get the dual frame (aka. reciprocal frame) of a frame. This is one example where the :class:`~sigmaepsilon.math.linalg.frame.CartesianFrame` is computationally more effective, as the dual frame of an orthonormal frame is itself." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "Array([[ 0.375, -0.125, -0.125],\n", " [-0.125, 0.375, -0.125],\n", " [-0.125, -0.125, 0.375]])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.dual()" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Vectors and Tensors" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "**Vectors** and **tensors** are fundamental concepts in mathematics and physics, representing quantities that can describe physical phenomena. A **vector** is a one-dimensional object characterized by both magnitude and direction, commonly used to represent physical quantities like velocity, force, and displacement. A **tensor** is a more generalized mathematical object that can describe relationships between vectors or other tensors and extends across multiple dimensions. Tensors are crucial in fields like continuum mechanics, electromagnetism, and general relativity, where they model complex interactions in various coordinate systems." ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### Vectors" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "As a basic example, consider two orthonormal frames in Euclidean space. The first one is the standard frame defined with the unit matrix $\\mathbf{I}$, the second is obtained by rotating $\\mathbf{I}$ around the $z$ axis with $90$ degrees. A vector is defined in $\\mathbf{I}$ (the source) and we want to know it's components in the rotated frame (the target)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We begin by creating the coordinates of the vector as an 1d array." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([0.8660254, 0.5 , 0. ])" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "\n", "# the coordinates of a vector in the source frame\n", "arr_source = np.array([3 ** 0.5 / 2, 0.5, 0])\n", "arr_source" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "In sigmaepsilon.math, the :class:`~sigmaepsilon.math.linalg.frame.ReferenceFrame` class provides the machinery to establish a connection between two frames. First, we create the frames as:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "from sigmaepsilon.math.linalg import ReferenceFrame\n", "\n", "source_frame = ReferenceFrame(dim=3) # this is a 3 by 3 identity matrix\n", "target_frame = source_frame.orient_new('Body', [0, 0, 90*np.pi/180], 'XYZ')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To transform the vector into the target frame, we can use the Vector object directly, and it handles everything in the background." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([ 0.5 , -0.8660254, 0. ])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sigmaepsilon.math.linalg import Vector\n", "\n", "arr_target = Vector(arr_source, frame=source_frame).show(target_frame)\n", "arr_target" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### Tensors" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The library also provides tools to deal with tensors through the :class:`~sigmaepsilon.math.linalg.tensor.Tensor` class and its derivative classes. The implementation has only three limitations:\n", "\n", "* The reference frames must be linear.\n", "* All the indices of a tensor must be either covariant or contravariant.\n", "* The sum total of indices used in any kind of operation must be less than or equal to 26. This is because the implementation relies on the `einsum` function in NumPy, and the english alphapet has 26 characters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a Tensor of order 6 in a frame with random components:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "from sigmaepsilon.math.linalg import Tensor, ReferenceFrame\n", "from numpy.random import rand\n", "\n", "frame = ReferenceFrame(dim=3)\n", "array = rand(3, 3, 3, 3, 3, 3)\n", "A = Tensor(array, frame=frame)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get the tensor in the dual frame:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "A_dual = A.dual()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create an other tensor, in this case a 5th-order one, and calculate their generalized dot product, which is a 9th-order tensor:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "from sigmaepsilon.math.linalg import dot\n", "\n", "array = rand(3, 3, 3, 3, 3)\n", "B = Tensor(array, frame=frame)\n", "C = dot(A, B, axes=[0, 0])\n", "assert C.rank == (A.rank + B.rank - 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The library also provides dedicated tensor classes for the most common quantities of engineering practice." ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "👉 :class:`~sigmaepsilon.math.linalg.tensor.Tensor2`\n", ":class:`~sigmaepsilon.math.linalg.tensor.Tensor2x3`\n", ":class:`~sigmaepsilon.math.linalg.tensor.Tensor4`\n", ":class:`~sigmaepsilon.math.linalg.tensor.Tensor4x3`" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### Objectivity of Tensorial Quantities" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Tensors are mathematical objects that generalize scalars, vectors, and higher-order quantities, and they are defined by how they transform between different coordinate systems. This property makes tensors essential in physics, engineering, and mathematics, where the description of physical quantities must remain consistent regardless of the chosen frame of reference. A remarkable feature of the library is that tensorial quantities keep their objectivity with respect to changes of their host reference frame. To show this, create a new vector." ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "arr = np.array([1, 0, 0])\n", "frame = ReferenceFrame(dim=3)\n", "vector = Vector(arr, frame=frame)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now transform the hosting frame of the vector inplace." ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "Array([[ 6.123234e-17, -1.000000e+00, 0.000000e+00],\n", " [ 1.000000e+00, 6.123234e-17, 0.000000e+00],\n", " [ 0.000000e+00, 0.000000e+00, 1.000000e+00]])" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.orient('Body', [0, 0, -90*np.pi/180], 'XYZ')\n", "frame.axes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you now check the coordinates of the vector, you will see that they have changed, due to the changes in the hosting frame." ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "Array([6.123234e-17, 1.000000e+00, 0.000000e+00])" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vector.array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, the vector stayed the same. You can check this, by asking for the coordinates of the vector in the ambient frame (which in this case is the frame it was defined in)." ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([1., 0., 0.])" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vector.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The point is that if a quantity is tensorial, it doesn't matter which frame you define the quantity in, since its meaning is independent from the hosting frame. Therefore, if the hosting frame changes, the coordinates of the tensorial quantity must change as well. A more trivial example is the following:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "Array([0.5, 0. , 0. ])" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr = np.array([1, 0, 0])\n", "frame = ReferenceFrame(dim=3)\n", "vector = Vector(arr, frame=frame)\n", "frame *= 2\n", "vector.array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It makes perfect sense. The length of the basis vectors doubled, so the coordinates shrunk to half of the original values. Resetting the basis vectors (axes) of the hosting frame, resets the coordinates of the vector as well." ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "Array([1., 0., 0.])" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.axes = np.eye(3)\n", "vector.array" ] }, { "cell_type": "raw", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Thanks to the flexibility of the :class:`~sigmaepsilon.math.linalg.frame.ReferenceFrame` class, you can assign any 3x3 positive definite matrix to the axes of the frame, and the meaning of the vector wouldn't change a thing. You can check this by asking for the coordinates of the vector in a fixed frame, which can be the ambient frame as well." ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([ 1.00000000e+00, -2.25984087e-14, -2.22766958e-14])" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sigmaepsilon.math.linalg import random_posdef_matrix\n", "\n", "random_axes = random_posdef_matrix(3) * 100\n", "frame.axes = random_axes\n", "vector.show()" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "(False, False)" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame.is_rectangular, frame.is_cartesian" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "It is clear that the coordinates of the vector are still the same in the ambient frame as they were at the beginning and the next cell shows that if we take the dot product of the vector with itself, the result is still $1.0$, even though the coordinates of it have changed and the frame of it is a random one." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "0.9999999999999374" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sigmaepsilon.math.linalg import dot\n", "\n", "dot(vector, vector)" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Again, this demonstrates that if a quantity is tensorial, the specific frame in which it is represented as an array is less important than ensuring it has a frame. This flexibility is valuable because, in many practical problems, it's easier to define a tensor's components in one particular frame. By having the freedom to choose the most advantageous frame, we can simplify the process of solving the problem while maintaining the tensor's properties. The next example shows that if we define two vectors, and we repeatedly alter their frames, their dot product doesn't change." ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Coordinates of vector A: [1. 0. 0.]\n", "Coordinates of vector B: [1. 0. 0.]\n", "dot(vectorA, vectorB): 1.0\n", "dot(arrayA, arrayB): 1.0\n", "---------------------------------------------------------------\n", "Coordinates of vector A: [ 0.14003901 0.11117086 -0.2200233 ]\n", "Coordinates of vector B: [ 0.04139533 -0.00550249 -0.03643786]\n", "dot(vectorA, vectorB): 1.000000000000006\n", "dot(arrayA, arrayB): 0.013202422509482573\n", "Coordinates of vector A: [ 4.00354869 -4.19724454 -0.14587514]\n", "Coordinates of vector B: [ 2.64746757 0.48450066 -1.39766909]\n", "dot(vectorA, vectorB): 1.0000000000021052\n", "dot(arrayA, arrayB): 8.769582746424097\n", "Coordinates of vector A: [ 2.06055094 -1.64386404 -1.55690927]\n", "Coordinates of vector B: [ 3.55832797 0.04023036 -3.09062927]\n", "dot(vectorA, vectorB): 0.9999999999802368\n", "dot(arrayA, arrayB): 12.077812168548878\n" ] } ], "source": [ "arr = np.array([1, 0, 0])\n", "frameA = ReferenceFrame(dim=3)\n", "vectorA = Vector(arr, frame=frameA)\n", "frameB = ReferenceFrame(dim=3)\n", "vectorB = Vector(arr, frame=frameB)\n", "\n", "print(\"Coordinates of vector A: \", vectorA.array)\n", "print(\"Coordinates of vector B: \", vectorB.array)\n", "print(\"dot(vectorA, vectorB): \", dot(vectorA, vectorB))\n", "print(\"dot(arrayA, arrayB): \", dot(vectorA.array, vectorB.array))\n", "print(63*\"-\")\n", "\n", "for _ in range(3):\n", " random_axes = random_posdef_matrix(3) * 100\n", " frameA.axes = random_axes\n", " random_axes = random_posdef_matrix(3) * 100\n", " frameB.axes = random_axes\n", " print(\"Coordinates of vector A: \", vectorA.array)\n", " print(\"Coordinates of vector B: \", vectorB.array)\n", " print(\"dot(vectorA, vectorB): \", dot(vectorA, vectorB))\n", " print(\"dot(arrayA, arrayB): \", dot(vectorA.array, vectorB.array))" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "This explains why it is incorrect to say that vectors are 1d arrays. They are only representations of vectors in the space of some basis vectors, but they are not the vectors themselves. You can also observe how the dot product of the two vectors remains the same, but the dot product of their respective arrays changes." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.11" }, "vscode": { "interpreter": { "hash": "4e251a336b180e3c877fd4b81be72acfad98293ac2abcf90f00390a06765d313" } } }, "nbformat": 4, "nbformat_minor": 4 }