User:Eml5526.s11.team5/sub6

= 6.1 [ Solving G1DM1.0/D1 using weak form with QLEBF (Quadratic Lagrangian Element Basis Functions)]  =

P.D.E

 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 6.1.1)
 * }

Boundary Conditions

 * {| style="width:100%" border="0"

$$u \ at \ {x=0} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 6.1.2)
 * }


 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 1 \ is \ \frac{12}{5} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 6.1.3)
 * }

Find: Approximate solution using weak form and its comparison to exact solution

 * Consider basis function QLEBF satisfying constraint basis solution


 * Compute $$ u_n^h (x = 0.5) \quad $$ for $$ n = 2, 4, 6..... \quad $$ until convergence of $$ u^h (0.5)\ to \ 0(10^{-6}) $$


 * error


 * {| style="width:100%" border="0"

$$ e_n(0.5) = u_n(0.5) - u_n^h(0.5) \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                       (Eq. 6.1.4)
 * }


 * Plot $$ e_n(0.5) \quad $$ vs. $$ n \quad $$

Solution
Strong form


 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                       (Eq. 6.1.5)
 * }

Boundary conditions


 * {| style="width:100%" border="0"

$$u \ at \ {x=0} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 6.1.6)
 * }


 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 1 \ is \ \frac{12}{5} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 6.1.7)
 * }

Finding of the Exact Solution
The exact solution of the displacement can be found as follows.

Integrating eqn 5.2.5, we get


 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = c_1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 6.1.8)
 * }

Substituting the natural B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ c_1 = 14.5 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.9)
 * }

Therefore eqn 5.1.8 becomes,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = 14.5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.10)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\frac{\partial u}{\partial x} = \frac{14.5-\frac{5x^2}{2}}{(2+3x)} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.11)
 * }

Integrating the above equation, we get,Reference : www.wolframalpha.com

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x + 482 log (2+3x) +60 ) + c_2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.12)
 * }

Substituting the Essential B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$c_2 = 0.3509542497 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.13)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x + 482 log (2+3x) +60 ) + 0.3509542497 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.14)
 * }

Finding of the Approximate Solution using Weak form with QLEBF
Weak Form

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \sum_{i} c_i \left [ \sum_{j} \left \{ \tilde{K_{ij}} \tilde{d_j} \right \} - \tilde{F_i} \right ] = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.15)
 * }

where,

<span id="(1)">
 * {| style="width:100%" border="0"

$$\tilde{K_{ij}} = \int_{0}^{1} b_i' \ a_2 \ b_j ' \ dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.16)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{F_i} = b_i (\alpha )h + \int_{0}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.17)
 * }

For Element 1
<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

Ist Basis Function: $$b_{1}^{(1)}=L_{1,3}=\frac{(x-x_{2}).(x-x_{3})}{(x_{1}-x_{2}).(x_{1}-x_{3})}=\frac{(x-0.25).(x-0.5)}{(0-0.25).(0-0.5)}=8x^{2}-6x+1 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.18)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IInd Basis Function: $$b_{2}^{(1)}=L_{2,3}=\frac{(x-x_{1}).(x-x_{3})}{(x_{2}-x_{1}).(x_{2}-x_{3})}=\frac{(x-0).(x-0.5)}{(0.25-0).(0.25-0.5)}=-16x^{2}+8x $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.19)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IIIrd Basis Function : $$b_{3}^{(1)}=L_{3,3}=\frac{(x-x_{1}).(x-x_{2})}{(x_{3}-x_{1}).(x_{3}-x_{2})}=\frac{(x-0).(x-0.25)}{(0.5-0).(0.5-0.25)}=8x^{2}-2x $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.20)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$${b_{1}^{(1)}}'=16x-6, \quad\quad {b_{2}^{(1)}}'=-32x+8 , \quad\quad {b_{3}^{(1)}}'=16x-2 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.21)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

As we know from the definition of K matrix <span id="(1)">
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{ij}^{(1)}}=\int_{\alpha }^{\beta }{b_{i}^{(1)}}'.a_{2}.{b_{j}^{(1)}}'dx $$
 * style="width:95%" |
 * style="width:95%" |

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.21)
 * }
 * {| style="width:100%" border="0"

So, the elements of K matrix are:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{11}^{(1)}}=\int_{0 }^{0.5 }(16x-6).(2+3x).(16x-6)dx=10.8333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.21)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{12}^{(1)}}=\mathbf{K_{21}^{(1)}}=\int_{0 }^{0.5 }(16x-6).(2+3x).(-32x+8)dx=-12.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.22)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{13}^{1}}=\mathbf{K_{31}^{(1)}}=\int_{0 }^{0.5 }(16x-6).(2+3x).(16x-2)dx=1.83333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.23)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{22}^{(1)}}=\int_{0 }^{0.5 }(-32x+8).(2+3x).(-32x+8)dx=29.3333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.24)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{23}^{(1)}}=\int_{0 }^{0.5 }(-32x+8).(2+3x).(16x-2)dx=-16.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.25)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{33}^{(1)}}=\int_{0 }^{0.5 }(16x-2).(2+3x).(16x-2)dx=14.83333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.26)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

S0, the Element Stiffness Matrix or Local Stiffness Matrix for the Element 1 is :
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{k^{(1)}}=\begin{bmatrix} 10.8333 & -12.6667 & 1.83333 \\ -12.6667 & 29.3333 & -16.6667 \\  1.83333 & -16.6667  & 14.83333 \end{bmatrix}

$$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.27)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

As we know from the definition of weak form the Force Matrix is given by:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{F_i} = b_i (\alpha )h + \int_{0}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.28)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the elements of force matrix are: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ f_{1}^{(1)}=\int_{0}^{0.5}(8x^{2}-6x+1).(5x)dx=0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{2}^{(1)}=\int_{0}^{0.5}(-16x^{2}+8x).(5x)dx=0.416667 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{3}^{(1)}=\int_{0}^{0.5}(8x^{2}-2x).(5x)dx=0.2083333 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the Element Force Matrix for Element 1 is:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{f^{(1)}}=\begin{bmatrix} 0\\ 0.416667\\ 0.208333 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.29)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The transformation from Local to Global is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

So, the global K and F matrices for Element 1 are given as:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{e}}= \mathbf{L^{e^{T}}}.\mathbf{k^{e}}.\mathbf{L^{e}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.30)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{F^{e}}=\mathbf{L^{e^{T}}}.\mathbf{f^{e}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.31)
 * }
 * {| style="width:100%" border="0"

For Element 1 $$ \mathbf{L^{e}}=\mathbf{L^{1}} $$.
 * style="width:95%" |
 * style="width:95%" |
 * }

$$ \mathbf{L^{e}} \quad $$ For Element 1
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{L^{1}}=\begin{bmatrix} 1 &0 &0  &0  &0 \\ 0 &1  &0  &0  &0 \\ 0 &0  &1  &0  &0 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.32)
 * }

Global K,F matrices For Element 1
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{1}}= \mathbf{L^{1^{T}}}.\mathbf{k^{1}}.\mathbf{L^{1}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{K^{1}}=\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} 10.8333 & -12.6667 & 1.83333 \\ -12.6667 & 29.3333 & -16.6667 \\  1.83333 & -16.6667  & 14.83333 \end{bmatrix}\begin{bmatrix} 1 &0 &0  &0  &0 \\ 0 &1  &0  &0  &0 \\ 0 &0  &1  &0  &0 \end{bmatrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{K^{1}}=\begin{bmatrix} 10.8333 &-12.6667 &1.8333    &0  &0 \\ -12.667 &29.33333  &-16.6667  &0  &0 \\ 1.8333  & -16.6667 & 14.83333 &0  &0 \\ 0 &0  &0  &0  &0 \\ 0 &0  &0  &0  &0 \end{bmatrix}$$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.32)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{F^{1}}=\mathbf{L^{1^{T}}}.\mathbf{f^{1}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{F^{1}}=\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} 0\\ 0.416667\\ 0.208333 \end{bmatrix} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{F^{1}}=\begin{bmatrix} 0\\ 0.416667\\ 0.208333\\ 0\\ 0 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.33)
 * }

For Element 2
<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

Ist Basis Function: $$b_{1}^{(2)}=\frac{(x-x_{4}).(x-x_{5})}{(x_{3}-x_{4}).(x_{3}-x_{5})}=\frac{(x-0.75).(x-1)}{(0.5-0.75).(0.5-1)}=8x^{2}-14x+6 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.34)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IInd Basis Function: $$b_{2}^{(2)}=\frac{(x-x_{3}).(x-x_{5})}{(x_{4}-x_{3}).(x_{4}-x_{5})}=\frac{(x-0.5).(x-1)}{(0.75-0.5).(0.75-1)}=-16x^{2}+24x-8 $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.35)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IIIrd Basis Function : $$b_{3}^{(2)}=\frac{(x-x_{3}).(x-x_{4})}{(x_{5}-x_{3}).(x_{5}-x_{4})}=\frac{(x-0.5).(x-0.75)}{(1-0.5).(1-0.75)}=8x^{2}-10x+3 $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.36)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$${b_{1}^{(2)}}'=16x-14, \quad\quad {b_{2}^{(2)}}'=-32x+24 , \quad\quad {b_{3}^{(2)}}'=16x-10 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.37)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

As we know from the definition of K matrix <span id="(1)">
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{ij}^{(2)}}=\int_{\alpha }^{\beta }{b_{i}^{(2)}}'.a_{2}.{b_{j}^{(2)}}'dx $$
 * style="width:95%" |
 * style="width:95%" |

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.38)
 * }
 * {| style="width:100%" border="0"

So, the elements of local k matrix (for Element 2) are:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{11}^{(2)}}=\int_{0.5 }^{1 }(16x-14).(2+3x).(16x-14)dx=17.8333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.39)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{k_{12}^{(2)}}=\mathbf{k_{21}^{(2)}}=\int_{0.5 }^{1}(16x-14).(2+3x).(-32x+24)dx=-20.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.40)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{13}^{(2)}}=\mathbf{k_{31}^{(2)}}=\int_{0.5 }^{1}(16x-14).(2+3x).(16x-10)dx=2.83333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.41)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{k_{22}^{(2)}}=\int_{0.5 }^{1}(-32x+24).(2+3x).(-32x+24)dx=45.3333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.42)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{23}^{(2)}}=\int_{0.5 }^{1 }(-32x+24).(2+3x).(16x-10)dx=-24.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.42)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{33}^{(2)}}=\int_{0.5 }^{1}(16x-10).(2+3x).(16x-10)dx=21.83333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.43)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

S0, the Element Stiffness Matrix or Local Stiffness Matrix for the Element 2 is :
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{k^{(2)}}=\begin{bmatrix} 17.8333 & -20.6667 & 2.83333 \\ -20.6667 & 45.3333 & -24.6667 \\  2.83333 & -24.6667  & 21.83333 \end{bmatrix}

$$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.44)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

As we know from the definition of weak form the Force Matrix is given by:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{f_i} = b_i (\alpha )h + \int_{0.5}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.45)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the elements of force matrix are: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ f_{1}^{(2)}=\int_{0.5}^{1}(8x^{2}-14x+6).(5x)dx=0.28333 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{2}^{(2)}=\int_{0.5}^{1}(-16x^{2}+24x-8).(5x)dx=1.25 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{3}^{(2)}=(1).(12)\int_{0.5}^{1}(8x^{2}-10x+3).(5x)dx=12.416667 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the Element Force Matrix for Element 2 is:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{f^{(2)}}=\begin{bmatrix} 0.208333\\ 1.25\\ 12.416667 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.46)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The transformation from Local to Global is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

Global K,F matrices For Element 2
So, the global K and F matrices for Element 2 are given as:
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{(e)}}= \mathbf{L^{(e)^{T}}}.\mathbf{k^{(e)}}.\mathbf{L^{(e)}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.47)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{F^{(e)}}=\mathbf{L^{(e)^{T}}}.\mathbf{f^{(e)}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.48)
 * }
 * {| style="width:100%" border="0"

For Element 2 $$ \mathbf{L^{(e)}}=\mathbf{L^{(2)}} $$.
 * style="width:95%" |
 * style="width:95%" |
 * }

$$ \mathbf{L^{(e)}} $$ For Element 2
<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{L^{(2)}}=\begin{bmatrix} 0 &0 &1  &0  &0 \\ 0 &0  &0  &1  &0 \\ 0 &0  &0  &0  &1 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.49)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{(2)}}= \mathbf{L^{(2)^{T}}}.\mathbf{k^{(2)}}.\mathbf{L^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{K^{(2)}}=\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 17.8333 & -20.6667 & 2.83333 \\ -20.6667 & 45.3333 & -24.6667 \\  2.83333 & -24.6667  & 21.83333 \end{bmatrix}\begin{bmatrix} 0 &0 &1  &0  &0 \\ 0 &0  &0  &1  &0 \\ 0 &0  &0  &0  &1 \end{bmatrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{K^{(2)}}=\begin{bmatrix} 0 &0 &0   &0  &0 \\ 0 &0  &0   &0  &0 \\ 0 &0  & 17.83333 &-20.6667  &2.83333 \\ 0 &0  &-20.6667  &45.3333  &-24.6667 \\ 0 &0  &2.83333  &-24.6667  &21.83333 \end{bmatrix}$$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.50)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{F^{(2)}}=\mathbf{L^{(2)^{T}}}.\mathbf{f^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{F^{(2)}}=\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 0.208333\\ 1.25\\ 12.46667 \end{bmatrix} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{F^{(2)}}=\begin{bmatrix} 0\\ 0\\ 0.208333\\ 1.25\\ 12.46667 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.51)
 * }

Global K matrix for 2 elements
<span id="(1)">
 * {| style="width:100%" border="0"

The Global K matrix for 2 elements is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{K}}=\sum_{e=1}^{2}\mathbf{K^{e}} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.52)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow \mathbf{\tilde{K}}=\mathbf{\tilde{K}^{(1)}}+\mathbf{\tilde{K}^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow \tilde{K}=\begin{bmatrix} 10.8333 &-12.6667 &1.8333    &0  &0 \\ -12.667 &29.33333  &-16.6667  &0  &0 \\ 1.83333 &-16.6667  & 32.6666 &-20.6667  &2.83333 \\ 0 &0  &-20.6667  &45.3333  &-24.6667 \\ 0 &0  &2.83333  &-24.6667  &21.83333 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.53)
 * }

Global F matrix for 2 elements
<span id="(1)">
 * {| style="width:100%" border="0"

The Global F matrix for 2 elements is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{F}}=\sum_{e=1}^{2}\mathbf{F^{e}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.54)
 * }
 * {| style="width:100%" border="0"

$$ \Rightarrow \mathbf{\tilde{F}}=\mathbf{\tilde{F}^{(1)}}+\mathbf{\tilde{F}^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow \tilde{F}=\begin{bmatrix} 0\\ 0.416667\\ 0.416667\\ 1.25\\ 12.41667 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.55)
 * }

Global d matrix for 2 elements
<span id="(1)">
 * {| style="width:100%" border="0"

As we know, <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{K}}.\mathbf{\tilde{d}}=\mathbf{\tilde{F}} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$\Rightarrow \mathbf{\tilde{d}}=\mathbf{\tilde{K}^{-1}}.\mathbf{\tilde{F}} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

So. the $$ \mathbf{\tilde{d}} $$ is :
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{d}}=\begin{bmatrix} 4\\ 5.5312\\ 6.6698\\ 7.5446\\ 8.2268 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.56)
 * }

At $$ x = 0.5 \quad $$ ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h(0.5) = 6.669855 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.57)
 * }

The exact solution can be found from the eqn. 6.1.14. It is given by ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u(0.5) = 6.671155 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Hence the error we get for $$ no. \quad of \quad elements \quad (n) = 2 $$ is ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$u^h(0.5) - u(0.5) = 0.0013 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.1.58)
 * }

MATLAB code with output
The K,d and F matrix for 2 elements only are as follows:

K =

10.8333 -12.6667    1.8333         0         0  -12.6667   29.3333  -16.6667         0         0    1.8333  -16.6667   32.6667  -20.6667    2.8333         0         0  -20.6667   45.3333  -24.6667         0         0    2.8333  -24.6667   21.8333

F =

[0   0.416666666666667    0.416666666666667    1.25000000000000    12.416666666666667]^T

d =

[4   5.53116172316385     6.66984463276837     7.54458967938842     8.22680850223224]^T

The outputs for convergence are as follows:

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=2 = 6.669845e+000

Error at n=2 is 1.311013e-003

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=4 = 6.671055e+000

Error at n=4 is 1.003161e-004

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=6 = 6.671135e+000

Error at n=6 is 2.076369e-005

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=8 = 6.671149e+000

Error at n=8 is 6.685488e-006

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=10 = 6.671153e+000

Error at n=10 is 2.761216e-006

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=12 = 6.671154e+000

Error at n=12 is 1.337704e-006

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=14 = 6.671155e+000

Error at n=14 is 7.240625e-007

u,u^h vs n Plot
Eml5526.s11.team5.Dahiya 16:45, 6 April 2011 (UTC)

u,u^h vs x Plot
Eml5526.s11.team5.Dahiya 16:43, 6 April 2011 (UTC)

