# Ritz method for a beam

**November, 2018**

We want to find a Ritz approximation of vibration frequencies/modes of a cantilever beam.
This is described by the following eigenvalue problem.

$$
\frac{\mathrm{d}^4w}{\mathrm{d}x^4} + \beta^4 w\, ,\quad
0 < x < L,\quad EI>0\, ,
$$

with

$$
w(0) = w'(0) = 0,\quad
\left(\frac{\mathrm{d}^2w}{\mathrm{d}x^2}\right)_{x=L} = 0,\quad
\left(\frac{\mathrm{d}^3 w}{\mathrm{d}x^3}\right)_{x=L} = 0\, ,
$$

and

$$\beta \equiv \left(\frac{\mu \omega^2}{EI}\right)^{1/4}\, .$$


```python
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
```


```python
%matplotlib notebook
init_printing()

# Graphics setup
gray = '#757575'
plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["text.color"] = gray
plt.rcParams["font.size"] = 12
plt.rcParams["xtick.color"] = gray
plt.rcParams["ytick.color"] = gray
plt.rcParams["axes.labelcolor"] = gray
plt.rcParams["axes.edgecolor"] = gray
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False
plt.rcParams["figure.figsize"] = 4, 3
```

The exact solution for this problem is

$$w_n(x) = A_1\left[(\cosh\beta_n x - \cos\beta_n x) +
\frac{\cos\beta_n L + \cosh\beta_n L}{\sin\beta_n L + \sinh\beta_n L}(\sin\beta_n x - \sinh\beta_n x)\right]\, ,$$

where $\beta_n$ is the $n$th root of $\cosh(\beta_n L)\cos(\beta_n L)  + 1$.


```python
x = symbols('x')
beta, L = symbols("beta L")
val1 = 1.8750909912
val2 = 4.6941049111
```


```python
def plot_expr(expr, x, rango=(0, 1), ax=None, linestyle="solid"):
    """Plot SymPy expressions of a single variable"""
    expr_num = lambdify(x, expr, "numpy")
    x0 = rango[0]
    x1 = rango[1]
    x_num = np.linspace(0, 1, 101)
    if ax is None:
        plt.figure()
        ax = plt.gca()
    ax.plot(x_num, expr_num(x_num), linestyle=linestyle)
```

## Conventional formulation

The quadratic functional for this problem is

$$J[u] = \int\limits_0^L \left[\left(\frac{\mathrm{d}^2 w}{\mathrm{d}x^2}\right)^2 + \beta^4 w^2\right]\mathrm{d}x\, ,$$

and the weak problem $B(v, w) = \beta^2 A(v, w)$, with

$$
B(v, u) = \int\limits_0^L \frac{\mathrm{d}^2 v}{\mathrm{d}x^2}\frac{\mathrm{d}^2 u}{\mathrm{d}x^2}\mathrm{d}x\, ,\quad
A(v, u) = \int\limits_0^L  vw\mathrm{d}x\, .
$$


```python
def quad_fun(x, w, L):
    F1 = diff(w, x, 2)**2
    F2 = w**2
    U = integrate(F1, (x, 0, L))
    T = integrate(F2, (x, 0, L))
    return U, T
    
```


```python
def ritz_conventional(x, L, nterms):
    a = symbols("a0:%i"%(nterms))
    w = sum(a[k]*x**(k + 2) for k in range(nterms))
    U, T = quad_fun(x, w, L)
    K = Matrix(nterms, nterms, lambda i, j: U.diff(a[i], a[j]))
    M = Matrix(nterms, nterms, lambda i, j: T.diff(a[i], a[j]))
    return K, M
```


```python
K, M = ritz_conventional(x, 1, 2)

```


```python
Kaux = M.inv() * K
vals = list(Kaux.eigenvals())
vals
```


```python
nvals = [N(val**0.25) for val in vals]
nvals
```

## Lagrange multiplier formulation

We can write the problem as minimizing the functional

$$J(\psi, w) = \int\limits_0^L\left[\left(\frac{\mathrm{d} \psi}{\mathrm{d}x}\right)^2 + \beta^4
w^2\right]\mathrm{d}x\, ,$$

subject to

$$G(\psi, w) \equiv \psi + \frac{\mathrm{d}w}{\mathrm{d}x} = 0\, .$$

The Lagrangian is given by

$$L(\psi, w, \lambda) = \int\limits_0^L\left[\left(\frac{\mathrm{d} \psi}{\mathrm{d}x}\right)^2 + \beta^4 w^2\right]\mathrm{d}x + \int\limits_0^L \lambda\left(\psi + \frac{\mathrm{d}w}{\mathrm{d}x}\right)\mathrm{d}x\, , $$

where $\lambda$ is the Lagrange multiplier, which in this case represents the shear force.


```python
def lagran(x, w, psi, lamda, L):
    F1 = diff(psi, x)**2 + lamda*(psi + diff(w, x))
    F2 = w**2
    U = integrate(F1, (x, 0, L))
    T = integrate(F2, (x, 0, L))
    return U, T
    
```


```python
def ritz_multiplier(x, L, nterms):
    a = symbols("a0:%i"%(nterms))
    b = symbols("b0:%i"%(nterms))
    c = symbols("c0:%i"%(nterms))
    var = a + b + c
    w = sum(a[k]*x**(k + 1) for k in range(nterms))
    psi = sum(b[k]*x**(k + 1) for k in range(nterms))
    lamda = sum(c[k]*x**k for k in range(nterms))
    U, T = lagran(x, w, psi, lamda, L)
    K = Matrix(3*nterms, 3*nterms, lambda i, j: U.diff(var[i], var[j]))
    M = Matrix(3*nterms, 3*nterms, lambda i, j: T.diff(var[i], var[j]))
    return K, M
```


```python
K, M = ritz_multiplier(x, 1, 2)
```


```python
K
```




$$\left[\begin{matrix}0 & 0 & 0 & 0 & 1 & \frac{1}{2}\\0 & 0 & 0 & 0 & 1 & \frac{2}{3}\\0 & 0 & 2 & 2 & \frac{1}{2} & \frac{1}{3}\\0 & 0 & 2 & \frac{8}{3} & \frac{1}{3} & \frac{1}{4}\\1 & 1 & \frac{1}{2} & \frac{1}{3} & 0 & 0\\\frac{1}{2} & \frac{2}{3} & \frac{1}{3} & \frac{1}{4} & 0 & 0\end{matrix}\right]$$




```python
M
```




