# Exercise: beam bending with the subdomain 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 subdomain method.


```python

```


```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
```

By domain-wise integration of the ODE with domain-wise weighting by unity
$$
EI w^\text{IV} - q_0 \frac{x}{L} = 0 \text{ with } EI = \text{const.}
$$
we arrive at the boundary flux-based formulation
$$
    0 = EI (w'''|_{x_{i+1}} - w'''|_{x_{i}}) - \frac{q_0}{2L}(x_{i+1}^2 - x_i^2)
$$
In our case, to determine $a_1$ and $a_2$, we define the two domains:


```python
domain_1 = [0,L/2]
domain_2 = [L/2,L]
```

Hence, we get the system of equations


```python
#get w'''
d3w_dx3 = deflection_ansatz().diff(x,3)
#flux balance in domain 1
flux_domain_1 = EI*(d3w_dx3.subs(x,domain_1[1]) - d3w_dx3.subs(x,domain_1[0])) - q/(2*L)*(domain_1[1]**2 - domain_1[0]**2)
#flux balance in domain 2
flux_domain_2 = EI*(d3w_dx3.subs(x,domain_2[1]) - d3w_dx3.subs(x,domain_2[0])) - q/(2*L)*(domain_2[1]**2 - domain_2[0]**2)
#combining both domains into an equation system                    
flux_conditions = (flux_domain_1,flux_domain_2)
```


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


```python
coefficients
```




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




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




$$\frac{L^{4} q_{0}}{128 \pi^{3} 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_subdomain.pdf')
```


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

    Maximum absolute error:  0.00155626045741949


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_subdomain_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.0177679314001814


*Tasks:*
- What happens if the interfal boundaries are shifted?
- How could you improve the solution?


```python

```

bba5bc1504eef39df551549757c38a516b0e42ee
ipynb
Jupyter Notebook
Beams/06_subdomain.ipynb
dominik-kern/Numerical_Methods_Introduction
09a0d6bd0ddbfc6e7f94b65516d9691766ed46ae
MIT
Beams/06_subdomain.ipynb
dominik-kern/Numerical_Methods_Introduction
09a0d6bd0ddbfc6e7f94b65516d9691766ed46ae
MIT
2022-01-04T19:02:05.000Z
2022-01-06T08:40:21.000Z
Beams/06_subdomain.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