= 6.2 [ Solving G1DM1.0/D1 using weak form with QLEBF (Quadratic Lagrangian Element Basis Functions)]  =

P.D.E
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.1)
 * }

Boundary Conditions
<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=1} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 0 \ is \ -6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.3)
 * }

Find: Approximate solution using weak form and its comparison to exact solution
 Find approximate solution using weak form and compare it to exact one


 * Consider QLEBF basis function satisfying constraint basis solution


 * Compute $$ u_n^h (x = 0.5) \quad $$ for $$ n = 2, 4, 6..... \quad $$ until convergence of $$ u^h (0.5)\ to \ 0(10^{-6}) $$


 * error

<span id="(1)">
 * {| style="width:100%" border="0"

$$ e_n(0.5) = u_n(0.5) - u_n^h(0.5) \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                       (Eq. 6.2.4)
 * }


 * Plot $$ e_n(0.5) \quad $$ vs. $$ n \quad $$

Solution
Strong form

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                       (Eq. 6.2.5)
 * }

Boundary conditions

<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=1} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.6)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 0 \ is \ -6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.7)
 * }

Finding of the Exact Solution
The exact solution of the displacement can be found as follows.

Differentiating eqn 5.2.5, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = c_1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.8)
 * }

Substituting the natural B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ c_1 = -12 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.9)
 * }

Therefore eqn 5.2.8 becomes,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = -12 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.10)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\frac{\partial u}{\partial x} = \frac{-12-\frac{5x^2}{2}}{(2+3x)} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.11)
 * }

Integrating the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x - 472 log (2+3x) +60 ) + c_2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.12)
 * }

Substituting the Essential B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$c_2 = 10.3394 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.13)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x - 472 log (2+3x) +60 ) + 10.3393953 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.14)
 * }

Finding of the Approximate Solution using Weak form with QLEBF
Weak Form

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \sum_{i} c_i \left [ \sum_{j} \left \{ \tilde{K_{ij}} \tilde{d_j} \right \} - \tilde{F_i} \right ] = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.15)
 * }

where,

<span id="(1)">
 * {| style="width:100%" border="0"

$$\tilde{K_{ij}} = \int_{0}^{1} b_i' \ a_2 \ b_j ' \ dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.16)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{F_i} = b_i (\alpha )h + \int_{0}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.17)
 * }

For Element 1
<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

Ist Basis Function: $$b_{1}^{(1)}=L_{1,3}=\frac{(x-x_{2}).(x-x_{3})}{(x_{1}-x_{2}).(x_{1}-x_{3})}=\frac{(x-0.25).(x-0.5)}{(0-0.25).(0-0.5)}=8x^{2}-6x+1 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.18)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IInd Basis Function: $$b_{2}^{(1)}=L_{2,3}=\frac{(x-x_{1}).(x-x_{3})}{(x_{2}-x_{1}).(x_{2}-x_{3})}=\frac{(x-0).(x-0.5)}{(0.25-0).(0.25-0.5)}=-16x^{2}+8x $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.19)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IIIrd Basis Function : $$b_{3}^{(1)}=L_{3,3}=\frac{(x-x_{1}).(x-x_{2})}{(x_{3}-x_{1}).(x_{3}-x_{2})}=\frac{(x-0).(x-0.25)}{(0.5-0).(0.5-0.25)}=8x^{2}-2x $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.20)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$${b_{1}^{(1)}}'=16x-6, \quad\quad {b_{2}^{(1)}}'=-32x+8 , \quad\quad {b_{3}^{(1)}}'=16x-2 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.21)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

As we know from the definition of K matrix <span id="(1)">
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{ij}^{(1)}}=\int_{\alpha }^{\beta }{b_{i}^{(1)}}'.a_{2}.{b_{j}^{(1)}}'dx $$
 * style="width:95%" |
 * style="width:95%" |

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.21)
 * }
 * {| style="width:100%" border="0"

So, the elements of K matrix are:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{11}^{(1)}}=\int_{0 }^{0.5 }(16x-6).(2+3x).(16x-6)dx=10.8333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.21)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{12}^{(1)}}=\mathbf{K_{21}^{(1)}}=\int_{0 }^{0.5 }(16x-6).(2+3x).(-32x+8)dx=-12.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.22)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{13}^{1}}=\mathbf{K_{31}^{(1)}}=\int_{0 }^{0.5 }(16x-6).(2+3x).(16x-2)dx=1.83333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.23)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{22}^{(1)}}=\int_{0 }^{0.5 }(-32x+8).(2+3x).(-32x+8)dx=29.3333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.24)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{23}^{(1)}}=\int_{0 }^{0.5 }(-32x+8).(2+3x).(16x-2)dx=-16.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.25)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K_{33}^{(1)}}=\int_{0 }^{0.5 }(16x-2).(2+3x).(16x-2)dx=14.83333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.26)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

S0, the Element Stiffness Matrix or Local Stiffness Matrix for the Element 1 is :
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{k^{(1)}}=\begin{bmatrix} 10.8333 & -12.6667 & 1.83333 \\ -12.6667 & 29.3333 & -16.6667 \\  1.83333 & -16.6667  & 14.83333 \end{bmatrix}

$$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.27)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

As we know from the definition of weak form the Force Matrix is given by:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{F_i} = b_i (\alpha )h + \int_{0}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.28)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the elements of force matrix are: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ f_{1}^{(1)}=(1).(12)+\int_{0}^{0.5}(8x^{2}-6x+1).(5x)dx=12 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{2}^{(1)}=\int_{0}^{0.5}(-16x^{2}+8x).(5x)dx=0.416667 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{3}^{(1)}=\int_{0}^{0.5}(8x^{2}-2x).(5x)dx=0.2083333 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the Element Force Matrix for Element 1 is:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{f^{(1)}}=\begin{bmatrix} 12\\ 0.416667\\ 0.208333 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.29)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The transformation from Local to Global is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

So, the global K and F matrices for Element 1 are given as:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{e}}= \mathbf{L^{e^{T}}}.\mathbf{k^{e}}.\mathbf{L^{e}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.30)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{F^{e}}=\mathbf{L^{e^{T}}}.\mathbf{f^{e}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.31)
 * }
 * {| style="width:100%" border="0"

For Element 1 $$ \mathbf{L^{e}}=\mathbf{L^{1}} $$.
 * style="width:95%" |
 * style="width:95%" |
 * }

$$ \mathbf{L^{e}} \quad $$ For Element 1
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{L^{1}}=\begin{bmatrix} 1 &0 &0  &0  &0 \\ 0 &1  &0  &0  &0 \\ 0 &0  &1  &0  &0 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.32)
 * }

Global K,F matrices For Element 1
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{1}}= \mathbf{L^{1^{T}}}.\mathbf{k^{1}}.\mathbf{L^{1}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{K^{1}}=\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} 10.8333 & -12.6667 & 1.83333 \\ -12.6667 & 29.3333 & -16.6667 \\  1.83333 & -16.6667  & 14.83333 \end{bmatrix}\begin{bmatrix} 1 &0 &0  &0  &0 \\ 0 &1  &0  &0  &0 \\ 0 &0  &1  &0  &0 \end{bmatrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{K^{1}}=\begin{bmatrix} 10.8333 &-12.6667 &1.8333    &0  &0 \\ -12.667 &29.33333  &-16.6667  &0  &0 \\ 1.8333  & -16.6667 & 14.83333 &0  &0 \\ 0 &0  &0  &0  &0 \\ 0 &0  &0  &0  &0 \end{bmatrix}$$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.32)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{F^{1}}=\mathbf{L^{1^{T}}}.\mathbf{f^{1}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{F^{1}}=\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} 12\\ 0.416667\\ 0.208333 \end{bmatrix} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{F^{1}}=\begin{bmatrix} 12\\ 0.416667\\ 0.208333\\ 0\\ 0 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.33)
 * }

For Element 2
<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

Ist Basis Function: $$b_{1}^{(2)}=\frac{(x-x_{4}).(x-x_{5})}{(x_{3}-x_{4}).(x_{3}-x_{5})}=\frac{(x-0.75).(x-1)}{(0.5-0.75).(0.5-1)}=8x^{2}-14x+6 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.34)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IInd Basis Function: $$b_{2}^{(2)}=\frac{(x-x_{3}).(x-x_{5})}{(x_{4}-x_{3}).(x_{4}-x_{5})}=\frac{(x-0.5).(x-1)}{(0.75-0.5).(0.75-1)}=-16x^{2}+24x-8 $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.35)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

IIIrd Basis Function : $$b_{3}^{(2)}=\frac{(x-x_{3}).(x-x_{4})}{(x_{5}-x_{3}).(x_{5}-x_{4})}=\frac{(x-0.5).(x-0.75)}{(1-0.5).(1-0.75)}=8x^{2}-10x+3 $$

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.36)
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$${b_{1}^{(2)}}'=16x-14, \quad\quad {b_{2}^{(2)}}'=-32x+24 , \quad\quad {b_{3}^{(2)}}'=16x-10 $$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.37)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

As we know from the definition of K matrix <span id="(1)">
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{K_{ij}^{(2)}}=\int_{\alpha }^{\beta }{b_{i}^{(2)}}'.a_{2}.{b_{j}^{(2)}}'dx $$
 * style="width:95%" |
 * style="width:95%" |

$$ <span id="(1)">
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.38)
 * }
 * {| style="width:100%" border="0"

So, the elements of local k matrix (for Element 2) are:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{11}^{(2)}}=\int_{0.5 }^{1 }(16x-14).(2+3x).(16x-14)dx=17.8333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.39)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{k_{12}^{(2)}}=\mathbf{k_{21}^{(2)}}=\int_{0.5 }^{1}(16x-14).(2+3x).(-32x+24)dx=-20.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.40)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{13}^{(2)}}=\mathbf{k_{31}^{(2)}}=\int_{0.5 }^{1}(16x-14).(2+3x).(16x-10)dx=2.83333 $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.41)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{k_{22}^{(2)}}=\int_{0.5 }^{1}(-32x+24).(2+3x).(-32x+24)dx=45.3333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.42)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{23}^{(2)}}=\int_{0.5 }^{1 }(-32x+24).(2+3x).(16x-10)dx=-24.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.42)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{k_{33}^{(2)}}=\int_{0.5 }^{1}(16x-10).(2+3x).(16x-10)dx=21.83333 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.43)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

S0, the Element Stiffness Matrix or Local Stiffness Matrix for the Element 2 is :
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{k^{(2)}}=\begin{bmatrix} 17.8333 & -20.6667 & 2.83333 \\ -20.6667 & 45.3333 & -24.6667 \\  2.83333 & -24.6667  & 21.83333 \end{bmatrix}

$$ $$ 2
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.44)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

As we know from the definition of weak form the Force Matrix is given by:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{f_i} = b_i (\alpha )h + \int_{0.5}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.45)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the elements of force matrix are: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ f_{1}^{(2)}=\int_{0.5}^{1}(8x^{2}-14x+6).(5x)dx=0.28333 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{2}^{(2)}=\int_{0.5}^{1}(-16x^{2}+24x-8).(5x)dx=1.25 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_{3}^{(2)}=\int_{0.5}^{1}(8x^{2}-10x+3).(5x)dx=0.416667 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So, the Element Force Matrix for Element 2 is:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{f^{(2)}}=\begin{bmatrix} 0.208333\\ 1.25\\ 0.416667 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.46)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The transformation from Local to Global is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

Global K,F matrices For Element 2
So, the global K and F matrices for Element 2 are given as:
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{(e)}}= \mathbf{L^{(e)^{T}}}.\mathbf{k^{(e)}}.\mathbf{L^{(e)}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.47)
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{F^{(e)}}=\mathbf{L^{(e)^{T}}}.\mathbf{f^{(e)}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.48)
 * }
 * {| style="width:100%" border="0"

For Element 2 $$ \mathbf{L^{(e)}}=\mathbf{L^{(2)}} $$.
 * style="width:95%" |
 * style="width:95%" |
 * }

$$ \mathbf{L^{(e)}} $$ For Element 2
<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \mathbf{L^{(2)}}=\begin{bmatrix} 0 &0 &1  &0  &0 \\ 0 &0  &0  &1  &0 \\ 0 &0  &0  &0  &1 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.49)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{K^{(2)}}= \mathbf{L^{(2)^{T}}}.\mathbf{k^{(2)}}.\mathbf{L^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{K^{(2)}}=\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 17.8333 & -20.6667 & 2.83333 \\ -20.6667 & 45.3333 & -24.6667 \\  2.83333 & -24.6667  & 21.83333 \end{bmatrix}\begin{bmatrix} 0 &0 &1  &0  &0 \\ 0 &0  &0  &1  &0 \\ 0 &0  &0  &0  &1 \end{bmatrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{K^{(2)}}=\begin{bmatrix} 0 &0 &0   &0  &0 \\ 0 &0  &0   &0  &0 \\ 0 &0  & 17.83333 &-20.6667  &2.83333 \\ 0 &0  &-20.6667  &45.3333  &-24.6667 \\ 0 &0  &2.83333  &-24.6667  &21.83333 \end{bmatrix}$$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.50)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{F^{(2)}}=\mathbf{L^{(2)^{T}}}.\mathbf{f^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\mathbf{F^{(2)}}=\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 0.208333\\ 1.25\\ 0.416667 \end{bmatrix} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |

$$ \Rightarrow \mathbf{F^{(2)}}=\begin{bmatrix} 0\\ 0\\ 0.208333\\ 1.25\\ 0.416667 \end{bmatrix} $$ $$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.51)
 * }

Global K matrix for 2 elements
<span id="(1)">
 * {| style="width:100%" border="0"

The Global K matrix for 2 elements is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{K}}=\sum_{e=1}^{2}\mathbf{K^{e}} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.52)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow \mathbf{\tilde{K}}=\mathbf{\tilde{K}^{(1)}}+\mathbf{\tilde{K}^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow \tilde{K}=\begin{bmatrix} 10.8333 &-12.6667 &1.8333    &0  &0 \\ -12.667 &29.33333  &-16.6667  &0  &0 \\ 1.83333 &-16.6667  & 32.6666 &-20.6667  &2.83333 \\ 0 &0  &-20.6667  &45.3333  &-24.6667 \\ 0 &0  &2.83333  &-24.6667  &21.83333 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.53)
 * }

Global F matrix for 2 elements
<span id="(1)">
 * {| style="width:100%" border="0"

The Global F matrix for 2 elements is given by: <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{F}}=\sum_{e=1}^{2}\mathbf{F^{e}} $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.54)
 * }
 * {| style="width:100%" border="0"

$$ \Rightarrow \mathbf{\tilde{F}}=\mathbf{\tilde{F}^{(1)}}+\mathbf{\tilde{F}^{(2)}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow \tilde{F}=\begin{bmatrix} 12\\ 0.416667\\ 0.416667\\ 1.25\\ 0.41667 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.55)
 * }

Global d matrix for 2 elements
<span id="(1)">
 * {| style="width:100%" border="0"

As we know, <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{K}}.\mathbf{\tilde{d}}=\mathbf{\tilde{F}} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$\Rightarrow \mathbf{\tilde{d}}=\mathbf{\tilde{K}^{-1}}.\mathbf{\tilde{F}} $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

So. the $$ \mathbf{\tilde{d}} $$ is :
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \mathbf{\tilde{d}}=\begin{bmatrix} 7.8642\\ 6.5882\\ 5.5934\\ 4.7540\\ 4 \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.56)
 * }

At $$ x = 0.5 \quad $$ ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h(0.5) = 5.5933908733 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.57)
 * }

The exact solution can be found from the eqn. 6.2.14. It is given by ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u(0.5) = 5.59352829 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Hence the error we get for $$ no. \quad of \quad elements \quad (n) = 2 $$ is ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$u^h(0.5) - u(0.5) = 0.00013827 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.2.58)
 * }

MATLAB code with output
The K,d and F matrix for 2 elements are as follows:

K =

10.8333 -12.6667    1.8333         0         0  -12.6667   29.3333  -16.6667         0         0    1.8333  -16.6667   32.6667  -20.6667    2.8333         0         0  -20.6667   45.3333  -24.6667         0         0    2.8333  -24.6667   21.8333

F =

[12 0.416666666666667 0.416666666666667 1.25000000000000 0.416666666666667]^T

d =

[7.8642  6.5882   5.5934   4.7540   4]^T

The outputs for two elements are as follows:

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=2 = 5.593386e+000

Error at n=2 is 1.380272e-004

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=4 = 5.593514e+000

Error at n=4 is 9.391941e-006

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=6 = 5.593522e+000

Error at n=6 is 1.871865e-006

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=8 = 5.593523e+000

Error at n=8 is 5.817644e-007

u,u^h vs x Plot
--Eml5526.s11.team5.srv 19:24, 6 April 2011 (UTC)

=6.3 Beam Tutorial for Calculix by Ty Beede=

Given: Tutorial by Beede on mechanicalhacks blog
A step by step tutorial explaining how to analyze a beam using Calculix can be found at

Find: Reproduce Beam Tutorial
Reproduce all steps in the tutorial by Ty Beede

Solution
Calculix is an open source Finite Element code (with ABAQUS like input) available at http://www.dhondt.de/

More details on installation can be found at Calculix Installation Procedure

