Vector operators in SymPy



In [1]:

    
from sympy import *



In [2]:

    
init_session()









    



IPython console for SymPy 1.4 (Python 3.6.9-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/1.4/



In [3]:

    
def scale_coeff(r_vec, coords):
    """
    Compyte scale coefficients for the vector
    tranform given by r_vec.
    
    Parameters
    =======
    r_vec : Matrix (3, 1)
        Transform vector (x, y, z) as a functoin of coordinates
        u1, u2, u3.
    coords : Tupl (3)
        Coordinates for the new reference system.
        
    Returns
    ========
    h_vec : Tuple (3)
        Scale coefficients.
    """
    if type(r_vec) == list:
        r_vec = Matrix(r_vec)
    u1, u2, u3 = coords
    h1 = simplify((r_vec.diff(u1)).norm())
    h2 = simplify((r_vec.diff(u2)).norm())
    h3 = simplify((r_vec.diff(u3)).norm())
    return h1, h2, h3

As an example, the tranformation for Spherical coordinates.



In [4]:

    
r, theta, phi = symbols("r theta phi", positive=True)
h_vec = scale_coeff([r*sin(theta)*cos(phi),
                     r*sin(theta)*sin(phi),
                     r*cos(theta)],
                    [r, theta, phi])



In [5]:

    
h_vec









    Out[5]:





$\displaystyle \left( 1, \  r, \  r \left|{\sin{\left(\theta \right)}}\right|\right)$

And, since $\theta$ is in $[0, \pi]$, the last component is $r \sin(\theta)$.



In [6]:

    
def grad(u, coords=(x, y, z), h_vec=(1, 1, 1)):
    """
    Compute the gradient of a scalara function phi.

    Parameters
    ==========
    u : SymPy expression
        Scalar function to compute the gradient from.
    coords : Tuple (3), optional
        Coordinates for the new reference system. This is an optional
        parameters, and it takes a cartesian (x, y, z), as default.
    h_vec : Tuple (3), optional
        Scale coefficients for the new coordinate system. It takes
        (1, 1, 1), as default.
        
    Returns
    ========
    gradient: Matrix (3, 1)
        Column vector with the components of the gradient.
    """
    return Matrix(3, 1, lambda i, j: u.diff(coords[i])/h_vec[j])

Let's see an example.



In [7]:

    
grad(-(cos(x)**2 + cos(y)**2)**2)









    Out[7]:





$\displaystyle \left[\begin{matrix}4 \left(\cos^{2}{\left(x \right)} + \cos^{2}{\left(y \right)}\right) \sin{\left(x \right)} \cos{\left(x \right)}\\4 \left(\cos^{2}{\left(x \right)} + \cos^{2}{\left(y \right)}\right) \sin{\left(y \right)} \cos{\left(y \right)}\\0\end{matrix}\right]$



In [8]:

    
def grad_vec(A, coords=(x, y, z), h_vec=(1, 1, 1)):
    """
    Gradient of a vector function A.
    
    Parameters
    ==========
    A : Matrix (3, 1), list
        Vector function to compute the gradient from.
    coords : Tuple (3), optional
        Coordinates for the new reference system. This is an optional parameter
        it takes (x, y, z) as default.
    h_vec : Tuple (3), optional
        Scale coefficients for the new coordinate system. It takes
        (1, 1, 1), as default.
        
    Returns
    =======
    gradient: Matrix (3, 3)
        Matrix with the components of the gradient. The position (i, j) has as components
        diff(A[i], coords[j].
    """ 
    return Matrix(3, 3, lambda i, j: (S(1)/h_vec[j])*A[i].diff(coords[j]))

Let's see a simple example.



In [9]:

    
grad_vec([x*y*z, x*y*z, x*y*z])









    Out[9]:





$\displaystyle \left[\begin{matrix}y z & x z & x y\\y z & x z & x y\\y z & x z & x y\end{matrix}\right]$



In [10]:

    
def div(A, coords=(x, y, z), h_vec=(1, 1, 1)):
    """
    Divergence of the vector function A.
    
    Parameters
    ==========
    A : Matrix, list
        Scalar function to compute the divergence from.
    coords : Tuple (3), optional
        Coordinates for the new reference system. This is an optional parameter
        it takes (x, y, z) as default.
    h_vec : Tuple (3), optional
        Scale coefficients for the new coordinate system. It takes
        (1, 1, 1), as default.
        
    Returns
    ========
    divergence: SymPy expression
        Divergence of A.
    """  
    h = h_vec[0]*h_vec[1]*h_vec[2]
    aux = simplify((S(1)/h)*sum(diff(A[k]*h/h_vec[k], coords[k])
                                for k in range(3)))
    return aux

Let's see a simple example.



In [11]:

    
div([x**2 + y*z, y**2+x*z, z**2 + x*y])









    Out[11]:





$\displaystyle 2 x + 2 y + 2 z$



In [12]:

    
def curl(A, coords=(x, y, z), h_vec=(1, 1, 1)):
    """
    Curl of a function vector A.
    
    Parameters
    ==========
    A : Matrix, List
        Vector function to compute the curl from.
    coords : Tuple (3), optional
        Coordinates for the new reference system. This is an optional parameter
        it takes (x, y, z) as default.
    h_vec : Tuple (3), optional
        Scale coefficients for the new coordinate system. It takes
        (1, 1, 1), as default.

    Returns
    =======
    curl : Matrix (3, 1)
        Column vector with the curl of A.
    """  
    perm = lambda i, j, k: (i - j)*(j - k)*(k - i)/S(2)
    h = h_vec[0]*h_vec[1]*h_vec[2]
    aux = [(S(1)/h)*sum(perm(i, j, k)*h_vec[i]*diff(A[k]*h_vec[k], coords[j]) for j in range(3) for k in range(3))
           for i in range(3)]
    return Matrix(aux)

Let's see an example.



In [13]:

    
curl([0, -x**2, 0])









    Out[13]:





$\displaystyle \left[\begin{matrix}0\\0\\- 2 x\end{matrix}\right]$



In [14]:

    
def lap(u, coords=(x, y, z), h_vec=(1, 1, 1)):
    """
    Laplacian of the scalar function phi.
    
    Parameters
    ==========
    u : SymPy expression
        Scalar function to compute the gradient from.
    coords : Tuple (3), optional
        Coordinates for the new reference system. This is an optional
        parameters, and it takes a cartesian (x, y, z), as default.
    h_vec : Tuple (3), optional
        Scale coefficients for the new coordinate system. It takes
        (1, 1, 1), as default.
        
    Returns
    =======
    laplacian: Sympy expression
        Laplacian of phi.
    """
    h = h_vec[0]*h_vec[1]*h_vec[2]
    return sum([1/h*diff(h/h_vec[k]**2*u.diff(coords[k]), coords[k]) for k in range(3)])

Hagamos un ejemplo simple



In [15]:

    
lap(x**2 +  y**2 + z**2)









    Out[15]:





$\displaystyle 6.0$



In [16]:

    
def lap_vec(A, coords=(x, y, z), h_vec=(1, 1, 1)):
    """
    Laplacian of a vector function A.
    
    Parameters
    ==========
    A : Matrix, List
        Vector function to compute the curl from.
    coords : Tuple (3), optional
        Coordinates for the new reference system. This is an optional parameter
        it takes (x, y, z) as default.
    h_vec : Tuple (3), optional
        Scale coefficients for the new coordinate system. It takes
        (1, 1, 1), as default.
        
    Devuelve
    ========
    laplacian : Matrix (3, 1)
        Column vector with the components of the Laplacian.
    """  
    return grad(div(A, coords=coords, h_vec=h_vec), coords=coords, h_vec=h_vec) -\
           curl(curl(A, coords=coords, h_vec=h_vec), coords=coords, h_vec=h_vec)

Let's check with an example.



In [17]:

    
lap_vec([x**2, y**2, z**2])









    Out[17]:





$\displaystyle \left[\begin{matrix}2\\2\\2\end{matrix}\right]$

Plate equation in polar coordinates

The equation of a plate in polar coordinates is given by $D\nabla^4 w(\mathbf{x}) = -q(\mathbf{x})$, with $\nabla^4 u = \nabla^2 \nabla^2 u$.



In [18]:

    
w = symbols("w", cls=Function)
rho, phi, z = symbols("rho phi z", positive=True)
D, q = symbols("D q")



In [19]:

    
coords = (rho, phi, z)
h_vec = [1, rho, 1]



In [20]:

    
L1 = lap(w(rho, phi), coords=coords, h_vec=h_vec)
L1









    Out[20]:





$\displaystyle \frac{\rho \frac{\partial^{2}}{\partial \rho^{2}} w{\left(\rho,\phi \right)} + \frac{\partial}{\partial \rho} w{\left(\rho,\phi \right)}}{\rho} + \frac{\frac{\partial^{2}}{\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{2}}$



In [21]:

    
L2 = lap(L1, coords=coords, h_vec=h_vec)
L2









    Out[21]:





$\displaystyle \frac{\rho \left(\frac{\rho \frac{\partial^{4}}{\partial \rho^{4}} w{\left(\rho,\phi \right)} + 3 \frac{\partial^{3}}{\partial \rho^{3}} w{\left(\rho,\phi \right)}}{\rho} - \frac{2 \left(\rho \frac{\partial^{3}}{\partial \rho^{3}} w{\left(\rho,\phi \right)} + 2 \frac{\partial^{2}}{\partial \rho^{2}} w{\left(\rho,\phi \right)}\right)}{\rho^{2}} + \frac{\frac{\partial^{4}}{\partial \rho^{2}\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{2}} + \frac{2 \left(\rho \frac{\partial^{2}}{\partial \rho^{2}} w{\left(\rho,\phi \right)} + \frac{\partial}{\partial \rho} w{\left(\rho,\phi \right)}\right)}{\rho^{3}} - \frac{4 \frac{\partial^{3}}{\partial \rho\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{3}} + \frac{6 \frac{\partial^{2}}{\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{4}}\right) + \frac{\rho \frac{\partial^{3}}{\partial \rho^{3}} w{\left(\rho,\phi \right)} + 2 \frac{\partial^{2}}{\partial \rho^{2}} w{\left(\rho,\phi \right)}}{\rho} - \frac{\rho \frac{\partial^{2}}{\partial \rho^{2}} w{\left(\rho,\phi \right)} + \frac{\partial}{\partial \rho} w{\left(\rho,\phi \right)}}{\rho^{2}} + \frac{\frac{\partial^{3}}{\partial \rho\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{2}} - \frac{2 \frac{\partial^{2}}{\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{3}}}{\rho} + \frac{\frac{\rho \frac{\partial^{4}}{\partial \rho^{2}\partial \phi^{2}} w{\left(\rho,\phi \right)} + \frac{\partial^{3}}{\partial \rho\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho} + \frac{\frac{\partial^{4}}{\partial \phi^{4}} w{\left(\rho,\phi \right)}}{\rho^{2}}}{\rho^{2}}$



In [22]:

    
display(expand(D*L2 + q))









    





$\displaystyle D \frac{\partial^{4}}{\partial \rho^{4}} w{\left(\rho,\phi \right)} + \frac{2 D \frac{\partial^{3}}{\partial \rho^{3}} w{\left(\rho,\phi \right)}}{\rho} - \frac{D \frac{\partial^{2}}{\partial \rho^{2}} w{\left(\rho,\phi \right)}}{\rho^{2}} + \frac{2 D \frac{\partial^{4}}{\partial \rho^{2}\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{2}} + \frac{D \frac{\partial}{\partial \rho} w{\left(\rho,\phi \right)}}{\rho^{3}} - \frac{2 D \frac{\partial^{3}}{\partial \rho\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{3}} + \frac{4 D \frac{\partial^{2}}{\partial \phi^{2}} w{\left(\rho,\phi \right)}}{\rho^{4}} + \frac{D \frac{\partial^{4}}{\partial \phi^{4}} w{\left(\rho,\phi \right)}}{\rho^{4}} + q$



In [ ]: