User talk:Eml5526.s11.team01.fercas

Hi team 01 =Problem 6.1: Quadratic Lagrangian Element Basis Function (QLEBF)= Solving a ordinary differential equation using Quadratic Lagrangian Element Basis Function (QLEBF). From Mtg.31

Given
The partial differential equation is:

The boundary conditions are:

$$ u\left(0\right) = 4$$

$$\frac{du}{dx}\left(1\right) = \frac{12}{5} $$

Find

 * 1. For nel = 2, compute $$\tilde{K}=\sum_{1}^{2}\tilde{K}^e$$ with $$\tilde{K}^e$$ given by equation (1) p.30-5
 * 2. Compute $$\tilde{k}^e$$ and $$L^e$$ for e = 1,2


 * 3. Compute $$\tilde{K}^e=L^{eT}k^eL^e$$ for e = 1,2. Compare to 1)
 * 4. Plot all the QLEBF for nel = 3
 * 5. Plot the exact solution, the numerical solution, and the error as a function of the number of elements

For nel = 2, compute $$\tilde{K}=\sum_{1}^{2}\tilde{K}^e$$ with $$\tilde{K}^e$$ given by equation (1) p.30-5
To calculate the stiffness matrix of the element we need the weak form; but, the weak form was already obtained in home work 5 which can be reach at http://en.wikiversity.org/wiki/User:Eml5526.s11.team01/Hwk5 Then, here we only copy that equation.

Now we use QLEBF to interpolate the solution $$U^h$$ and the weight funtion w. The aproximated soution can be interpolated as:


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

$$U^h=\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}u_1+\frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)}u_2+\frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}u_3 $$
 * $$\displaystyle (6.1.2) $$
 * }
 * }

Where $$u_{1}, u_{2}, u_{3}$$ are the nodal values of the funtion at the nodes 1, 2, and 3 respectively. The last equation can be written as:


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

$$U^h=\left [ N \right ]\left [ u \right ]^{e} $$
 * $$\displaystyle (6.1.3) $$
 * }
 * }

where $$\left [ u \right ]^{e}=\begin{bmatrix} u_1\\ u_2\\

u_3\end{bmatrix} $$

In the same way we can interpolate the weight funtion w to obtain:


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

$$w=\left [ N \right ]\left [ w \right ]^{e} $$
 * $$\displaystyle (6.1.4) $$
 * }
 * }

where $$\left [ w \right ]^{e}$$ is the vector with the nodal values of w

The derivative of $$U^h$$ and w are:


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

$$\frac{dU^h}{dx}=\begin{bmatrix} \frac{2x-x_2-x_3}{\left (x_1-x_2 \right )\left (  x_1-x_3\right )} & \frac{2x-x_1-x_3}{\left (x_2-x_1  \right )\left (  x_2-x_3\right )}  & \frac{2x-x_1-x_2}{\left (x_3-x_1  \right )\left (  x_3-x_2\right )} \end{bmatrix}\left [ u \right ]^e $$ The latter equation can be written as:
 * $$\displaystyle (6.1.5) $$
 * }
 * }


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

$$\frac{dU^{h}}{dx}=\left[B\right]\left [ u \right ]^{e} $$
 * $$\displaystyle (6.1.6) $$
 * }
 * }

and in similar way w can be written as:


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

$$\frac{dw}{dx}=\left[B\right]\left [ w \right ]^{e} =\left [ B \right ]^T\left [ w \right ]^{eT}$$ Replacing equations (6.1.3), (6.1.4), (6.1.6), and (6.1.7) in the weak form equation (5.7.4) we obtaine:
 * $$\displaystyle (6.1.7) $$
 * }
 * }


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

$$\left [ w \right ]^{eT}\int_{0}^{l}\left [ B \right ]^{T}\left ( 2+3x \right )\left [ B \right ]dx\left [ u \right ]^{e}=5\left [ w\right ]^{eT}\int_{0}^{l}\left [ B \right ]^{T}x\times dx-30w^{1}$$
 * $$\displaystyle (6.1.8) $$
 * }
 * }

The last equation is the weak form for any element; that is why the integration was done from 0 to l where l is the lenght of the element. Also in the last equation the factor $$w^{1}$$ in the nodal value of w at x = 1; therefore that factor is present in the weak form for the last element. Cancelling the factor $$\left [ w \right ]^{eT}$$ in equation (6.1.8)we obtain an equation of the form:$$Kd=f$$; therefore, we can extract from that equation the stiffness matrix which is:


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

$$k^{e}=\int_{0}^{l}\left [ B \right ]^{T}\left ( 2+3x \right )\left [ B \right ]dx$$
 * $$\displaystyle (6.1.9) $$
 * }
 * }

We developed a Matlab code (shown below) to calculate the global stiffness matrix $$\tilde{K}$$ using equation (1) p.30-5. The obtained matrix was:

Compute $$\tilde{k}^e$$ and $$L^e$$ for e = 1,2
Using the weak form developed above we calculate the stiffness matrix of the elements 1 and 2. As each elemet has 3 degrees of freedom and the total number of degrees of freedom is 7, the Location matrix is 3X5. The stiffness matrix are:


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