Analysis using Calculix is divided into 4 steps:


 * 1) Buliding Geometry and Meshing


 * 2) Exporting Mesh, Loads, and Boundary Conditions


 * 3) Writing an Input File for the CCX Solver


 * 4) Post Processing Results in CGX

In this tutorial we analyze a cantilever beam subjected to pressure load at one end and fixed at the other end.

Step 1: Buliding Geometry and Meshing
Create a new file 'Beam.fbd' in the SciTE editor. To construct the 3D geometry of the beam we do the following:

(Note: code used to generate the geometry and commands required to display the geometry at each step are embedded at the bottom of each screen shot)

a) Define points and join them to create lines

b) Sweep the lines to create a surface

c) Sweep the surface to create a volume

d) Assign Element Types and mesh the geometry

Step 2: Exporting Mesh, Loads, and Boundary Conditions
In the earlier section we created the geometry and generated a mesh.

The 'send beam abq' command used at the end of the 'beam.fbd' file, generates the node & element information upon preprocessing. This data is saved in 'beam.msh' file in the same directory.

The mesh data can be found here:beam.msh

After the mesh has been created we need to Apply the Boundary conditions and Loads:

Applying the Boundary Conditions
-> Use the command 'plot n all' to display all the nodes

-> Zoom in to the left end of the beam where Boundary condition is to be applied.

-> Type 'qadd fixed' and hit enter. The qadd command provides a graphical selection method that can be used interactively with the mouse. The qadd command will stay active until the 'q' button is pressed on the keyboard.

At the tip of the mouse pointer is a little rectangle outlined in black. With the 'qadd' command active the mouse is used to select items like nodes, elements, faces, points, lines, and etc. The items desired for selection must fall within the bounds of the little black rectangle shown above. When the rectangle is hovering over the desired object press the keyboard button corresponding to the item of interest in order to add it to the set. In this case the 'n' button will be used to select nodes within the bounds of the rectangle.

For adding items in bulk it will be necessary to change the size of the selection rectangle.The 'r' key on the keyboard is used to define two points of a rectangle which will become the new selection area for qadd. Simply press the 'r' button and the current cursor position becomes the new location for a corner of the selection area. Move the mouse in the horizontal and vertical directions such that it is offset from the previous point. Press the 'r' button again and a new selection rectangle has been defined.

In order to select all the nodes of the left hand edge the selection area will be made large enough to encompass the entire edge. In addition there are two selection modes used by the qadd command. The default selection mode adds one item each time the corresponding keyboard button is pressed. This works good for selecting one of a particular item with a small selection rectangle. If the desire is to bulk add all the items located within the selection rectangle the mode:a should be used. While the 'qadd' command is active press the 'a' button on the keyboard to active bulk add.

-> after mode:a has been selected hit 'n' to select all the nodes inside the rectangle to apply the fixed boundary condition.

The output window looks like this upon selecting the nodes to be fixed:

Now that the nodes have been added to the set fixed the qadd command can be quit using the q button.

-> In order to visualize the nodes contained in the set fixed, Change the color of the nodes and toggle the background color

Background can be toggled to black using the menu system command 'Viewing->Toggle Background Color'

Color of the node list can be changed by using the following command 'plus n fixed g'

The visualization screens appears as follows:

-> The boundary condition is saved to a separate file(.bou) using the command 'send fixed abq spc 123'

The generated file containing Boundary conditions, the command used to generate the file and the message displayed on the output window are shown:

The boundary conditions data can be found at: fixed_123.bou

Applying the loads
In this example we apply a distributed pressure load on an element towards the right end in the Y-direction, to create a moment along the Z-direction.

-> Plot the faces of the beam using 'plot f beam' command, then use 'view elem' command.

-> Zoom in to the right end and type 'qadd load' and hit enter. Then move the mouse over the element where load needs to be applied and hit 'f' to select that face.



-> Then type 'send load abq pres 10' and hit enter, this command applies a load of 10 on the selected face and saves the load data file.



Step 3: Writing an Input File for the CCX Solver
Open SciTE and create a new file named 'beam.inp'. This document should be located with the files from the previous articles. Add the following code:

In this example the input file contains the files names of the node, element data, boundary conditions, loads. and other information like material properties are defined in this file. (Note: It is not necessary that node data, boundary conditions, load data be declared in separate files. If the data is simple & small it can directly be typed in the input file)

After the input file is ready, we can submit the job for solving ('Ctrl+F10' or 'Tools>Solve').

a screen shot of the input code for ccx and the output window after solving the code are shown:



--Eml5526.s11.team5.srv 19:19, 6 April 2011 (UTC)

Step 4: Post Processing Results in CGX
->After the solving is done, post-processing can be done by clicking on 'Tools>Post Process' or 'Shift + F10'.

->In this simulation there are two fields available for visualization. They are displacement and stress. Clicking on one of the menu items at this level will select the results field of interest.Then the particular information from the field will need to be specified. This is done by choosing an option from within the -Entity- submenu will result in a colored plot of the selected field information as shown below.











-> The deformed shape can be visualized by selecting 'Viewing>ToggleAdd-Displacement'

To scale up the displacement use the keyboard command scal d 10000.



->Another option to visualize the shape of the deformed structure is through animation. This feature will animate the model’s deformation with user specified scale. Use the menu system and select Animate->Start to view the animation. The image below shows the beam as it oscillates between -100% and 100% of the deformation amplitude.



->The cgx post-processing mode allows the user to specify a cut plane. The plane will cut through the model and show the field results for the intersection of the model and plane. The cut plane is defined by three nodes. The nodes can be specified from the menu system using the Cut menu.



The following image shows the result of cutting the beam top to bottom. The field shown is the normal stress in the x direction. The results are as we would expect, the top of the beam is in tension and the bottom is in compression.



->In addition values from the different result fields can be plotted to a graph for interpretation. The graph command is accessible from the menu system under Graph. The graph is also accessible through the command line.

The nodes are selected three at a time, then the 'tra l 10' command is used to pan the model to the left by translation of 10 units. At which point another few nodes are selected. This process is repeated as the translation command moves towards the free end. Then the last nodes are selected and the right mouse button is used to finish the selection process. Finally, the following graph is promptly displayed.



->While viewing the colored plot of a results field it is also possible to add the nodal values. The following images shows top side of the beam’s free end. The nodal values are printed in yellow using the commmand 'plus nv all y'. These values agree with the results of the graph.



--Eml5526.s11.team5.vijay 20:53, 6 April 2011 (UTC)

= 6.4 Divergence Theorem =

Given
1) A vector field  $$ q_x =-y^2, q_y=-2xy \quad $$ on a domain

2) A vector field $$ q_x =3x^2y + y^3, q_y=3x + y^3 \quad $$ on a domain

3) $$ \oint_\Gamma \mathbf{n} d \Gamma =0 \quad $$

$$ \int_\Omega \vec \nabla \cdot \vec q d \Omega =\oint_\Gamma \vec q \cdot \vec n d \Gamma \quad $$ or  $$ \int_\Omega \mathbf{\nabla}^T  \mathbf{q} d \Omega=\oint_\Gamma \mathbf{q}^T \mathbf{n} d \Gamma \quad $$ <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.4.1) $$

Find
Verify the divergence theorem

Solution
1) To verify the divergence theorem, first evaluate the integrand on the left hand side of Eq. 6.4.1

Integrating the above over the problem domain gives,

Evaluating the boundary integral counterclockwise gives

Thus, the divergence theorem is verified.

2) To verify the divergence theorem, first evaluate the integrand on the left hand side of Eq. 6.4.1

Integrating the above over the problem domain gives,

Evaluating the boundary integral counterclockwise gives

Eq. 6.4.10 equals Eq. 6.4.8 Thus the divergence theory is verified.

3) To prove using the divergence theorem, first let, Using the divergence theorem,

And since,

Then,

Thus,

--Eml5526.s11.team5.smith 15:43, 6 April 2011 (UTC)

= 6.5 =

Given: 2D Strong Form and BCs
<span id="(1)">
 * {| style="width:100%" border="0"

$$\Omega=\bar{w}=\Box$$ biunit square $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.5.1)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$div(\underline{K} \cdot grad (u))+f=\rho c \frac{\partial u}{\partial t}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.5.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \underline{K}= \underline{I}$$ (identity matrix) $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.5.3)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$f=0 \quad$$,  $$\frac{\partial u}{\partial t}= 0$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.5.4)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Essential boundary condition:
 * style="width:95%" |
 * style="width:95%" |
 * $$g=2 \quad$$ on  $$\partial\Omega \quad$$

$$
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 6.5.5)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Natural boundary condition:
 * style="width:95%" |
 * style="width:95%" |
 * none

$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