$$\left[\begin{matrix}\frac{2}{3} & \frac{1}{2} & 0 & 0 & 0 & 0\\\frac{1}{2} & \frac{2}{5} & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$




```python
Maux = K.inv() * M
vals = list(Maux.eigenvals())
nvals = [N(1/val**0.25) for val in vals if val != 0]
nvals
```

# The penalty function formulation

The augmented functional for this formulation is given by

$$P_K (\psi, w) = J(\psi, w) + \frac{K}{2}\int\limits_0^L \left(\psi + \frac{\mathrm{d}w}{\mathrm{d}x}\right)^2\mathrm{d}x\, ,$$

where $K$ is the penalty parameter.



```python
def augmented(x, w, psi, K, L):
    F1 = diff(psi, x)**2 + S(K)/2*(psi + diff(w, x))**2
    F2 = w**2
    U = integrate(F1, (x, 0, L))
    T = integrate(F2, (x, 0, L))
    return U, T
```


```python
def ritz_penalty(x, K, L, nterms):
    a = symbols("a0:%i"%(nterms))
    b = symbols("b0:%i"%(nterms))
    var = a + b
    w = sum(a[k]*x**(k + 1) for k in range(nterms))
    psi = sum(b[k]*x**(k + 1) for k in range(nterms))
    U, T = augmented(x, w, psi, K, L)
    K = Matrix(2*nterms, 2*nterms, lambda i, j: U.diff(var[i], var[j]))
    M = Matrix(2*nterms, 2*nterms, lambda i, j: T.diff(var[i], var[j]))
    return K, M
```


```python
K, M = ritz_penalty(x, 100, 1, 2)
```


```python
K
```




$$\left[\begin{matrix}100 & 100 & 50 & \frac{100}{3}\\100 & \frac{400}{3} & \frac{200}{3} & 50\\50 & \frac{200}{3} & \frac{106}{3} & 27\\\frac{100}{3} & 50 & 27 & \frac{68}{3}\end{matrix}\right]$$




```python
M
```




$$\left[\begin{matrix}\frac{2}{3} & \frac{1}{2} & 0 & 0\\\frac{1}{2} & \frac{2}{5} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]$$




```python
Maux = K.inv() * M
vals = list(Maux.eigenvals())
nvals = [re(N(1/val**0.25)) for val in vals if val != 0]
nvals
```

## Mixed formulation

The mixed formulation involves rewriting a given higher order equation as a pair of lower
order equations by introducing secondary dependent variables. The original equation can be
decomposed into

$$
M(x)= \frac{\mathrm{d}^2 w}{\mathrm{d}x^2}\, ,\quad
\frac{\mathrm{d}^2M(x)}{\mathrm{d}x^2} = -\beta^4 w\, ,\quad 0<x<L\, .
$$

The functional in this case is

$$
I(w, M) = \int\limits_0^L\left(\frac{\mathrm{d}w}{\mathrm{d}x}\frac{\mathrm{d}M}{\mathrm{d}x}
+ \frac{M^2}{2} - \beta^4 w^2\right)\mathrm{d}x
$$


```python
def mixed_fun(x, w, M, L):
    F1 = diff(w, x)*diff(M, x) + M**2/2
    F2 = -w**2
    U = integrate(F1, (x, 0, L))
    T = integrate(F2, (x, 0, L))
    return U, T
    
```


```python
def ritz_mixed(x, L, nterms):
    a = symbols("a0:%i"%(nterms))
    b = symbols("b0:%i"%(nterms))
    var = a + b
    w = sum(a[k]*x**(k + 1) for k in range(nterms))
    M = sum(b[k]*(x - L)**(k + 1) for k in range(nterms))
    display(w, M)
    U, T = mixed_fun(x, w, M, L)
    Kmat = Matrix(2*nterms, 2*nterms, lambda i, j: U.diff(var[i], var[j]))
    Mmat = Matrix(2*nterms, 2*nterms, lambda i, j: T.diff(var[i], var[j]))
    return Kmat, Mmat
```


```python
K, M = ritz_mixed(x, 1, 2)
```


```python
K
```




$$\left[\begin{matrix}0 & 0 & 1 & -1\\0 & 0 & 1 & - \frac{2}{3}\\1 & 1 & \frac{1}{3} & - \frac{1}{4}\\-1 & - \frac{2}{3} & - \frac{1}{4} & \frac{1}{5}\end{matrix}\right]$$




```python
M
```




$$\left[\begin{matrix}- \frac{2}{3} & - \frac{1}{2} & 0 & 0\\- \frac{1}{2} & - \frac{2}{5} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]$$




```python
Maux = K.inv() * M
vals = list(Maux.eigenvals())
nvals = [re(N(1/val**0.25)) for val in vals if val != 0]
nvals
```


```python

```

8f30562a8ee9823bdfb96482cd5c5cc3e754f17b
ipynb
Jupyter Notebook
variational/ritz_vibration_beam.ipynb
nicoguaro/FEM_resources
32f032a4e096fdfd2870e0e9b5269046dd555aee
MIT
2015-11-06T16:59:39.000Z
2022-02-25T18:18:49.000Z
variational/ritz_vibration_beam.ipynb
oldninja/FEM_resources
e44f315be217fd78ba95c09e3c94b1693773c047
MIT
variational/ritz_vibration_beam.ipynb
oldninja/FEM_resources
e44f315be217fd78ba95c09e3c94b1693773c047
MIT
2018-06-24T22:12:00.000Z
2022-01-12T15:57:37.000Z
Qwen/Qwen-72B
1. YES
2. YES