# Exercise: beam bending with the least squares method



The differential equation for beam bending
$$
(EI w'')'' - q = 0
$$
can be integrated analytically by specifying the linear loading as $q = q_0\frac{x}{L}$ and the boundary conditions
$$
w(0) = 0\\
w(L) = 0\\
w''(0)=0\\
w''(L) = 0
$$
to yield the deflection solution ("Biegelinie"):
$$
w(x) = \frac{q_0 L^4}{360 EI} \left[ 3\left(\frac{x}{L}\right)^5 - 10 \left(\frac{x}{L}\right)^3 + 7 \left(\frac{x}{L}\right)\right]
$$
We now seek to approximate this solution numerically using the least squares method.


```python
import numpy as np #numerical methods
import sympy as sp #symbolic operations
import matplotlib.pyplot as plt #plotting
sp.init_printing(use_latex='mathjax') #makes sympy output look nice

#Some plot settings
plt.style.use('seaborn-deep')
plt.rcParams['lines.linewidth']= 2.0
plt.rcParams['lines.color']= 'black'
plt.rcParams['legend.frameon']=True
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.family'] = 'serif'
plt.rcParams['legend.fontsize']=14
plt.rcParams['font.size'] = 14
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.left'] = True
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.axisbelow'] = True
```


```python
#Defining the geometric an material characteristics as symbolic quantities
L,q,EI,x = sp.symbols('L q_0 EI x')
```


```python
#Analytical solution to deflection
def deflection_analytical():
    a = x/L
    f = q*L**4/(360 * EI)
    return f*(3*a**5 - 10*a**3 + 7*a)
```


```python
deflection_analytical() #check definition
```




$$\frac{L^{4} q_{0}}{360 EI} \left(\frac{7 x}{L} - \frac{10 x^{3}}{L^{3}} + \frac{3 x^{5}}{L^{5}}\right)$$



Now, let's plot the analytical solution. For that purpose, we use some Python magic ("lambdify"). We sample the analytical solution for $x \in [0,L]$ at 100 points and plot the dimensionless deflection over the dimensionless length.


```python
lam_x = sp.lambdify(x, deflection_analytical(), modules=['numpy'])
#For the variable x the function deflection_analytical() will obtain something

x_vals = np.linspace(0, 1, 100)*L #This something is x from 0 to L
analytical = lam_x(x_vals) #We calculate the solution by passing x = [0,...,L] to deflection_analytical

plt.plot(x_vals/L, analytical/(L**4*q)*EI)
plt.xlabel('$x / L$')
plt.ylabel('$w / L^4 q_0  EI^{-1}$')
```

## Trigonometric Ansatz

Let's try the approximation
$$
    \tilde{w} = a_1 \sin \left(\pi \frac{x}{L}\right) + a_2 \sin \left(2\pi\frac{x}{L}\right)
$$


```python
a1, a2 = sp.symbols('a_1 a_2')#Defining the free values as new symbols
```


```python
def deflection_ansatz():
    return a1*sp.sin(sp.pi/L*x) + a2*sp.sin(2*sp.pi/L*x) #defining the approximate solution with the unknown coefficients
```

Now we substitute this solution into the fourth-order ODE for beam bending by differentiating it twice and adding the distributed loading:
$$
EI w^\text{IV} - q_0 \frac{x}{L} = 0 \text{ with } EI = \text{const.}
$$


```python
r = EI * deflection_ansatz().diff(x,4) - q * (x/L) #residual
r
```




$$\frac{\pi^{4} EI}{L^{4}} \left(a_{1} \sin{\left (\frac{\pi x}{L} \right )} + 16 a_{2} \sin{\left (\frac{2 \pi}{L} x \right )}\right) - \frac{q_{0} x}{L}$$



Now we perform the derivatives with respect to the Ansatz free values for obtaining the stationarity of $\int 1/2 r^2 \text{d} x$:


```python
dr_da1 = r.diff(a1)
dr_da2 = r.diff(a2)
dr_da1, dr_da2
```




$$\left ( \frac{\pi^{4} EI}{L^{4}} \sin{\left (\frac{\pi x}{L} \right )}, \quad \frac{16 EI}{L^{4}} \pi^{4} \sin{\left (\frac{2 \pi}{L} x \right )}\right )$$



This yields the two equations to be solved:


```python
stationarity_conditions = (sp.integrate(dr_da1*r,(x,0,L)),sp.integrate(dr_da2*r,(x,0,L)))
```


```python
coefficients = sp.solve(stationarity_conditions,a1,a2)
```


```python
coefficients
```




$$\left \{ a_{1} : \frac{2 L^{4} q_{0}}{\pi^{5} EI}, \quad a_{2} : - \frac{L^{4} q_{0}}{16 \pi^{5} EI}\right \}$$




```python
deflection_ansatz().subs([(a1,coefficients[a1]),(a2,coefficients[a2])]).simplify()
```




$$\frac{L^{4} q_{0}}{16 \pi^{5} EI} \left(32 \sin{\left (\frac{\pi x}{L} \right )} - \sin{\left (\frac{2 \pi}{L} x \right )}\right)$$



Now we're ready to plot the result and compare it to the analytical solution. 


```python
#We first substite the now known Ansatz free values into our Ansatz
z = sp.symbols('z')
w_numerical = deflection_ansatz().subs([(a1,coefficients[a1]),(a2,coefficients[a2]),(x,z*L)])
#We also made the coordinate dimensionless (x/L --> z) because of sympy problems

lam_x_num = sp.lambdify(z, w_numerical, modules=['numpy'])
#For the variable x the expression w_numerical will be given something
z_vals = np.linspace(0, 1,100) #This something is z from 0 to 1
numerical = lam_x_num(z_vals) #We calculate the solution by passing x = [0,...,L] to deflection_analytical
plt.plot(x_vals/L, analytical/(L**4*q)*EI,label='analytical')
plt.plot(z_vals, numerical/(L**4*q)*EI,label='numerical')
plt.legend()
plt.xlabel('$x\ /\ L$')
plt.ylabel('$w\ /\ L^4 q_0  EI^{-1}$')
plt.tight_layout()
plt.savefig('beam_least_squares.pdf')
```


```python
print("Maximum absolute error: ", np.max(np.abs(analytical/(L**4*q)*EI - numerical/(L**4*q)*EI)))
```

    Maximum absolute error:  3.41648400857485e-5


We can also plot and compare the bending moment. Let's first find the analytical expression by symbolically differentiating the deflection expression twice to obtain $M(x) = -EI w''(x)$:


```python
#analytical bending moment
moment_analytical = -deflection_analytical().diff(x,2)*EI
#numerical bending moment
moment_numerical = -deflection_ansatz().subs([(a1,coefficients[a1]),(a2,coefficients[a2])]).diff(x,2)*EI

#create lambdas for plotting along dimensionless length z
lam_x_analyt = sp.lambdify(z, moment_analytical.subs(x,z*L), modules=['numpy'])
lam_x_num = sp.lambdify(z, moment_numerical.subs(x,z*L), modules=['numpy'])
z_vals = np.linspace(0, 1,100)
analytical = lam_x_analyt(z_vals)
numerical = lam_x_num(z_vals)

#plot
plt.plot(x_vals/L, analytical/(L**2*q),label='analytical')
plt.plot(z_vals, numerical/(L**2*q),label='numerical')
plt.xlabel('$x\ /\ L$')
plt.ylabel('$M\ /\ L^2 q_0$')
plt.legend()
plt.tight_layout()
plt.savefig('beam_least_squares_moment.pdf')
```


```python
print("Maximum absolute error in bending moment: ", np.max(np.abs(analytical/(L**2*q) - numerical/(L**2*q))))
```

    Maximum absolute error in bending moment:  0.00384875067588690



```python

```

d8eaa44918008c163ae5536bedd6fc6f0e4731e8
ipynb
Jupyter Notebook
Beams/05_least_squares.ipynb
dominik-kern/Numerical_Methods_Introduction
09a0d6bd0ddbfc6e7f94b65516d9691766ed46ae
MIT
Beams/05_least_squares.ipynb
dominik-kern/Numerical_Methods_Introduction
09a0d6bd0ddbfc6e7f94b65516d9691766ed46ae
MIT
2022-01-04T19:02:05.000Z
2022-01-06T08:40:21.000Z
Beams/05_least_squares.ipynb
dominik-kern/Numerical_Methods_Introduction
09a0d6bd0ddbfc6e7f94b65516d9691766ed46ae
MIT
2020-12-03T13:01:55.000Z
2022-03-16T14:07:04.000Z
Qwen/Qwen-72B
1. YES
2. YES