Find: Temperature inside and error at (0,0)
<span id="(1)">
 * {| style="width:100%" border="0"

Use 2D Lagrange interpolate basis functions (LIBF) with n=m=2,4,6,8,... to solve for $$u^h$$ until error is less then $$10^{-6}$$ at center (x,y)=(0,0) $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Solution
$$ K_{ij}^e = \int_{\Omega^e} \triangledown{b^e_i}'(x, y) \kappa(x, y) \triangledown{b^e_j}'(x, y) d \omega_x \ $$

, where $$ x = x^e \ $$, where $$x^e \ $$ represents local node.

For linear 2D case with boundary being square

Bases, $$ b_i \ $$ are defined to be $$ L_{[I,J]} (x, y) = N^e_{[I,J]} (x, y) \ $$

where $$ N^e_{[I,J]}(x, y) = N_I^e(x)N_J^e(y) \ $$ and $$ I, J \ $$ corresponds to numbering of nodes on element scale.



also noting that it's a rectangular shape, such that the distance from x: 1-to-2 = 1-to-3 and node 1 = node 4 in 'x' coordiantes, y: 1-to-4 = 1-to-3 and node 1 = node 2 in 'x' coordinates, we substitute $$ y^e_4 $$ for $$ y^e_3 \ $$ and $$ y^e_1 $$ for $$ y^e_2 \ $$ and the same for other coordinates

$$ N^e_{[1,1]} (x, y) = \frac{(x^e-x^e_2)(x^e-x^e_2)(x^e-x^e_1)}{(x^e_1 - x^e_2)(x^e_1 - x^e_2)} \frac{(y^e-y^e_1)(y^e-y^e_4)(y^e-y^e_4)}{(y^e_1 - y^e_4)(y^e_1 - y^e_4)} = \frac{1}{(A^e)^2}(x-x^e_2)^2(x-x^e_1)(y - y^e_4)^2(y - y^e_1) = b_1^e (x,y) \ $$

$$ N^e_{[2,1]} (x, y) = \frac{(x^e-x^e_1)(x^e-x^e_2)(x^e-x^e_1)}{(x^e_2 - x^e_1)(x^e_2 - x^e_1)} \frac{(y^e-y^e_1)(y^e-y^e_4)(y^e-y^e_4)}{(y^e_1 - y^e_4)(y^e_1 - y^e_4)} = \frac{1}{(A^e)^2}(x-x^e_1)^2(x-x^e_2)(y - y^e_4)^2(y - y^e_1) = b_1^e (x,y) \ $$

$$ N^e_{[2,2]} (x, y) = \frac{(x^e-x^e_1)(x^e-x^e_2)(x^e-x^e_1)}{(x^e_2 - x^e_1)(x^e_2 - x^e_1)} \frac{(y^e-y^e_1)(y^e-y^e_4)(y^e-y^e_4)}{(y^e_1 - y^e_4)(y^e_1 - y^e_4)} = \frac{1}{(A^e)^2}(x-x^e_1)^2(x-x^e_2)(y - y^e_4)^2(y - y^e_1) = b_1^e (x,y) \ $$

$$  N^e_{[1,1]} (x, y) = \frac{(x^e-x^e_2)(x^e-x^e_2)(x^e-x^e_1)}{(x^e_1 - x^e_2)(x^e_1 - x^e_2)} \frac{(y^e-y^e_1)(y^e-y^e_4)(y^e-y^e_4)}{(y^e_1 - y^e_4)(y^e_1 - y^e_4)} = \frac{1}{(A^e)^2}(x-x^e_1)^2(x-x^e_2)(y - y^e_4)^2(y - y^e_1) = b_1^e (x,y) \ $$

Next we will convert the square element to same size parent element of bi-unit square

$$ x^e = x_1^e N_1^e(\xi) + x_2^e N_2^e(\xi) = x_1^e \frac{1-\xi}{2} + x_2^e \frac{1+\xi}{2}, \xi \in [-1 +1] \ $$

$$ y^e = y_1^e N_1^e(\eta) + y_4^e N_4^e(\eta) = y_1^e \frac{1-\eta}{2} + y_4^e \frac{1+\eta}{2}, \eta \in [-1 +1] \ $$

Then by noting which side/value of 'x' will be greater in defined rectangular shown. Moreover, since we will be using 4 square elements (in order to have at least 1 element that is not forced by essential boundary condition, we had to partition the Global element into 4 sub-elements).

Therefore, we convert element coordinates to parent coordinates

$$ N^e_{[1,1]} (\xi, \eta) = \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} + x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right)^2 \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right)^2 \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2}\right)= b^e_1(\xi, \eta) \ $$

$$ N^e_{[2,1]} (\xi, \eta) = \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} + x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right) \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)^2\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right)^2 \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2}\right)= b^e_1(\xi, \eta) \ $$

$$ N^e_{[2,2]} (\xi, \eta) = \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} + x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right) \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)^2\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right) \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2})^2\right)= b^e_1(\xi, \eta) \ $$

$$ N^e_{[1,2]} (\xi, \eta) = \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} + x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right)^2 \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right) \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2})^2 \right)= b^e_1(\xi, \eta) \ $$

Jacobian Matrix:

$$ \mathbf{J}^e =

\begin{bmatrix} \frac{\partial x}{\partial \xi} &  \frac{\partial y}{\partial \xi}    \\ \\ \frac{\partial x}{\partial \eta} &  \frac{\partial y}{\partial \eta} \end{bmatrix}

$$

taking the $$ x^e \ $$ and $$ y^e \ $$ values $$ \mathbf{J}^e =

\begin{bmatrix} \frac{\partial x}{\partial \xi} &  \frac{\partial y}{\partial \xi}    \\ \\ \frac{\partial x}{\partial \eta} &  \frac{\partial y}{\partial \eta} \end{bmatrix} =

\begin{bmatrix} \frac{1}{2} (x_2^e - x_1^e) & 0    \\ \\ 0 &  \frac{1}{2} (y_4^e - y_1^e) \end{bmatrix}

$$

From given

$$ d \omega_x = det(\mathbf{J}^e) d \omega_\xi = d \omega_x = \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \ $$

or

$$ d \omega_\xi = \frac{ d \omega_x}{det(\mathbf{J}^e)} = d \omega_\xi = \frac{ d \omega_x}{\frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e)} $$

With this last information we have a completely defined stiffness element matrix $$ K^e \ $$. Also note that since we have, 4 degrees of freedom in 1 element node, we do expect the gram matrix $$ \Gamma \ $$ to be 4x4, therefore, $$ \kappa_{4x4} \ $$

$$ K_{ij}^e = \int_{\Omega^e} \triangledown{b^e_i}'(x, y) \kappa(x, y) \triangledown{b^e_j}'(x, y) d \omega_x \ $$

$$ K_{ij}^e = \int_{-1}^{+1}

\begin{bmatrix} \triangledown{b^e_1} \triangledown{b^e_1} & \triangledown{b^e_1} \triangledown{b^e_2} & \triangledown{b^e_1} \triangledown{b^e_3} & \triangledown{b^e_1} \triangledown{b^e_4}      \\ & \triangledown{b^e_2} \triangledown{b^e_2} & \triangledown{b^e_2} \triangledown{b^e_3} & \triangledown{b^e_2}\triangledown{b^e_4}      \\ & &  \triangledown{b^e_3} \triangledown{b^e_3} & \triangledown{b^e_3} \triangledown{b^e_4}      \\ symetry & &   & \triangledown{b^e_4} \triangledown{b^e_4} \end{bmatrix}

\begin{bmatrix} 1 & 0 & 0 & 0      \\ 0 & 1 &  0 & 0      \\ 0 & 0 &  1 & 0      \\ 0 & 0 &  0 & 1      \end{bmatrix}

\left ( d \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \right ) \ $$

note, to save some space in typing the formula (else it stretched out of the page), the independent variables for $$ b_i^e (\xi, \eta) \ $$ were avoided

for the initial square of length of '1', we divided into 4 smaller squares each 0.5 length, the element nodes stiffness matrix came out to be:

K(:,:,1) = [  9/286720, 0, 0, -3/1638400] [          0, 0, 0,          0] [          0, 0, 0,          0] [ -3/1638400, 0, 0,   9/286720] K(:,:,2) = [   9/286720,  31/45875200,   1/13107200,  -3/1638400] [ 31/45875200, 129/16056320,  -5/19267584,  1/13107200] [  1/13107200,  -5/19267584, 129/16056320, 31/45875200] [  -3/1638400,   1/13107200,  31/45875200,    9/286720]

this answer shows 2 elements stiffness matrices located along 'x' axes. Noting that because of the symmetry of each elements, the equivalent spacing of each node in elements, and because the material matrix $$ \kappa \ $$ is identity matrix, all stiffness element matrices will be the same. The only difference in element stiffness matrix will occur when we will change the number of stiffness.