$$k^{1}=\begin{bmatrix} 10.83 & -12.66 &1.83 \\ -12.66 & 29.33 & -16.66\\ 1.83 & -16.66 & 14.83 \end{bmatrix} $$
 * $$\displaystyle (6.1.11) $$
 * }
 * }


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

$$k^{2}=\begin{bmatrix} 17.83 & -20.66 &2.83 \\ -20.66 & 45.33 & -24.66\\ 2.83 & -24.66 & 21.83 \end{bmatrix} $$ And the location matrix are:
 * $$\displaystyle (6.1.12) $$
 * }
 * }


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

$$L^{2}=\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & \\ 0 & 1 & 0 &0 & 0 & 0 & \\ 0 & 0 & 1 & 0 &0  & 0 & \end{bmatrix} $$
 * $$\displaystyle (6.1.13) $$
 * }
 * }


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

$$L^{2}=\begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 & \\ 0 & 0 & 0 &0 & 1 & 0 & \\ 0 & 0 & 0 & 0 &0  & 1 & \end{bmatrix} $$
 * $$\displaystyle (6.1.14) $$
 * }
 * }

Now the global matrix are calculated by using equation (2) p.31-1 which is repeated here:


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

$$\tilde{k}^{e}=L^{eT}k^{e}L^{e} $$
 * $$\displaystyle (6.1.15) $$
 * }
 * }

Therefore, for the element 1 whe have:


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

$$\tilde{k}^{1}=\begin{bmatrix} 1 & 0 &0 \\ 0 & 1 &0 \\ 0 & 0 &1 \\ 0 &0 &0 \\ 0 & 0 &0 \\ 0 & 0 & 0\\  0& 0 & 0 \end{bmatrix}\begin{bmatrix} 10.83 & -12.66 & 1.83\\ -12.66 & 29.33 &-16.66 \\ 1.83 & -16.66 & 14.83 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0& 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{bmatrix}=\begin{bmatrix} 10.83 & -12.66 & 1.83 & 0 & 0\\ -12.66 & 29.33 & -16.66 & 0 &0 \\ 1.83& -16.66 & 14.84 & 0 & 0\\  0& 0 & 0 & 0 & 0\\  0& 0 & 0 &0  & 0 \end{bmatrix} $$
 * $$\displaystyle (6.1.16) $$
 * }
 * }

In similar way for the elment 2 we have:


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

$$\tilde{k}^{2}=\begin{bmatrix} 0 & 0 &0 \\ 0 & 0 &0 \\ 0 & 0 &0 \\ 0 &0 &0 \\ 1 & 0 &0 \\ 0 & 1 & 0\\  0& 0 & 1 \end{bmatrix}\begin{bmatrix} 17.83 & -20.66 & 2.83\\ -20.66 & 45.33 &-24.66 \\ 2.83 & -24.66 & 21.83 \end{bmatrix}\begin{bmatrix} 0 & 0 & 0& 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 &0 \\ 0& 0& 17.83 & -20.66 & 2.83\\  0& 0 & -20.66 & 45.33 & -24.66\\  0& 0 &2.830 &-24.66  & 21.83 \end{bmatrix}$$
 * $$\displaystyle (6.1.17) $$
 * }
 * }

Performing the sum we obtain:

We can observ that this equation is equal to equation (6.1.10) because it corresponds to the same system with the same weak form. The location matrix is only a methodology for asembly the global matrix and therefore it can not change the final result.

Plot all the QLEBF for nel = 3
Qubic Lagrange Element Basis Function for every element are given by:


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

$$L_{1}=\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}; L_{2} =\frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)}; L_{3} = \frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}u_3 $$
 * $$\displaystyle (6.1.19) $$
 * }
 * }

Where $$x_{1}, x_{2}, and x_{3}$$, are the coordinates of the nodes 1, 2, and 3 respectively. This funtions are repeated on all the elements. Figure 6.1.1 shows the QLEBF for 3 elements, where it can be observed that every basis funtion has the value of 1 at the respective node and 0 on the others, of course the basis function are valid only inside their respective element. In the Figure, Li,j corresponds to j the basis funtion of the element i



Plot the exact solution, the numerical solution, and the error as a function of the number of elements
The analitical solution was developed before in other homeworks; therefore,here we only copy it

Figures 6.1.2 shows the analitical and numerical solution obtained with different number of elements. Although good aproximation was achieved even whit low number of elements, the convergence with a tolerance of $$10^{-6}$$ was not achieved even with more than 100 elements. With a tolerance of $$10^{-4}$$ convergence was achieved with 54 elements. Figure 6.1.3 shows the error as a function of the number of elements. In that figure, it can be observed that initially the error decreases fast but from aproximately 25 elements the relation between the number of elements and the error in logaritmic sacale tend to be linear and to achieve a tolerance of $$10^{-6}$$ more than 200 elements could be needed; but, is probalby that with higher number of elements the error decrease even slower because with more element the round-off error increases.





Appendix
 Matlab Code: