We consider the eigenvalue problem associated to the longitudinal vibrations of an elastic rod
\begin{align} -\partial_{xx} u = \omega^2 u, \quad \Omega, \quad u_{|\partial \Omega} = 0 \end{align}where $\Omega = (0, 1)$.
The weak formulation, for $u \in \mathcal{V}_h$, is given by
\begin{align} a(u,v) = \omega^2 b(u,v), \quad \forall v \in \mathcal{V}_h \label{eq:vibration-wf} \end{align}where $a(u,v) = \int_{\Omega} \nabla v \cdot \nabla u~d\Omega$ and $b(u,v) = \int_{\Omega} v u~d\Omega$.
The bilinear forms $a$ and $b$ give the Stiffness and Mass matrices after discretization. Their associated GLT symbols are respectivaly $\mathfrak{s}_p$ and $\mathfrak{m}_p$. Therefor, the eigenvalues will be approximated by a uniform sampling of $\frac{\mathfrak{s}_p}{\mathfrak{m}_p}$.
from sympde.calculus import grad, dot
from sympde.topology import ScalarFunctionSpace
from sympde.topology import Domain
from sympde.topology import elements_of
from sympde.expr import BilinearForm
from sympde.expr import integral
DIM = 1
domain = Domain('\Omega', dim=DIM)
V = ScalarFunctionSpace('V', domain)
u,v = elements_of(V, names='u,v')
laplace = BilinearForm((u,v),
integral(domain, dot(grad(v), grad(u))))
mass = BilinearForm((u,v),
integral(domain, u*v))
# imports from gelato
from gelato.expr import gelatize
from gelato.printing import latex
We shall need to defined the following symbols
# imports from sympy
from sympy import Symbol
from sympy import ratsimp
from sympy import expand
tx = Symbol("tx")
nx = Symbol("nx", integer=True)
N = Symbol("N", integer=True)
wl = Symbol('\omega_l')
We compute the GLT expression by calling the function gelatize for a given degree $p$
p = 3
sl = gelatize(laplace, p)
sm = gelatize(mass, p)
expr = sl/sm
Print the latex expression using the latex function
from IPython.display import Math
Math(latex(ratsimp(expr/nx**2)))
expr = expr.subs({tx: wl / N})
expr = expr.subs({nx: N-p}) #N = nx + p
expr = expand(expr)
Let us print again our expression in Latex,
Math(latex(expr))
Well this can lead to a complicated expression. Let us take the expansion of the eigenvalues with respect to $\frac{1}{N}$.
# here we consider an expansion up to the forth order
order = 4
expr = expr.series(1/N, 0, order)
Now let us check again the eigenvalue expression
Math(latex(expr))
# css style
from IPython.core.display import HTML
def css_styling():
styles = open("../styles/custom.css", "r").read()
return HTML(styles)
css_styling()