,where

$$ b^e_1(\xi, \eta) = \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} +  x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right)^2 \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right)^2 \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2}\right)\ $$

$$ b_2^e(\xi, \eta) = \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} +  x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right) \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)^2\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right)^2 \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2}\right)\ $$

$$ b_3^e (\xi, \eta)  =  \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} +  x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right) \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)^2\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right) \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2})^2\right) \ $$

$$  b_4^e (\xi, \eta) =   \frac{1}{(A^e)^2}\left(x_1^e \frac{1-\xi}{2} +  x^e_2\left[ \frac{1+\xi}{2}- 1 \right]\right)^2 \left(x_1^e \left[\frac{1-\xi}{2}-1\right] +  x^e_2\frac{1+\xi}{2}\right)\left(y_1^e \frac{1-\eta}{2} + y^e_4\left[ \frac{1+\eta}{2} - 1 \right] \right) \left(y_1^e \left[ \frac{1+\eta}{2} - 1 \right] + y^e_4\frac{1+\eta}{2})^2 \right)  \ $$

force
From prof. Vu Quoc's notes for

$$ f^e(w) = - \int_{\Gamma_h} wh d {\Gamma_h} + \int_\omega wf d \omega \ $$

$$ f^e(w) = - \int_{\Gamma_h} wh d {\Gamma_h} + \int_\omega w0 d \omega = - \int_{\Gamma_h} wh d {\Gamma_h} \ $$

converting to local coordinates

$$ f^e(w) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

F(:,:,1) = [ 1/1024, 1/1024] [      0,       0] [      0,       0] [ 1/1024, -1/1024] F(:,:,2) = [ 1/1024, 1/1024] [ 1/2048, -1/2560] [ 1/2048,  1/2560] [ 1/1024, -1/1024]

KGlobal = 9/4480        0              0             -3/25600        0              0              0              0                0              0              9/4480         0              0             -3/25600        0              0              0              0              0              0              0              0              0              0              0              0              0             -3/25600        0              0              9/2240         0              0             -3/25600        0              0              0             -3/25600        0              0              9/2240         0              0             -3/25600        0              0              0              0              0              0              0              0              0              0              0              0              0             -3/25600        0              0              9/4480         0              0              0              0              0              0             -3/25600        0              0              9/4480         0              0              0              0              0              0              0              0              0              0

Kchange =

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              0              0              0              0              0              0              0              0              0              0              0              0              9/2240         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              0              0              0              0              0              0              0              0              0              0              0              0

And the answer came out to be '2'

u = 2

We run few more basis function increments and we got an answer of '2' as well.

Discussion
Since the element was bounded by essential BCs of the same temperature, namely '2' and no heat generation inside the element, we would expect the temperature, 'u', to be uniform and the same as given BCs. The results from increasing bases function degree, that should have had increased the accuracy on convergence, still output as '2'. There was no more accurate results available.

Eml5526.s11.team5.savery


 * }

Discussion
This answer is as it was expected, since no heat was generated and the value at (0,0) also should stay on 2 (essential BC). This was exactly what was demonstrated in our HW. We got consistent answers of '2' inside the nodes of the mesh.

--Eml5526.s11.team5.JA 20:51, 6 April 2011 (UTC)

Solution to (3): WolframAlpha Syntax Explained
<span id="(1)">
 * {| style="width:100%" border="0"

As seen in Meeting 36, page 1, some interesting discrepancies exist within the syntax of WolframAlpha. The following will explain those discrepancies. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

According to, the Mathematica Guide for Mathematica syntax as part of the Wolfram Mathematica Documentation Center, "Basic Syntax" explains that any function argument is to be augmented with square brackets; i.e. $$\boldsymbol{f}[x,y]$$. For example, the function "transpose" would call for square brackets around the data desired to be transposed. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Also of note, according to basic syntax, curly brackets are used only when defining a list of elements, or a series of elements. If the function "transpose" is used augmented with curly brackets, the data set inside those curly brackets will first be transformed into a series of elements; then transposed. This will yield an incorrect answer if attempting to transpose a matrix and accordingly in matrix algebra this can prove problematic. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

For instance, Mathematica syntax defines that $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\left \{ \left \{ 1,1,1,1 \right \},\left \{ 2,-1,3,2 \right \},\left \{ 3,2,6,3 \right \} \right \}$$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.25)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

will define an appropriate matrix: a series of 3 elements (each in 3 sepparate rows), with each of those 3 rows containing 4 elements of its own (each in 4 sepparate columns). The columns are sepparated by commas. Seen in this link is an example which results in the following: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\begin{bmatrix} 1 &1 &1  &1 \\ 2 &-1  &3  &2 \\ 3 &2  &6  &3 \end{bmatrix}$$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.26)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Now, if an additional set of curly brackets are made to surround the previous matrix (as defined in syntax), further bounded by a function such as "transpose" (seen below) $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

transpose \left \{\left \{\left \{ 1,1,1,1 \right \},\left \{ 2,-1,3,2 \right \},\left \{ 3,2,6,3 \right \} \right \}\right \}$$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.27)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

the matrix before executing the function will be reformatted to account of the additional curly brackets which call for the data to be set in series. This is an erroneous and unnecessary step, as it redefines the matrix as a series itself of its own rows. The transpose of this outcome results in the following and can be seen in the following link : $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

(\left \{ 1,1,1,1 \right \},\left \{ 2,-1,3,2 \right \},\left \{ 3,2,6,3 \right \})^{T}=\begin{pmatrix} \left \{ 1,1,1,1 \right \}\\ \left \{ 2,-1,3,2 \right \}\\ \left \{ 3,2,6,3 \right \} \end{pmatrix}$$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.28)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Obviously, this is not the outcome being sought. If square brackets are used instead, the superfluously created series set would be avoided. In fact, in the case of Equation 6.7.25, additional brackets would not even be required as utilizing the command "transpose" with Equation 6.7.25 is acceptable because the matrix is defined sufficiently already. Brackets here would be superfluous and in the case of using curly brackets, would conjure incorrect results as explained above. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

It should be noted that the function "invert" does however recognize the matrix syntax form correctly by utilizing curly brackets (square brackets also yield the same result for the function "invert"). The function "invert" automatically can recognize the correct form of a matrix without first performing any re-definition of the matrix itself. This allows the function "invert" to utilize either bracket type. However, for continuity and accuracy across all functions utilized within Mathematica syntax, it is recommended to adhere to the Basic Syntax and implement square brackets. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

--Eml5526.s11.team5.mcdaniel 20:45, 6 April 2011 (UTC)

>>gauss([1 0 1],0,4) ans = 25.3333
 * }

>> gauss([1 0 2 0 0],-1,1) ans = 1.7333

Given
and

Find
Use three-point Gauss quadrature to evaluate the following integrals. Compare to the analytical integral.

Solution
For 3 points the points and weights are:


 * $$x_1=0 \quad$$
 * $$x_2=\frac{\sqrt{15}}{5}$$
 * $$x_3=\frac{-\sqrt{15}}{5}$$
 * $$W_1=\frac{8}{9}$$
 * $$W_2=\frac{5}{9}$$
 * $$W_3=\frac{5}{9}$$

a)
 * Next we need to map the physical domain [a,b] to the parent domain [-1,1] where a=-1 b=-1 for this case:


 * $$\xi=\frac{1}{2}(a+b)+\frac{1}{2}x(b-a)=x$$


 * $$f(x)=\frac{x}{x^2+1} \quad$$


 * Now we can use the weights and points to obtain a solution to eqn 6.8.11:


 * $$I=\frac{L}{2}W_1f(x_1)+\frac{L}{2}W_2f(x_2)+\frac{L}{2}W_3f(x_3)$$


 * $$I=\frac{b-a}{2}W_1\left [\frac{x_1}{x_1^2+1}\right]+\frac{b-a}{2}W_2\left [\frac{x_2}{x_2^2+1}\right]+\frac{b-a}{2}W_3\left [\frac{x_3}{x_3^2+1}\right]$$


 * $$I=\frac{1--1}{2}*\frac{8}{9}\left [\frac{0}{0^2+1}\right]+\frac{1--1}{2}*\frac{5}{9}\left [\frac{\left(\frac{\sqrt{15}}{5}\right)}{\left(\frac{\sqrt{15}}{5}\right)^2+1}\right]+\frac{1--1}{2}*\frac{5}{9}\left [\frac{\left(\frac{-\sqrt{15}}{5}\right)}{\left(\frac{-\sqrt{15}}{5}\right)^2+1}\right]$$

$$I=0\quad$$


 * Analytical solution:

$$\int_{-1}^{1}\frac{\xi}{\xi^2+1}\, dx=\left[\frac{log\left(\xi^2+1\right)}{2}\right]_{-1}^{1}=\left(\frac{log\left(1^2+1\right)}{2}\right)-\left(\frac{log\left((-1)^2+1\right)}{2}\right)=0$$

b)
 * We need to map the physical domain [a,b] to the parent domain [-1,1] where a=-1 b=-1 for this case:


 * $$\xi=\frac{1}{2}(a+b)+\frac{1}{2}x(b-a)=x$$


 * $$f(x)=\cos^2(\pi x)\quad$$


 * Now we can use the weights and points to obtain a solution to eqn 6.8.12:


 * $$I=\frac{L}{2}W_1f(x_1)+\frac{L}{2}W_2f(x_2)+\frac{L}{2}W_3f(x_3)$$


 * $$I=\frac{b-a}{2}W_1\left [\cos^2(\pi x_1)\right]+\frac{b-a}{2}W_2\left [\cos^2(\pi x_2)\right]+\frac{b-a}{2}W_3\left [\cos^2(\pi x_3)\right]$$


 * $$I=\frac{1--1}{2}*\frac{8}{9}\left [\cos^2(\pi*0)\right]+\frac{1--1}{2}*\frac{5}{9}\left [\cos^2\left(\pi \left(\frac{\sqrt{15}}{5}\right)\right)\right]+\frac{1--1}{2}*\frac{5}{9}\left [\cos^2\left(\pi \left(\frac{-\sqrt{15}}{5}\right)\right)\right]$$

$$I=1.52996\quad$$


 * Analytical solution:

$$\int_{-1}^{1}\cos^2(\pi \xi)\, dx=\left[\frac{\sin{\pi\xi}\cos{\pi\xi}}{2\pi}\right]_{-1}^{1}=\left(\frac{\sin{\pi*1}\cos{\pi*1}}{2\pi}\right)-\left(\frac{\sin{\pi*-1}\cos{\pi*-1}}{2\pi}\right)=1$$


 * We can see that the analytical method varies from the computational method and that the computational is not exact.

Using the gauss.m matlab code developed above and setting the number of points to 3 we get the same answers as we did manually:

I = 0

I = 1.5300

Given
$$\int_{-1}^{1}3\xi^3+2\, d\xi$$

Find
The integral above can be integrated exactly using two-point Gauss quadrature. How is the accuracy affected if

a. one-point quadrature is employed;

b. three-point quadrature is employed.

Solution
Using the Matlab code developed above we first set N=2 to see what the exact answer is:

>> gauss([ 3 0 0 2],-1,1) N = 2 ans = 4.0000

We then set N=1 to see what one-point Gauss quadrature results in:

>> gauss([ 3 0 0 2],-1,1) N = 1 ans = 4

And then N=3 for three-point Gauss quadrature:

>> gauss([ 3 0 0 2],-1,1) N = 3 ans = 4.0000

As you can see for this integral we get the same exact answer for all 3 so the accuracy is not affected.

Eml5526.s11.team5.savery 20:53, 6 April 2011 (UTC)

= 6.9 Heat conduction in 2D =

Example 8.3
Consider the heat conduction problem given in example 8.1 modeled with 16 quadrilateral elements. Solve this problem using finite element coding given in Fish and Belyschko section 12.5.

Given
(i) The conductivity of the plate is isotropic and $$ k = 5 W^0 C^{-1} \quad $$

(ii) Essential boundary condition T = 0 along the bottom edge and left edge.

(iii) Natural boundary condition $$ \bar{q} = 0 \ W m^{-1} $$ and  $$ \bar{q} = 20 \ W m^{-1} $$ along right edge and top edge respectively.

(iv) A constant heat source $$ s = 6 \ W m^{-2} $$

Find
Reproduce the post processing results obtained from the finite element code goven in section 12.5 in Fish and Belyschko book.

Solution
Using the matlab codes provided in section 12.5 of the book, we obtain the following results.

The Matlab code to analyze the problem 8.3 can be downloaded from here:Matlab code for Problem 8.3

This code can be altered to solve 2D heat conduction problems. We need to modify the 'input_file_xxxx.m' according to the boundary conditions, mesh size and the 'mesh2d.m' code also needs to be changed according to the given geometry.

For each problem, mesh configuration we create a separate 'input_file' which is called from 'preprocessor.m'.

The input files for each case are given below.

Given


(i) The conductivity of the concrete (inner material) is $$ k = 2 W^0 C^{-1} \quad $$ and that of the brick (outer material) is $$ k = 0.9 W^0 C^{-1} \quad $$

(ii) Essential boundary condition T = 10 along the top edge and T = 140 along the bottom edge.

(iii) Natural boundary condition $$ \bar{q} = 0 \ W m^{-1} $$ and  $$ \bar{q} = 20 \ W m^{-1} $$ along right edge and left edge respectively.

Find
Find the temperature distribution for different mesh densities:

(i) 2 x 2 quadrilateral element

(ii)4 x 4 quadrilateral element

(iii)8 x 8 quadrilateral element

Solution
Since the geometry has changed in this problem, we need to edit the mesh2d code to create a proper mesh. The code used for mesh creation for this geometry is shown below: % Mesh2D for 1 element

Matlab code for heat conduction in 2D plate using 64 elements
--Eml5526.s11.team5.vijay 20:55, 6 April 2011 (UTC)

=6.10 Newton's Law of Cooling=

Given
$$ \vec \nabla \cdot \vec q - s=0   on ~  \Omega ,\quad $$ <p style="text-align:right"> $$ \displaystyle    (Eq. 6.10.1) $$ $$ q_n=\vec q \cdot \vec n= \bar q   on   \Gamma_ {q} ,\quad $$ <p style="text-align:right"> $$ \displaystyle    (Eq. 6.10.2) $$ $$ q_n=\vec q \cdot \vec n= h(T-T_\infty )  on   \Gamma_h ,\quad $$ <p style="text-align:right"> $$ \displaystyle    (Eq. 6.10.3) $$ $$ T=\bar T  on   \Gamma_T ,\quad $$ <p style="text-align:right"> $$ \displaystyle    (Eq. 6.10.4) $$ $$ \Gamma= \Gamma_T \cup \Gamma_q \cup \Gamma_h,\quad $$ <p style="text-align:right"> $$ \displaystyle    (Eq. 6.10.5) $$

Find
1) The weak form for heat conduction in 2-D with boundary convection 2) Develop semi-discrete equations (ODEs). Give detailed expression of all quantities. 3) Solve G2DM2.0/D1 using 2D LIBF till $$\displaystyle 10^{-6} $$ accuracy at the center. Plot solution in 3D with contour.

Solution
1)Multiply energy balance equation (Eqn 6.10.1) and the natural boundary conditions (Eqn 6.10.2 and Eqn 6.10.3) by the weight function w

Applying Green's formula to the first term in Eqn. 6.10.6 yields,

Inserting Eq 6.10.9 into 6.10.6 yields

Where the first integral on the RHS of Eq. 6.10.10 was subdivided into the prescribed temperature and flux boundaries, which is possible because of Eq. 6.10.5. Substituting Eqs. 6.10.6 and 6.10.7 into the integral on $$ \Gamma_q $$ yields,

Now set w=0 on the essential boundary so the weight function will vanish on that portion of the boundary. Therefore, the integral on $$ \Gamma_T $$ vanishes and the weak form is given by

2) Let $$ u=[ N]{{u}_{e}}\quad $$ and $$ \quad w=[ N ]{{w}_{e}},$$ where $$\displaystyle [N] $$ is the Lagrange polynomial and $${u}_{e}\quad and \quad {w}_{e} $$ are the nodal values of $$ u \quad $$ and $$w \quad $$ respectively.

Then take the derivative of basis functions with respect to x1 and x2

The semi-discrete form becomes,

The heat source vector, conductivity matrix, and the capacitance matrix are

$$\tilde{f}(w)= -\int\limits_{h{{[N]}^{T}}d{{\Gamma }_{h}}}+\int\limits_{H{{[N]}^{T}}{{u}_{\infty }}d{{\Gamma }_{h}}}+\int\limits_{\Omega }{{{[N]}^{T}}f(x,t)d\Omega } $$

$$\tilde{k}(w,u)= \int\limits_{\Omega }{{{[B]}^{T}}[B]d\Omega }+\int\limits_{H{{[N]}^{T}}[N]d{{\Gamma }_{h}}} $$

$$\tilde{m}(w,u)= \int\limits_{\Omega }{\rho c{{[N]}^{T}}[N]\frac{\partial u}{\partial t}d\Omega } $$

--Eml5526.s11.team5.smith 15:46, 6 April 2011 (UTC)

= References =

<references