Handling matrix of differential operator when using the Ritz method for an extremum problem The internal energy (or strain energy) of a typical structural member under a displacement field $\{u\}$ can be represented as:
$$
U = \frac{1}{2}\int_{V}[u]^T[B]^T[F]^T[B][u]dV
$$
where $V$ is the volume, $[F]$ the material constitutive matrix and $[B]$ the matrix of differential operators that appears when we apply the kinematic equations (strains as a function of displacements). In the Ritz method we assume a trial function that must satisfy only the geometric boundary conditions, then the displacement field can be written as:
$$
\{u\}=[g]\{c\}
$$
with $[g]$ containing the trial functions and $\{c\}$ the Ritz constants. The internal energy can be rewritten for the approximated displacement field:
$$
U = \frac{1}{2}\{c\}^T\left(\int_{V}[g]^T[B]^T[F]^T[B][g]dV\right)\{c\}
$$
This is usually error-prone when solved by hand, since there are numerous differential operations to be done. I am currently working on a differential operator in the Python module Sympy, but this should be much improved to avoid doing the differentiation before the matrices multiplications.
I would like to hear more approaches to handle matrix $[B]$.
 A: The following approach has been adopted to handle this type of equation involving a matrix of differential operators $[B]$. Go here to see the implementation in Python module SymPy.
Among many tests that have been performed, the following real example wil be used to explain the concept. These equations are used to calculate the buckling loads and the buckling modes of conical structures (taken from Shadmehri et al. 2012), where $D(x)=\partial/\partial{x}$ and $D(\theta)=\partial/\partial{\theta}$.
$$
[B] = \begin{bmatrix}
    D(x) &       0       &    0     &    0     &        0 \\
\frac1rsin\alpha &   \frac1rD(\theta)    & \frac1rcos\alpha &    0     &        0 \\
\frac1r D(\theta) & D(x)-\frac1rsin\alpha &    0     &    0     &        0 \\
    0    &       0       &    0     &   D(x)   &        0 \\
    0    &       0       &    0     & \frac1rsin\alpha & \frac1rD(\theta) \\
    0    &       0       &    0     & \frac1rD(\theta) & D(x)-\frac1x \\
    0    & -\frac1rcos\alpha     & \frac1rD(\theta) &    0     &        1 \\
    0    &       0       &   D(x)   &    1     &        0 \\
\end{bmatrix}
$$
with $\alpha$ a constant (the cone semi-angle). The $[g]$ matrix is a set of trial functions, for instance:
$$
[g]=\begin{bmatrix}
u_{11}  &0&0&0&0& \cdots & u_{21}& \cdots & u_{mn}\\
0& v_{11} &0&0&0& \cdots & v_{21}& \cdots & v_{mn}\\
0&0& w_{11} &0&0& \cdots & w_{21}& \cdots & w_{mn}\\
0&0&0& \phi_{x11} &0& \cdots & \phi_{x21}& \cdots & \phi_{xmn}\\
0&0&0&0&  \phi_{\theta 11}& \cdots &  \phi_{\theta 21}& \cdots & \phi_{\theta mn}\\
\end{bmatrix} 
\\
$$
with:
$$
\begin{matrix}
u_{ij} = sin(ia)cos(j\theta) & v_{ij} = cos(ia)sin(j\theta) &w_{ij}=sin(ia)sin(j\theta)\\
\phi_{xij} = cos(ia)sin(j\theta) & \phi_{\theta ij} = sin(ia)cos(j\theta)
\end{matrix}
$$
and:
$$
a=\pi (x-x_t)/(x_b-x_t)
$$
with $x_t$ and $x_b$ constants. Matrix $[F]$ is a constant $[8x8]$ matrix.
The following steps were successfully used to solve this problem:


*

*all functions involved in the multiplication $[B][g]$ (in this case $sin$ and $cos$) must be non-commutative

*the same for the symbols involved (in this case $x$, $\theta$, $r$, $alpha$, $x_t$ and $x_b$) 

*the differential operators $D()$ cannot be evaluated imediately when multiplying $[B][g]$. They should hold until the last expressions are obtained

*after obtaining the last expressions an evaluator function is applied, which search where $D()$ is encountered from the right to the left and applied the differentiation to the corresponding right-hand side

*the operator $D()$ must be programmed such that: $D(x)\cdot D(x)=D(x,x)$ (or $\partial/\partial{x} \cdot \partial/\partial{x}=\partial^2/\partial{x^2}$) to avoid incorrect evaluations or loosing parts of the expressions
