User:Eml5526.s11.team2.penultimate yay

=Problem 6.1: Solving G1DM1.0/D1b using Quad Element Lagrangian Basis Functions with Uniform Discretization =

Given: Strong Form and Data Set with Boundary Conditions
The Strong Form from lecture 9-1 and Fish and Belytschko p. 73:

1. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K}} $$ for $$ n_{el} = 2 $$
For $$ \displaystyle n_{el} = 2 $$, compute:

Where $$\displaystyle \mathbf{\tilde{ K^{e}}} $$ is defined by

3. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K^{e}}} $$ By Virtue of the Elemental Stiffness Matrix and Scattering Matrix for $$ n_{el} = 2 $$
For $$ \displaystyle e = 1, 2 $$ and compare these results to (1).

5. Plot the Approximate Solution and Analytic Solution and Error

 * For $$ \displaystyle n_{el} = 2,4,6,8,... $$ until convergence of the order $$ \displaystyle 10^{-6} $$ is achieved:
 * Plot $$ \displaystyle {u^h}_\tilde{n} $$ and $$ \displaystyle u $$ for the provided domain.
 * Plot the log of the error defined by, $$ \displaystyle E = log|{u^h}_\tilde{n}(0.5) - u(0.5)| $$ versus the number of degrees of freedom for the quadratic interpolation.

1. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K}} $$ for $$ n_{el} = 2 $$
Nodes 1, 2, 3, 4 and 5 values are $$ \displaystyle x_1 = 0 $$, $$ \displaystyle x_2 = \frac{1}{4} $$ , $$ \displaystyle x_3 = \frac{1}{4} $$ , $$ \displaystyle x_4 = \frac{3}{4} $$ and  $$ \displaystyle x_5 =1 $$  respectively as denoted in the  figure 6.1.1

To compute $$ \displaystyle \tilde{K^1} $$ and $$ \displaystyle \tilde{K^2} $$ it is necessary to establish the element stiffness matrix based on nodal positions and global elements. These take the following form and the indices values are defined by (6.1.7)

The basis function needs to be defined in terms of the Quadratic Lagrange Basis Function which has the following representation

Using $$ \displaystyle x_1 = 0 $$, $$ \displaystyle x_2 = \frac{1}{4} $$ , $$ \displaystyle x_3 = \frac{1}{4} $$ , $$ \displaystyle x_4 = \frac{3}{4} $$ and  $$ \displaystyle x_5 =1 $$ the first & second element basis function have the following form:

Taking derivatives of the basis functions provides the following equations

Having defined the derivatives of the basis functions the weak form can be used to create the stiffness matrix for each element using (6.1.7) which is defined as $$ \displaystyle {K^e}_{ij} = {\int\limits_ {{b^e}_i} ^\prime }({x^e}){a_2}({x^e}){b^e}{_j^\prime }({x^e})d{x^e} $$ which leads to the following values for (6.1.10) and (6.1.11)

Expanding each element stiffness matrix $$ \tilde{K}^1 $$ and $$ \tilde{K}^2 $$ to a 5x5 matrix so they can be summed as indicated in (6.1.6) provides the following Global stiffness matrix.

2. Determine the Elemental Stiffness Matrix $$ \mathbf{k^{e}} $$ and Scattering Matrix $$ \mathbf{L^{e}} $$ for $$ n_{el} = 2 $$
First, the scatter matrices can be computed for both elements 1 and 2:

Then the transpose of those matrices can be computed:

Making us of (6.1.7) and the basis functions derivatives $$ \displaystyle k^1 $$ and $$ \displaystyle k^2 $$ can be constructed for each index element

3. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K^{e}}} $$ By Virtue of the Elemental Stiffness Matrix and Scattering Matrix for $$ n_{el} = 2 $$
Performing matrix algebra on the scatter matrices scatter matrix transpose and $$ \displaystyle K^1 $$ and $$ \displaystyle  K^2 $$ provides the following

Comparing (6.1.26) and (6.1.25) to (6.1.18) it can be seen that the sum of the matrices is equivalent to the global stiffness matrix therefore the element stiffness matrices are equivalent for $$ \displaystyle e =1,2 $$.

Matlab Code
=Problem 6.2: Solving G1DM1.0/D1 using Quad Element Lagrangian Basis Functions with Uniform Discretization =

Given: Strong Form and Data Set with Boundary Conditions
The Strong Form from lecture 9-1 and Fish and Belytschko p. 73:

1. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K}} $$ for $$ n_{el} = 2 $$
For $$ \displaystyle n_{el} = 2 $$, compute:

Where $$\displaystyle \mathbf{\tilde{ K^{e}}} $$ is defined by

3. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K^{e}}} $$ By Virtue of the Elemental Stiffness Matrix and Scattering Matrix for $$ n_{el} = 2 $$
For $$ \displaystyle e = 1, 2 $$ and compare these results to (1).

5. Plot the Approximate Solution and Analytic Solution and Error

 * For $$ \displaystyle n_{el} = 2,4,6,8,... $$ until convergence of the order $$ \displaystyle 10^{-6} $$ is achieved:
 * Plot $$ \displaystyle {u^h}_\tilde{n} $$ and $$ \displaystyle u $$ for the provided domain.
 * Plot the log of the error defined by, $$ \displaystyle E = log|{u^h}_\tilde{n}(0.5) - u(0.5)| $$ versus the number of degrees of freedom for the quadratic interpolation.

1. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K}} $$ for $$ n_{el} = 2 $$
Nodes 1, 2, 3, 4 and 5 values are $$ \displaystyle x_1 = 0 $$, $$ \displaystyle x_2 = \frac{1}{4} $$ , $$ \displaystyle x_3 = \frac{1}{4} $$ , $$ \displaystyle x_4 = \frac{3}{4} $$ and  $$ \displaystyle x_5 =1 $$  respectively as denoted in the  figure 6.1.1

To compute $$ \displaystyle \tilde{K^1} $$ and $$ \displaystyle \tilde{K^2} $$ it is necessary to establish the element stiffness matrix based on nodal positions and global elements. These take the following form and the indices values are defined by (6.1.7)

The basis function needs to be defined in terms of the Quadratic Lagrange Basis Function which has the following representation

Using $$ \displaystyle x_1 = 0 $$, $$ \displaystyle x_2 = \frac{1}{4} $$ , $$ \displaystyle x_3 = \frac{1}{4} $$ , $$ \displaystyle x_4 = \frac{3}{4} $$ and  $$ \displaystyle x_5 =1 $$ the first & second element basis function have the following form:

Taking derivatives of the basis functions provides the following equations

Having defined the derivatives of the basis functions the weak form can be used to create the stiffness matrix for each element using (6.1.7) which is defined as $$ \displaystyle {K^e}_{ij} = {\int\limits_ {{b^e}_i} ^\prime }({x^e}){a_2}({x^e}){b^e}{_j^\prime }({x^e})d{x^e} $$ which leads to the following values for (6.1.10) and (6.1.11)

Expanding each element stiffness matrix $$ \tilde{K}^1 $$ and $$ \tilde{K}^2 $$ to a 5x5 matrix so they can be summed as indicated in (6.1.6) provides the following Global stiffness matrix.

2. Determine the Elemental Stiffness Matrix $$ \mathbf{k^{e}} $$ and Scattering Matrix $$ \mathbf{L^{e}} $$ for $$ n_{el} = 2 $$
First, the scatter matrices can be computed for both elements 1 and 2:

Then the transpose of those matrices can be computed:

Making us of (6.1.7) and the basis functions derivatives $$ \displaystyle k^1 $$ and $$ \displaystyle k^2 $$ can be constructed for each index element

3. Determine the Global Stiffness Matrix $$ \mathbf{\tilde{K^{e}}} $$ By Virtue of the Elemental Stiffness Matrix and Scattering Matrix for $$ n_{el} = 2 $$
Performing matrix algebra on the scatter matrices scatter matrix transpose and $$ \displaystyle K^1 $$ and $$ \displaystyle  K^2 $$ provides the following

Comparing (6.1.26) and (6.1.25) to (6.1.18) it can be seen that the sum of the matrices is equivalent to the global stiffness matrix therefore the element stiffness matrices are equivalent for $$ \displaystyle e =1,2 $$.

Matlab Code
=Problem 6.3: Using Calculix to Reproduce Beede's Models=

Calculix is an open source program that we can solve our FEA problems. Bconverged is a really good provider of calculix that would offer engineers work between different platforms. After downloading calculix software from bconverged we will have two modules of calculix and text editor called SciTE. These two modules that belongs to calculix are $$ \displaystyle\mathit{cgx} $$ and $$\displaystyle\mathit{ccx} $$. In the first module of calculix we can design our solid models and define boundary conditions on it. The solver,crunchix (ccx), is solver and post processor.Lets start building models that are presented by Beede's.Beede's examples. Also this cite and more information is provided by Dr. Vu-Quoc in our EML5526 FEA class page. Class page

FEA BEAM PART 1: Getting Started
As a first step we have to open SciTE text editor that is a place where we type our commands and view our models. To run it


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

$$  \displaystyle [Start] \Rightarrow [All Programs]\Rightarrow [bConverged]\Rightarrow [SciTE] $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Here are some information about SciTE editor.



This is very important to understand that we will use two typing methods to type our commands. This i a very beautiful property of calculix. We can type in text editor and that 'run ' all of the codes in it or we can type directly by clicking on the graph screen and execute the code by hitting enter.


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

$$  \displaystyle [open SciTE]\Rightarrow [File]\Rightarrow [Save As]\Rightarrow [beam.fbd] $$
 * style="width:95%" |
 * style="width:95%" |
 * }

The step presented above is a essential step. Because you can not run the codes without saving!!!


 * Type these codes into your SciTE editor. Then hit F10 key or click tools--post processor to run.
 * pnt p1 0 0 0
 * pnt p2 100 0 0
 * line l1 p1 p2 25



Let's get familiar with our codes.


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

$$  \displaystyle \begin{matrix} pnt \Rightarrow name & x coord. & y coord. & z coord. & \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }


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

$$  \displaystyle \begin{matrix} line \Rightarrow name & start point & end point & division & \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }


 * Proceed with expanding commands in text editor and creating four lines and five nodes. We are partially creating this big line because if we had constructed the line in one step the mesh would have been arbitrary. So to avoid arbitrariness we are doing more steps. We want equal divisions among mesh elements.


 * Here how we can view our model in different aspects. In the area surrounded by thin black surface use this properties of mouse.


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

$$  \displaystyle \begin{matrix} LMB & rotate \\ MMB & zoom\\ RMB & pan \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }


 * While cgx window is active it means that we are working in 'interactive mode.' Type these commands while cgx window is active.
 * plot pa all
 * plus la all


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

$$  \displaystyle \begin{matrix} plot & n(nodes) & e(elements) &f(faces) & p(points) & l(lines) & s(surfaces) & b(bodies) & all(set name) \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }


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

$$  \displaystyle \begin{matrix} plus & n(nodes) & e(elements) &f(faces) & p(points) & l(lines) & s(surfaces) & b(bodies) & all(set name) \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Commands are same for plot and plus. But there is a slight difference between them. Plot command plots only one type of entity removes other entities from graph screen but plus command adds entities to the other entities that are already visual in graph screen.


 * Type seta lines l l1 l2 l3 l4
 * swep lines sweplines tra 0 10 0 10
 * seta surfaces s A001 A002 A003 A004
 * swep surfaces swepsurfaces tra 0 0 1 1


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

$$  \displaystyle \begin{matrix} seta & type & elements. in.set & \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }


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

$$  \displaystyle \begin{matrix} swep & set & fx &fy & fz &  tra& dx & dy & dz \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }



As we see above we have created body from swepping surfaces, surfaces from swepping lines. When we create a high order entity lower order entities will be kept in order to define the high order one. For example if we create surface lines and points will be kept.

After constructing the body now we can mesh our body.We have different kinds of elements for 1D,2D,and 3D meshes. Type these commands in interactive mode.
 * elty all he8
 * mesh all
 * plot m all


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

$$  \displaystyle \begin{matrix} elty & set & meshtype \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }


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

$$  \displaystyle \begin{matrix} mesh & set \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }



By using different types of views we can get different aspects of our mesh and other entities.

FEA BEAM PART 2: Exporting Meshes, Loads and Boundary Conditions
In this section again we start with our text editor and open beam.fbd. Now we are updating written codes in it. New codes
 * seto beam
 * pnt commands
 * line commands
 * sets and swep codes
 * setc beam
 * defining and setting up mesh

Between seto and setc commands we have defined set named beam and we put all of the entities in this general set. Also we have defined other sets just before we did in the last section.If we hit the F10 key and run our new codes beam.msh will be created automatically.Opening beam.msh yields following



As we see above two sets defined in this document. Node and element sets. Looking at element sets we can see that elements are 3D elements with 8 nodes. The style of presenting data and storing is calculix convenience.



Here we have defined another set named fixed to assign fixed nodes in it. View the appropriate section of the beam and type quadd fixed while in interactive mode. Then press r and move the cursor down and press r to create select box. Be careful to create select box that can fit all tha nodes that we are interested in. Press a to active mode a then press n to select nodes in the box. You will see the node properties in the output section and by pressing q we will have finished our job.

As mentioned in the article we have two different boundary conditions. We have defined fixed dots where we will apply restriction on their displacement. Also turn the view to see faces and viewing elements,we found element '751'.And we have applied pressure force its upward surface. And we applied two force vectors on its two nodes.

We can show element correspondence with their families.

FEA BEAM PART 3
In this section as you see in the figure we created new text file but with an extension .inp. This means that now we are using ccx commands. Generally start with '*'. First we are taking input parameter from *.msh and *.123 documents. Between STEP and ENDSTEP we are also defining material properties and selensting Solve and Post process tabs from the tools bar we can find stresses along the beam.



FEA BEAM PART 4
In this section we are interested in how we can visualize the results and interpret them easily. We can not see how much the bar deformed in the first glance but after changing scale factor we can be more sure about deformations.
 * {| style="width:100%" border="0"

$$  \displaystyle Deformed Shape $$
 * style="width:95%" |
 * style="width:95%" |
 * }




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

$$  \displaystyle \begin{matrix} Part & Of  & Animation \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Another way to visualize the results is using animation window.Just go to animate bar and select start to. The deformation will be between 100 percent in both directions but one can change its limits.


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

$$  \displaystyle \begin{matrix} Y stress & At  & Portion \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

=Problem 6.4: Verifying the divergence theorem=

1. The Following Vector Field (FB 6.1)
Defined by the components,
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \mathbf{q_x} &= - {y^2}\mathbf{i} \\ \mathbf{q_y} &= -2xy\mathbf{j} \end{align} $$     (6.4.1)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Vector Field and Domain

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

Figure 1: Vector Field and Domain.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Mathematica Code:
Mathematica code used to generate Figure 2: Grid[Partition[ VectorPlot[{3*x^2*y + y^3, 3*x + y^3}, {x, -2, 2}, {y, 0, 4},VectorColorFunction -> Hue, VectorScale -> #] & /@ {Medium, Large}, 2]]|undefined

2. The Following Vector Field (FB 6.2)

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

$$ \displaystyle \begin{align} \mathbf{q_x} &= \left( 3{x^2}y + {y^3} \right) \mathbf{i} \\ \mathbf{q_y} &= \left( 3x+{y^3} \right) \mathbf{j} \end{align} $$     (6.4.2)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Vector Field and Domain

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

Figure 2: Vector Field and Domain.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Mathematica Code:
Mathematica code used to generate Figure 1: VectorPlot[{-y^2, -2*x*y}, {x, -1, 1}, {y, -1, 1}, VectorColorFunction -> Hue, VectorScale -> Medium]

3. The Divergence Theorem Proof (FB 6.3)
Which relates the area integral of the divergence of a vector field to the contour integral of a vector field and is mathematically expressed as,
 * {| style="width:100%" border="0"

$$ \displaystyle \displaystyle \int_\gamma \mathbf{\nabla q} d\gamma = \oint_\Gamma {\mathbf{n}d\Gamma} $$     (6.4.3)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

3. Using the divergence theorem
Prove:

$$ \displaystyle \oint_\Gamma  {\mathbf{n}d\Gamma} = 0 $$

Solution
The divergence theorem states:
 * {| style="width:100%" border="0"

$$ \displaystyle \int\limits_\Omega {{\nabla ^T}{\mathbf{q}}d\Omega }  = \oint_\Gamma  {{{\mathbf{q}}^T}{\mathbf{n}}d\Gamma } $$     (6.4.4)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

1. The Vector Field (FB 6.1)
Defined by the components,
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \mathbf{q_x} &= - {y^2}\mathbf{i} \\ \mathbf{q_y} &= -2xy\mathbf{j} \end{align} $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

For the given vector field the left hand side of the divergence theorem is as follows


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

$$ \displaystyle \begin{align} \int\limits_\Omega {\nabla {\mathbf{q}}d\Omega }  &= \int\limits_{ - 1}^1 {\int\limits_{ - 1}^1 {\left( {\frac + \frac} \right)dydx} } \\ &=\int\limits_{ - 1}^1 {\int\limits_{ - 1}^1 {\left( { - 2x} \right)dxdy = \int_{ - 1}^1 {\left. { - {x^2}} \right|_{ - 1}^1dy} } } \\ &=\int_{ - 1}^1 { 0dy} \\ &=0 \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.5)
 * 
 * }

For the given vector field the right hand side of the divergence theorem is


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

$$ \displaystyle \oint\limits_\Gamma {{{\mathbf{q}}^T}{\mathbf{n}}d\Gamma }  = \int\limits_{AB} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(1)}}d\Gamma }  + \int\limits_{BC} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma }  + \int\limits_{CD} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma }  + \int\limits_{DA} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(4)}}d\Gamma } $$     (6.4.6)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Where $${\mathbf{n}^{(1)}},{\mathbf{n}^{(2)}},{\mathbf{n}^{(3)}},{\mathbf{n}^{(4)}}$$ can be found in Fig. 1 and


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

$$ \displaystyle \begin{align} \int\limits_{AB} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(1)}}d\Gamma } &= \int_{ - 1}^1 {\left( { - {y^2}{\mathbf{i}} - 2xy{\mathbf{j}}} \right) \cdot \left( { - {\mathbf{j}}} \right)dx} \\ &=\int_{ - 1}^1 {2x\underbrace y_{y = - 1}dx} =\int_{ - 1}^1 {-2xdx} \\ &=0 \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.7)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \int\limits_{BC} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma } &= \int_{ - 1}^1 {\left( { - {y^2}{\mathbf{i}} - 2xy{\mathbf{j}}} \right) \cdot \left( {\mathbf{i}} \right)dy}  \\ &=\int_{ - 1}^1 { - {y^2}dy} \\ &=-\frac{2}{3} \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.8)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \int\limits_{AB} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(3)}}d\Gamma } &= \int_{ - 1}^1 {\left( { - {y^2}{\mathbf{i}} - 2xy{\mathbf{j}}} \right) \cdot \left( {  {\mathbf{j}}} \right)dx} \\ &=\int_{ - 1}^1 {-2x\underbrace y_{y =  1}dx} =\int_{ - 1}^1 {-2xdx} \\ &=0 \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.9)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \int\limits_{BC} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(4)}}d\Gamma } &= \int_{ - 1}^1 {\left( { - {y^2}{\mathbf{i}} - 2xy{\mathbf{j}}} \right) \cdot \left( -{\mathbf{i}} \right)dy}  \\ &=\int_{ - 1}^1 { {y^2}dy}  \\ &=\frac{2}{3} \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.10)
 * 
 * }

Continuing from equation 4.6.3


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

$$ \displaystyle \begin{align} \oint\limits_\Gamma {{{\mathbf{q}}^T}{\mathbf{n}}d\Gamma }  &= \int\limits_{AB} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(1)}}d\Gamma }  + \int\limits_{BC} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma }  + \int\limits_{CD} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma }  + \int\limits_{DA} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(4)}}d\Gamma } \\ &=0+\left( { - \frac{2}{3}} \right)+0+\frac{2}{3}\\ &=0 \end{align} $$     (6.4.10)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

From the equivalence of equations 4.6.2 and 4.6.8 we see that the divergence theorem holds for $${{\mathbf{q}}_{(x,y)}} = \left( { - {y^2}} \right){\mathbf{i}} + \left( { - 2xy} \right){\mathbf{j}}$$

2. The Vector Field (FB 6.2)

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

$$ \displaystyle \begin{align} \mathbf{q_x} &= \left( 3{x^2}y + {y^3} \right) \mathbf{i} \\ \mathbf{q_y} &= \left( 3x+{y^3} \right) \mathbf{j} \end{align} $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

For the given vector field the left hand side of the divergence theorem is as follows


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

$$ \displaystyle \begin{align} \int\limits_\Omega {\nabla {\mathbf{q}}d\Omega }  &= \int_{ - 2}^2 {\int_0^{ - {x^2} + 4} {\left( {6xy + 3{y^2}} \right)dydx} }  \\ &=\int_{ - 2}^2 {\left. {\left( {3x{y^2} + {y^3}} \right)} \right|} _0^{y = - {x^2} + 4}dx \\ &=\int_{ - 2}^2 {\left[ {3x\left( { - {x^2} + 4} \right) + {{\left( { - {x^2} + 4} \right)}^3}} \right]} dx \\ &=\int_{ - 2}^2 {\left( { - {x^6} + 3{x^5} + 12{x^4} - 24{x^3} - 48{x^2} + 48x + 64} \right)} dx\\ &=\frac{4096}{35} \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.11)
 * 
 * }

For the given vector field the right hand side of the divergence theorem is


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

$$ \displaystyle \oint\limits_\Gamma {{{\mathbf{q}}^T}{\mathbf{n}}d\Gamma }  = \int\limits_{AB} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(1)}}d\Gamma }  + \int\limits_{BA} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma } $$     (6.4.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Where $${\mathbf{n}^{(1)}}$$ and $${\mathbf{n}^{(2)}}$$ can be found in Fig 2 above and


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

$$ \displaystyle \begin{align} \int\limits_{AB} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(1)}}d\Gamma } &= \int_{ - 2}^2 {\left[ {\left( {3{x^2}y + {y^3}} \right){\mathbf{i}} + \left( {3x + {y^3}} \right){\mathbf{j}}} \right]}  \cdot \left( { - {\mathbf{j}}} \right)dx \\ &=\int_{ - 2}^2 {\underbrace {\left( { - 3x - {y^3}} \right)}_{y = 0\forall x \in \left[ { - 2,2} \right]}} dx \\ &= \int_{ - 2}^2 { - 3xdx}\\ &=0 \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.13)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \int\limits_{BA} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma } &=\int_{ - 2}^2 {\left[ {\left( {3{x^2}y + {y^3}} \right){\mathbf{i}} + \left( {3x + {y^3}} \right){\mathbf{j}}} \right]} \cdot \left( 2x{\mathbf{i}} + {\mathbf{j}} \right)\left( { - dx} \right)\\ &= - \int_{ - 2}^2 {\underbrace {\left[ {2x\left( {3{x^2}y + {y^3}} \right) + \left( {3x + {y^3}} \right)} \right]}_{y = - {x^2} + 4}} dx\\ &= - \int_{ - 2}^2 {\left[ {2x\left( {3{x^2}\left( { - {x^2} + 4} \right) + {{\left( { - {x^2} + 4} \right)}^3}} \right) + \left( {3x + {{\left( { - {x^2} + 4} \right)}^3}} \right)} \right]} dx\\ &=-\int_{ - 2}^2 {\left( { - 2{x^7} - {x^6} + 18{x^5} + 12{x^4} - 72{x^3} - 48{x^2} + 131x + 64} \right)} dx\\ &=\frac \end{align}
 * style="width:95%" |
 * style="width:95%" |

$$     (6.4.14)
 * <p style="text-align:right">
 * }

Continuing from equation 4.6.10


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

$$ \displaystyle \begin{align} \oint\limits_\Gamma {{{\mathbf{q}}^T}{\mathbf{n}}d\Gamma }  &= \int\limits_{AB} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(1)}}d\Gamma }  + \int\limits_{BA} {{{\mathbf{q}}^T}{{\mathbf{n}}^{(2)}}d\Gamma }\\ &= 0+\frac{4096}{35} \end{align} $$     (6.4.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

From the equivalence of equations 4.6.9 and 4.6.13 we have shown that the divergence theorem holds for $${{\mathbf{q}}_{(x,y)}} = \left( {3{x^2}y + {y^3}} \right){\mathbf{i}} + \left( {3x + {y^3}} \right){\mathbf{j}}$$

3. The Divergence Theorem Proof (FB 6.3)
The given equation to prove is


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

$$ \displaystyle \oint\limits_\Gamma {{\mathbf{n}}d\Gamma }  = 0$$ (6.4.16) From equation 6.4.1, divergence theorem, we see that for 6.4.14 to be true the following must be true
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \int\limits_\Omega {{\nabla ^T}{\mathbf{q}}d\Omega }  = 0$$ (6.4.17) From 6.4.14 we know that $$\mathbf{q}$$ must equal 1. Therefore equation 6.4.15 becomes
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \int\limits_\Omega {{\nabla ^T}{\mathbf{q}}d\Omega }  &= \int\limits_\Omega  {{\nabla ^T}(1)d\Omega }\\ &=\int\limits_\Omega {\left( {\frac + \frac +  \ldots } \right)d\Omega } \\ &= \int\limits_\Omega {0d\Omega } \\ &=0 \end{align} $$     (6.4.18) Then from the divergence theorem equation 6.4.14 must be true
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

=Problem 6.5: Two Dimensional Lagrange Interpolation Basis Function Problem=

Given: The Following 2D Data Set
Given the data set provided below.

Domain: Bi-Unit Square
1) The following domain information $$ \Omega = \bar \omega  = \square $$ which is for a bi-unit square

Equation: PDE, K, f, and Initial Condition
2) The partial differential equation of the form $$ \displaystyle \operatorname{div} \left ( K  \cdot\nabla{u} \right ) + f = \rho c\frac $$ 2.a) $$  {K}= {I}=I_n = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix} $$

2.b) $$ f=0, \frac = 0 $$

Find: 2-D Lagrange Interpolation Basis Functions (LIBF) and Solve the PDE for $$ m = n = 2,4,6,8 $$
Use 2-dimensional Lagrange Interpolation Basis Functions (LIBF) with $$ \displaystyle m = n = 2,4,6,8 $$ nodes per element to solve the data set above. Ensure the accuracy of the approximate solution $$ \displaystyle u^h(x = 0, y = 0) = 10^{-6} $$ at the center of the element $$ \displaystyle (x,y) = (0,0) $$

The 2-D LIBF can be defined as in (6.5.1)

Where $$ \displaystyle L_{i,m}(x) $$ and  $$ \displaystyle L_{j,n}(y) $$ represent Lagrange Interpolation Basis Functions

Solution
Per class lecture notes Lecture (3) 35-4 the Weak Form for the 2-D element can be defined as follows, given the steady-state nature and lack of body forces & natural boundary condition on the element and kappa being defined as an identity matrix

m = n = 2
The element for m = n = 2 consist of 4 nodes. The representation of the element is indicated below in the figure.

The element is a 2-D element and therefore the weak form consist of a double integral with limits of integration from $$ \displaystyle x = -1 $$ to $$ \displaystyle x = 1 $$ and $$ \displaystyle y = -1 $$ to $$ \displaystyle y = 1 $$. This leads to (6.5.3) being updated in the following matter.

The following step is to modify (6.5.4) to use the proper basis functions. For this 2-D element the 2-D LIBF has been defined in (6.5.1). Modifying the weak form to represent the proper basis functions results in the following.

Now converting the gradient operator to vector form and expanding the identity matrix to a 2 by 2 matrix.

Now it is necessary to define the LIBFs for the case where $$ \displaystyle  m = n = 2 $$ this provides the following values for $$ \displaystyle  N_I  $$ where $$ \displaystyle I = 1 ,2 ,3 ,4 $$

Solving for the Element Stiffness Matrix using (6.5.6) and the relation between Element Nodes and indices of the Element Stiffness Matrix where $$ \displaystyle I = i + (j-1)*m $$ leads to the following

Now establishing the system of equations for the element provides the following.

Solving provides the following values.

Based on the uniform essential boundary condition of $$ \displaystyle g = 2 $$ all displacements are 2 for all m = n.

Matlab Code

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

$$  \displaystyle \begin{matrix} d & matrix \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }




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

$$  \displaystyle \begin{matrix} Error \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }



=Problem 6.6: Solve a Given Data Set Using 2D Linear Lagrangian Element Basis Function=

Given: The Following 2D Data Set, G2DM1.0/D1
Given the data set provided below.

Domain: Bi-Unit Square
1) The following domain information $$ \Omega = \bar \omega  = \square $$ which is for a bi-unit square

Equation: PDE, K, f, and Initial Condition
2) The partial differential equation of the form $$ \displaystyle \operatorname{div} \left ( K  \cdot\nabla{u} \right ) + f = \rho c\frac $$ 2.a) $$  {K}= {I}=I_n = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix} $$

2.b) $$ f=0, \frac = 0 $$

Find: Stiffness Matrix, Approximate Solution $$ u_h $$, and Compare to u for GLQ
1) Solve G2DM1.0/D1 with 2D LLEBF until an accuracy of $$ 10^{-6} $$ 2)  Integrate  $$  \underline{k^{e}} $$ in parent coordinate and find u exactly with GLQ (Gauss-Legendre Quadrature) 3) Use distorted quad meshes, compare results with uniform meshes (CGX)

LLEBF
Per class lecture notes Lecture (3) 35-4 the Weak Form for the 2-D element can be defined as follows, given the steady-state nature and lack of body forces & natural boundary condition on the element and kappa being defined as an identity matrix

m = n = 2
The element for m = n = 2 consist of 4 nodes. The representation of the element is indicated below in the figure.

The element is a 2-D element and therefore the weak form consist of a double integral with limits of integration from $$ \displaystyle x = -1 $$ to $$ \displaystyle x = 1 $$ and $$ \displaystyle y = -1 $$ to $$ \displaystyle y = 1 $$. This leads to (6.5.3) being updated in the following matter.

The following step is to modify (6.5.4) to use the proper basis functions. For this 2-D element the 2-D LIBF has been defined in (6.6.1). Modifying the weak form to represent the proper basis functions results in the following.

Now converting the gradient operator to vector form and expanding the identity matrix to a 2 by 2 matrix.

Now it is necessary to define the LIBFs for the case where $$ \displaystyle  m = n = 2 $$ this provides the following values for $$ \displaystyle  N_I  $$ where $$ \displaystyle I = 1 ,2 ,3 ,4 $$

Solving for the Element Stiffness Matrix using (6.5.6) and the relation between Element Nodes and indices of the Element Stiffness Matrix where $$ \displaystyle I = i + (j-1)*m $$ leads to the following

Now establishing the system of equations for the element provides the following.

Solving provides the following values.

Based on the uniform essential boundary condition of $$ \displaystyle g = 2 $$ all displacements are 2 for all m = n.

m = n = 3
The element for m = n = 2 consist of 4 nodes. The representation of the element is indicated below in the figure.

The element is a 2-D element and therefore the weak form consist of a double integral with limits of integration from $$ \displaystyle x = -1 $$ to $$ \displaystyle x = 1 $$ and $$ \displaystyle y = -1 $$ to $$ \displaystyle y = 1 $$. This leads to (6.6.3) being updated in the following matter.

The following step is to modify (6.6.4) to use the proper basis functions. For this 2-D element the 2-D LLEBF has been defined in (6.6.1). Modifying the weak form to represent the proper basis functions results in the following.

Now converting the gradient operator to vector form and expanding the identity matrix to a 3 by 3 matrix.

Now it is necessary to define the LLEBFs for the case where $$ \displaystyle  m = n = 3 $$ this provides the following values for $$ \displaystyle  N_I  $$ where $$ \displaystyle I = 1 ,2 ,4 ,5 $$ for Element 1.

Solving for the Element Stiffness Matrix using (6.5.6) and the relation between Element Nodes and indices of the Element Stiffness Matrix where $$ \displaystyle I = i + (j-1)*m $$ leads to the following

6.6 MATLAB Code:


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

$$  \displaystyle \begin{matrix} Error \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }



=Problem 6.7: Proof of Matrix algebra and Wolframalpha syntax=

Given: Two Matrices A and B

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

$$ \displaystyle A = \left[ {\begin{array}{*{20}{c}} 1&1&1\\ 2&{ - 1}&3\\ 3&2&6 \end{array}} \right] $$     (6.7.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle B = \left[ {\begin{array}{*{20}{c}} 1&3&5\\ 1&{ - 4}&1\\ 2&5&8 \end{array}} \right] $$     (6.7.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Find: Prove the validity of matrix operations
Verify that 1)$${\left( \right)^T} = {{\mathbf{B}}^T}{{\mathbf{A}}^T}$$ 2)$${\left(  \right)^T} = {\left(  \right)^{ - 1}} = {{\mathbf{A}}^{ - T}}$$

Also determine why Wolfram Alpha uses the syntax described in [[media:Fe1.s11.mtg36.djvu|Mtg 36-1,2]]

Solution to 6.7.1

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

$$ \displaystyle\begin{align} \mathbf{AB}&=\left[ {\begin{array}{*{20}{c}} 1&1&1\\ 2&{ - 1}&3\\ 3&2&6 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&3&5\\ 1&{ - 4}&1\\ 2&5&8 \end{array}} \right]\\ &=\left[ {\begin{array}{*{20}{c}} 4&4&{14}\\ 7&{25}&{33}\\ {17}&{31}&{65} \end{array}} \right]\end{align} $$     (6.7.3) Then
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle {\left( \right)^T}=\left[ {\begin{array}{*{20}{c}} 4&7&{17}\\ 4&{25}&{31}\\ {14}&{33}&{65} \end{array}} \right] $$     (6.7.4) Next we must show what the right hand side equals
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle\begin{align} {{\mathbf{A}}^T}&=\left[ {\begin{array}{*{20}{c}} 1&2&3\\ 1&{ - 1}&2\\ 1&3&6 \end{array}} \right]\\ {{\mathbf{B}}^T}&=\left[ {\begin{array}{*{20}{c}} 1&1&2\\ 3&{ - 4}&5\\ 5&1&8 \end{array}} \right] \end{align} $$     (6.7.5) So then
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle\begin{align} {{\mathbf{B}}^T}{{\mathbf{A}}^T}&=\left[ {\begin{array}{*{20}{c}} 1&2&3\\ 1&{ - 1}&2\\ 1&3&6 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&1&2\\ 3&{ - 4}&5\\ 5&1&8 \end{array}} \right]\\ &=\left[ {\begin{array}{*{20}{c}} 4&7&{17}\\ 4&{25}&{31}\\ {14}&{33}&{65} \end{array}} \right]\end{align} $$     (6.7.6)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Comparing equations 6.7.4 and 6.7.6 we can see that the two are equivalent, proving our statement.

Solution to 6.7.2
Starting with $${\left( \right)^T}$$, we have


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

$$ \displaystyle {{\mathbf{A}}^{ - 1}}=\left[ {\begin{array}{*{20}{c}} {1.5}&{.5}&{ - .5}\\ {.375}&{ - .375}&{.125}\\ { - .875}&{ - .125}&{.375} \end{array}} \right] $$     (6.7.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle {\left( \right)^T}=\left[ {\begin{array}{*{20}{c}} {1.5}&{.375}&{ - .875}\\ {.5}&{ - .375}&{ - .125}\\ { - .5}&{.125}&{.375} \end{array}} \right] $$     (6.7.8)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we will examine $${\left( \right)^{ - 1}}$$. From eq 6.7.5 we have $$\mathbf{A}^T$$ taking the inverse of this we get
 * {| style="width:100%" border="0"

$$ \displaystyle {\left( \right)^{ - 1}}=\left[ {\begin{array}{*{20}{c}} {1.5}&{.375}&{ - .875}\\ {.5}&{ - .375}&{ - .125}\\ { - .5}&{.125}&{.375} \end{array}} \right] $$     (6.7.9) Comparing equations 6.7.8 and 6.7.9 we see that the equivalence is upheld. Next we must verify the last portion of the equivalence, $${{\mathbf{A}}^{ - T}}$$.Using Matlab we see that
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle {{\mathbf{A}}^{ - T}}=\left[ {\begin{array}{*{20}{c}} {1.5}&{.375}&{ - .875}\\ {.5}&{ - .375}&{ - .125}\\ { - .5}&{.125}&{.375} \end{array}} \right] $$     (6.7.10) Comparing equations 6.7.8 through 6.7.10 we see that they are all equivalent. Proving our statement.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Solution to 6.7.3:Syntax
Wolfram Alpha, [WA], shows that the use of parentheses is for grouping; the use of braces or curly brackets, {}, is for lists; and brackets are used to denote the arguments of a function. Matrices are inputted into Wolfram Alpha as nested lists, hence the use of braces. So a simple 2x2 matrix would be written as
 * {| style="width:100%" border="0"

$$ \displaystyle \{ \{ 1,2\} ,\{ 3,4\} \} $$     (6.7.11) Multiplying two 2x2 matrices would be written as
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \{ \{ 1,2\} ,\{ 3,4\} \} *\{ \{ 1,2\} ,\{ 3,4\} \} $$      (6.7.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

From above, to take the transpose of the matrix you would expect the syntax to be


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

$$ \displaystyle transpose\left[ {\{ \{ 1,2\} ,\{ 3,4\} \} } \right] $$     (6.7.13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle transpose {\{ \{ 1,2\} ,\{ 3,4\} \} } $$     (6.7.14) While this will work, you can also leave off the brackets, as in eq 6.7.14, and the syntax will work. The issue comes when you use nested functions. If you wanted to take the inverse of the transpose or the transpose of the inverse you would have to write
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle invert\left[ {transpose\{ \{ 1,2\} ,\{ 3,4\} \} } \right] $$     (6.7.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle invert\left\{ {transpose\{ \{ 1,2\} ,\{ 3,4\} \} } \right\} $$     (6.7.16) Because there are nested functions the brackets are needed. Simply using braces will not work. Braces would actually be redundant. We already have a specified list from the original braces. Transposing the list is still a list. So by trying to use braces again as in eq 6.7.16, would be invalid syntax. Therefore the brackets are needed for this function as opposed to eq 6.7.13
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

=Problem 6.8: Gaussian Quadrature=

Given: The Legendre Polynomials
The Legendre Polynomians $$ P_{n}(x) $$ defined here
 * {| style="width:100%" border="0"

$$ \displaystyle {P}_{n}(x) = \sum_{i=0}^{\lfloor n/2 \rfloor} (-1)^{i} \cdot \frac {(2n-2i)! \cdot x^{n-2i}} {2^{n}\cdot i! \cdot (n-i)! \cdot (n-2i)!} $$     (8.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

1. Verify The Gaussian Quadrature
The points $$ \displaystyle x_i $$ and weights $$\displaystyle w_i $$ for $$\displaystyle n=1,2,3,4,5 $$.

2. Solve the Following:
The following problems can be found in Fish and Belytschko pg. 91-92:

Problem 4.6: Gauss Quadrature to Obtain Numerical Integration
Use the Gauss quadrature to obtain exact values for the following integrals. Also verify by analytic integration:

(c) Write MATLAB Code
That utilizes the function gauss.m and performs Gauss integration. Check the solutions in (a) and (b) against the MATLAB code.

Problem 4.7: Three-Point Gauss Quadrature to Obtain Numerical Integration
Use the three-point Gauss quadrature to evaluate the following integrals. Compare the results to the analytic integral.

(c) Write MATLAB Code
That utilizes the function gauss.m and performs Gauss integration. Check the solutions in (a) and (b) against the MATLAB code

Problem 4.8: n-Point Gauss Quadrature Comparison for $$ n=1,2,3 $$
The integral $$ \displaystyle \int_{-1}^{1}{\left(3\xi^3+2\right)d\xi} $$ can be integrated exactly using the two-point Gauss quadrature. How is the accuracy effected if

(c) Use MATLAB Code
Check the solutions in (a) and (b) against the MATLAB code

1. Verify The Gaussian Quadrature
As was discussed in class, $$ \displaystyle x_i $$ is the root of the $$ n-th $$ order Legendre Polynomial $$ \displaystyle P_{n}(x) $$ can be easily obtained for $$ n=1,2,3,4,5 $$

For n = 1
The root of the Legendre Polynomial $$ P_1(x)=0 $$ is found by solving

$$ \displaystyle P_1(x)=x=0 $$

And the weight, $$ w_{1} $$ is

$$ \displaystyle \begin{align} w_{1} & = \frac{-2}{(1+1)[P_{1}^{'}({{x}_{1}}){{P}_{2}}({{x}_{1}})]}  \\ & = \frac{-2}{(1+1)\left[ 1\cdot \frac{3{x_1}^{2}-1}{2}\right] }  \\ & = \frac{-2}{(1+1)\left[ 1\cdot \frac{3\cdot{0}^{2}-1}{2}\right] } \\ & = 2         \end{align} $$

For n = 2
The root of the Legendre Polynomial $$ P_2(x)=0 $$ are found by solving

$$ \displaystyle P_2(x)=(3x^2-1)=0 $$

Where the weights, $$ w_{1} $$ and $$ w_{2} $$ are

$$ \displaystyle \begin{align} w_{1} & = \frac{-2}{(2+1)[P_2^{'}(x_1)P_3(x_1)]}  \\ & = \frac{-2}{(2+1)\left[ 3x_1 \cdot \frac{5{x_1}^{3}-3x_1}{2}\right] }  \\ & = \frac{-2}{(2+1)\left[ 3\left( {\frac{1}{\sqrt 3 }}\right)^3 \cdot \frac{5\left( {\frac{1}{\sqrt 3 }}\right)^3 -3\left(\frac{1}{\sqrt 3 }\right) }{2}\right] } \\ & = 1         \end{align} $$

$$ \displaystyle \begin{align} w_{2} & = \frac{-2}{(2+1)[P_2^{'}(x_2)P_3(x_2)]}  \\ & = \frac{-2}{(2+1)\left[ 3x_2 \cdot \frac{5{x_2}^{3}-3x_2}{2}\right] }  \\ & = \frac{-2}{(2+1)\left[ 3\left( {-\frac{1}{\sqrt 3 }}\right)^3 \cdot \frac{5\left( {-\frac{1}{\sqrt 3 }}\right)^3 -3\left( -\frac{1}{\sqrt 3 }\right) }{2}\right] } \\ & = 1         \end{align} $$

For n = 3
The root of the Legendre Polynomial $$ P_3(x)=0 $$ are found by solving

$$ \displaystyle P_3(x)=(5x^3-3x)=0 $$

Where the weights, $$ w_{1}, w_{2}, $$ and $$ w_{3} $$ are

$${{W}_{1}}=\frac{-2}{(3+1)[P_{3}^{'}({{x}_{1}}){{P}_{4}}({{x}_{1}})]}=\frac{8}{9} $$

$${{W}_{2}}=\frac{-2}{(3+1)[P_{3}^{'}({{x}_{2}}){{P}_{4}}({{x}_{2}})]}=\frac{5}{9} $$

$${{W}_{3}}=\frac{-2}{(3+1)[P_{3}^{'}({{x}_{3}}){{P}_{4}}({{x}_{3}})]}=\frac{5}{9}$$

For n = 4
The root of the Legendre Polynomial $$ P_4(x)=0 $$ are found by solving

$$ \displaystyle P_4(x)=(35x^4-30x^2+3)=0 $$

Where the weights, $$ w_{1}, w_{2}, w_{3}$$ and $$ w_{4} $$ are

$${{W}_{1}}=\frac{-2}{(4+1)[P_{4}^{'}({{x}_{1}}){{P}_{5}}({{x}_{1}})]}=(18-\sqrt{30})/36 $$

$${{W}_{2}}=\frac{-2}{(4+1)[P_{4}^{'}({{x}_{2}}){{P}_{5}}({{x}_{2}})]}=(18-\sqrt{30})/36 $$

$${{W}_{3}}=\frac{-2}{(4+1)[P_{4}^{'}({{x}_{3}}){{P}_{5}}({{x}_{3}})]}=(18+\sqrt{30})/36 $$

$${{W}_{4}}=\frac{-2}{(4+1)[P_{4}^{'}({{x}_{4}}){{P}_{5}}({{x}_{4}})]}=(18+\sqrt{30})/36 $$

For n = 5
The root of the Legendre Polynomial $$ P_5(x)=0 $$ are found by solving

$$ \displaystyle P_5(x)=(63x^5-70x^3+15x)=0 $$

Where the weights, $$ w_{1}, w_{2}, w_{3}, w_{4} $$ and $$ w_{5} $$ are

$${{W}_{1}}=\frac{-2}{(5+1)[P_{5}^{'}({{x}_{1}}){{P}_{6}}({{x}_{1}})]}=\frac{128}{225}$$

$${{W}_{2}}=\frac{-2}{(5+1)[P_{5}^{'}({{x}_{2}}){{P}_{6}}({{x}_{2}})]}=\frac{322-13\sqrt{70}}{900}$$

$${{W}_{3}}=\frac{-2}{(5+1)[P_{5}^{'}({{x}_{3}}){{P}_{6}}({{x}_{3}})]}=\frac{322-13\sqrt{70}}{900}$$

$${{W}_{4}}=\frac{-2}{(5+1)[P_{5}^{'}({{x}_{4}}){{P}_{6}}({{x}_{4}})]}=\frac{322+13\sqrt{70}}{900}$$

$${{W}_{5}}=\frac{-2}{(5+1)[P_{5}^{'}({{x}_{5}}){{P}_{6}}({{x}_{5}})]}=\frac{322+13\sqrt{70}}{900}$$

As can be seen from the values computed here, the table present in lecture has been verified by myself and is in agreement with the NIST publication.

2. Solve the Following:
The following problems can be found in Fish and Belytschko pg. 91-92:

Problem 4.6: Gauss Quadrature to Obtain Numerical Integration
Use the Gauss quadrature to obtain exact values for the following integrals. Also verify by analytic integration:

(a) Integrate the following:

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

$$ \displaystyle \int_{0}^{4}{\left(x^2+1\right)dx} $$ (6.8.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the two-point Gauss quadrature where $$\displaystyle n_{gp} = 2, a=0, b=4 $$ we can determine the roots $$\displaystyle \xi_1 $$ and $$ \displaystyle \xi_2 $$ along with the weights $$\displaystyle W_1 $$ and $$\displaystyle W_2 $$ as shown above or by using the following system of equations. Either method is valid and for completeness we have decided to demonstrate both methods.


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

$$ \displaystyle \begin{align} \left[ {\begin{array}{*{20}{c}} 1&1\\ \xi_1 & \xi_2 \\ \xi_1^2 & \xi_2^2 \\ \xi_1^3 & \xi_2^3 \end{array}} \right]
 * style="width:95%" |
 * style="width:95%" |

\left[ {\begin{array}{*{20}{c}} W_1\\ W_2 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 2\\ 0\\ \frac{2}{3}\\ 0 \end{array}} \right] \end{align} $$ (6.8.13)
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad W_1 = W_2 = 1 \quad and \quad \xi_1=-\frac{\sqrt{3}}{3} \quad, \quad \xi_2 = \frac{\sqrt{3}}{3} $$ (6.8.14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we must consider the interval that $$\displaystyle x $$ is defined on. In its current form the Gauss quadrature is defined for $$\displaystyle x \in [-1,1] $$. To correct for this we will define the following relation between $$\displaystyle \xi $$ and $$ x $$.


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

$$ \displaystyle x= \frac{1}{2}(a+b)+\frac{1}{2} \xi (b-a) $$ (6.8.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle \quad x &= \frac{1}{2}(0+4)+\frac{1}{2} \xi (4-0) \\ &= \frac{1}{2}(0+4)+\frac{1}{2} \xi (4-0) \\ &= 2+ {2} \xi \end{align} $$ (6.8.16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad f(\xi)= (2+ {2} \xi)^2+1 $$ (6.8.17)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle I &= \int_{0}^{4}{\left(x^2+1\right) \underbrace{dx}_{{\color{red}\frac{1}{2}(b-a)d \xi } }} \\ &= \frac{1}{2}(b-a) \int_{-1}^{1}{\left( 2+ {2} \xi)^2+1 \right) d\xi} \\ &= \frac{1}{2}(b-a) \left[ \cancelto{w_1}((2+ {2} \cancelto{\xi_1})^2+1) + \cancelto{w_2}((2+ {2} \cancelto{\xi_2})^2+1) \right] \\  &= \frac{76}{3} \end{align} $$ (6.8.18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In this case, the approximate solutions happens to be equal to the exact solution. This is proven by the analytical solution provided by wolfframalpha

(b) $$\int_{-1}^{1}{\left(\xi^4+2\xi^2\right)d\xi} $$
To keep this consistent with other problems, we will replace $$ \xi $$ with x making this integral: $$\int_{-1}^{1}{\left(x^4+2x^2\right)dx} $$ Using the two-point Gauss quadrature where $$\displaystyle n_{gp} = 2, a=-1, b=1 $$ we can determine the roots $$\displaystyle \xi_1 $$ and $$ \displaystyle \xi_2 $$ along with the weights $$\displaystyle W_1 $$ and $$\displaystyle W_2 $$ as shown above or by using the following system of equations. Either method is valid and for completeness we have decided to demonstrate both methods.


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

$$ \displaystyle \begin{align} \left[ {\begin{array}{*{20}{c}} 1&1\\ \xi_1 & \xi_2 \\ \xi_1^2 & \xi_2^2 \\ \xi_1^3 & \xi_2^3 \end{array}} \right]
 * style="width:95%" |
 * style="width:95%" |

\left[ {\begin{array}{*{20}{c}} W_1\\ W_2 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 2\\ 0\\ \frac{2}{3}\\ 0 \end{array}} \right] \end{align} $$ (6.8.19)
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad W_1 = W_2 = 1 \quad and \quad \xi_1=-\frac{\sqrt{3}}{3} \quad, \quad \xi_2 = \frac{\sqrt{3}}{3} $$ (6.8.20)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we must consider the interval that $$\displaystyle x $$ is defined on. In its current form the Gauss quadrature is defined for $$\displaystyle x \in [-1,1] $$. To correct for this we will define the following relation between $$\displaystyle \xi $$ and $$ x $$.


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

$$ \displaystyle x= \frac{1}{2}(a+b)+\frac{1}{2} \xi (b-a) $$ (6.8.21)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle \quad x &= \frac{1}{2}(-1+1)+\frac{1}{2} \xi (1-(-1)) \\ &= \frac{1}{2}(0)+\frac{1}{2} \xi (2) \\ &= 0+ {1} \xi \end{align} $$ (6.8.22)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad f(\xi)= \xi^4+{2}\xi^2 $$ (6.8.23)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle I &= \int_{-1}^{1}{\left(x^4+2x^2\right) \underbrace{dx}_{{\color{red}\frac{1}{2}(b-a)d \xi } }} \\ &= \frac{1}{2}(b-a) \int_{-1}^{1}{\left( \xi^4 + {2} \xi^2 \right) d\xi} \\ &= \frac{1}{2}(b-a) \left[ {w_1}(\xi_1^{4}+ {2} \xi_1^{2}) + {w_2}(\xi_2^{4}+ {2} \xi_2^{2}) \right] \\ &= \frac{14}{9} \end{align} $$ (6.8.24)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In this case, the approximate solution has an error of 0.177778

This is proven by the analytical solution provided by wolfframalpha

(c) Write MATLAB Code
That utilizes the function gauss.m and performs Gauss integration. Check the solutions in (a) and (b) against the MATLAB code.

4.6a MATLAB Code:

4.6b MATLAB Code:

Problem 4.7: Three-Point Gauss Quadrature to Obtain Numerical Integration
Use the three-point Gauss quadrature to evaluate the following integrals. Compare the results to the analytic integral.

(a) Integrate the following:
$$\displaystyle \int_{-1}^{1}{\frac{\xi}{\xi^2+1}dx} $$

Using the three-point Gauss quadrature where $$\displaystyle n_{gp} = 3, a=-1, b=1 $$ we can determine the roots $$\displaystyle \xi_1 $$, $$ \displaystyle \xi_2 $$, and $$ \displaystyle \xi_3 $$ along with the weights $$\displaystyle W_1 $$, $$\displaystyle W_2 $$, and $$\displaystyle W_3 $$ as shown above.


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

$$ \displaystyle \therefore \quad W_1 = W_2 = \frac{5}{9}, W_3= \frac {8}{9} \quad and \quad \xi_1=-\sqrt{\frac{3}{5}} \quad, \quad \xi_2 = \sqrt{\frac{3}{5}}, \quad \xi_3 = 0 $$ (6.8.25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we must consider the interval that $$\displaystyle x $$ is defined on. In its current form the Gauss quadrature is defined for $$\displaystyle x \in [-1,1] $$. To correct for this we will define the following relation between $$\displaystyle \xi $$ and $$ x $$.


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

$$ \displaystyle x= \frac{1}{2}(a+b)+\frac{1}{2} \xi (b-a) $$ (6.8.26)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle \quad x &= \frac{1}{2}(-1+1)+\frac{1}{2} \xi (1-(-1)) \\ &= \frac{1}{2}(0)+\frac{1}{2} \xi (2) \\ &= 0+ {1} \xi \end{align} $$ (6.8.27)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad f(\xi)= \frac{\xi}{\xi^2+1} $$ (6.8.28)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle I &= \int_{-1}^{1}{\left(\frac{x}{x^2+1}\right) \underbrace{dx}_{{\color{red}\frac{1}{2}(b-a)d \xi } }} \\ &= \frac{1}{2}(b-a) \int_{-1}^{1}{\left( \frac{\xi}{\xi^2 + 1} \right) d\xi} \\ & = \frac{1}{2}\left( {b - a} \right)\left[ {{W_1}\left( {\frac} \right) + {W_2}\left( {\frac} \right) + {W_3}\left( {\frac} \right)} \right] \\ & = \frac{1}{2}\left( {1 - (-1)} \right)\left[ {{\frac{5}{9}}\left( {\frac} \right) + {\frac{5}{9}}\left( {\frac} \right) + {\frac{8}{9}}\left( {\frac} \right)} \right] \\ &= 0 \end{align} $$ (6.8.29)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In this case, the approximate solution has an error of 0

This is proven by the analytical solution provided by wolfframalpha

(b) Integrate the following:
$$\displaystyle \int_{-1}^{1}{cos^2\pi\zeta d\xi} $$

Using the three-point Gauss quadrature where $$\displaystyle n_{gp} = 3, a=-1, b=1 $$ we can determine the roots $$\displaystyle \xi_1 $$, $$ \displaystyle \xi_2 $$, and $$ \displaystyle \xi_3 $$ along with the weights $$\displaystyle W_1 $$, $$\displaystyle W_2 $$, and $$\displaystyle W_3 $$ as shown above.


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

$$ \displaystyle \therefore \quad W_1 = W_2 = \frac{5}{9}, W_3= \frac {8}{9} \quad and \quad \xi_1=-\sqrt{\frac{3}{5}} \quad, \quad \xi_2 = \sqrt{\frac{3}{5}}, \quad \xi_3 = 0 $$ (6.8.30)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we must consider the interval that $$\displaystyle x $$ is defined on. In its current form the Gauss quadrature is defined for $$\displaystyle x \in [-1,1] $$. To correct for this we will define the following relation between $$\displaystyle \xi $$ and $$ x $$.


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

$$ \displaystyle x= \frac{1}{2}(a+b)+\frac{1}{2} \xi (b-a) $$ (6.8.31)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle \quad x &= \frac{1}{2}(-1+1)+\frac{1}{2} \xi (1-(-1)) \\ &= \frac{1}{2}(0)+\frac{1}{2} \xi (2) \\ &= 0+ {1} \xi \end{align} $$ (6.8.32)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad f(\xi)= {\cos ^2}(\pi \xi ) $$ (6.8.33)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle I &= \int_{-1}^{1}{\left({\cos ^2}(\pi \xi )\right) \underbrace{dx}_{{\color{red}\frac{1}{2}(b-a)d \xi } }} \\ &= \frac{1}{2}(b-a) \int_{-1}^{1}{\left( {\cos ^2}(\pi \xi ) \right) d\xi} \\ & = \frac{1}{2}\left( {b - a} \right)\left[ {{W_1}\left( {{{\cos }^2}\left( {\pi {\xi _1}} \right)} \right) + {W_2}\left( {{{\cos }^2}\left( {\pi {\xi _2}} \right)} \right) + {W_3}\left( {{{\cos }^2}\left( {\pi {\xi _3}} \right)} \right)} \right] \\ & = \frac{1}{2}\left( {1 - (-1)} \right)\left[ {\frac{5}{9}\left( {{{\cos }^2}\left( {\pi \cdot  - \sqrt {\frac{3}{5}} } \right)} \right) + \frac{5}{9}\left( {{{\cos }^2}\left( {\pi \sqrt {\frac{3}{5}} } \right)} \right) + \frac{8}{9}\left( {{{\cos }^2}\left( {\pi 0} \right)} \right)} \right] \\ &= 1.52996 \end{align} $$ (6.8.34)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In this case, the approximate solution has an error of 0.52996

This is proven by the analytical solution provided by wolfframalpha

(c) Write MATLAB Code
That utilizes the function gauss.m and performs Gauss integration. Check the solutions in (a) and (b) against the MATLAB code

4.7a MATLAB Code:

4.7b MATLAB Code:

Problem 4.8: n-Point Gauss Quadrature Comparison for $$ n=1,2,3 $$
The integral $$ \displaystyle \int_{-1}^{1}{\left(3\xi^3+2\right)d\xi} $$ can be integrated exactly using the two-point Gauss quadrature. How is the accuracy effected if

(a) One-Point Quadrature is Employed
Using the one-point Gauss quadrature where $$\displaystyle n_{gp} = 1, a=-1, b=1 $$ we can determine the roots $$\displaystyle \xi_1 $$ along with the weights $$\displaystyle W_1 $$ as shown above.


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

$$ \displaystyle \therefore \quad W_1 = 2 \quad and \quad \xi_1=0 \quad $$ (6.8.35)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we must consider the interval that $$\displaystyle x $$ is defined on. In its current form the Gauss quadrature is defined for $$\displaystyle x \in [-1,1] $$. To correct for this we will define the following relation between $$\displaystyle \xi $$ and $$ x $$.


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

$$ \displaystyle x= \frac{1}{2}(a+b)+\frac{1}{2} \xi (b-a) $$ (6.8.36)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle \quad x &= \frac{1}{2}(-1+1)+\frac{1}{2} \xi (1-(-1)) \\ &= \frac{1}{2}(0)+\frac{1}{2} \xi (2) \\ &= 0+ {1} \xi \end{align} $$ (6.8.37)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad f(\xi)= (3\xi^{3}+2)$$ (6.8.38)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle I &= \int_{-1}^{1}{\left(3\xi^{3}+2\right) \underbrace{dx}_{{\color{red}\frac{1}{2}(b-a)d \xi } }} \\ &= \frac{1}{2}(b-a) \int_{-1}^{1}{\left( 3\xi^{3}+2 \right) d\xi} \\ & = \frac{1}{2}\left( {b - a} \right)\left[ {{W_1}\left( {3{\xi _1}^3 + 2} \right)} \right] \\ & = \frac{1}{2}\left( {1 - (-1)} \right)\left[ {2\left( {3{{(0)}^3} + 2} \right)} \right] \\ &= 4 \end{align} $$ (6.8.39)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In this case, the approximate solution has an error of 0

This is proven by the analytical solution provided by wolfframalpha

(b) Three-Point Quadrature is Employed
Using the three-point Gauss quadrature where $$\displaystyle n_{gp} = 3, a=-1, b=1 $$ we can determine the roots $$\displaystyle \xi_1 $$, $$ \displaystyle \xi_2 $$, and $$ \displaystyle \xi_3 $$ along with the weights $$\displaystyle W_1 $$, $$\displaystyle W_2 $$, and $$\displaystyle W_3 $$ as shown above.


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

$$ \displaystyle \therefore \quad W_1 = W_2 = \frac{5}{9}, W_3= \frac {8}{9} \quad and \quad \xi_1=-\sqrt{\frac{3}{5}} \quad, \quad \xi_2 = \sqrt{\frac{3}{5}}, \quad \xi_3 = 0 $$ (6.8.40)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we must consider the interval that $$\displaystyle x $$ is defined on. In its current form the Gauss quadrature is defined for $$\displaystyle x \in [-1,1] $$. To correct for this we will define the following relation between $$\displaystyle \xi $$ and $$ x $$.


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

$$ \displaystyle x= \frac{1}{2}(a+b)+\frac{1}{2} \xi (b-a) $$ (6.8.41)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle \quad x &= \frac{1}{2}(-1+1)+\frac{1}{2} \xi (1-(-1)) \\ &= \frac{1}{2}(0)+\frac{1}{2} \xi (2) \\ &= 0+ {1} \xi \end{align} $$ (6.8.42)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \therefore \quad f(\xi)= (3\xi^{3}+2)$$ (6.8.43)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \begin{align} \displaystyle I &= \int_{-1}^{1}{\left(3\xi^{3}+2\right) \underbrace{dx}_{{\color{red}\frac{1}{2}(b-a)d \xi } }} \\ &= \frac{1}{2}(b-a) \int_{-1}^{1}{\left( 3\xi^{3}+2 \right) d\xi} \\ & = \frac{1}{2}\left( {b - a} \right)\left[ {{W_1}\left( {3{\xi _1}^3 + 2} \right) + {W_2}\left( {3{\xi _2}^3 + 2} \right) + {W_3}\left( {3{\xi _3}^3 + 2} \right)} \right] \\ & = \frac{1}{2}\left( {1 - (-1)} \right)\left[ {\frac{5}{9}\left( {3{{\left( { - \sqrt {\frac{3}{5}} } \right)}^3} + 2} \right) + \frac{5}{9}\left( {3{{\left( {\sqrt {\frac{3}{5}} } \right)}^3} + 2} \right) + \frac{8}{9}\left( {3{{(0)}^3} + 2} \right)} \right] \\ &= 4 \end{align} $$ (6.8.44)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In this case, the approximate solution has an error of 0

This is proven by the analytical solution provided by wolfframalpha

(c) Use MATLAB Code
Check the solutions in (a) and (b) against the MATLAB code

4.8a MATLAB Code:

4.8b MATLAB Code:

=Problem 6.9 Fish and Belytschko Chapter 8 Code=

3. Solve the following heat conduction problem
Consider a chimney constructed of two isotropic materials: dense concrete $$ (k = 2.0 W^{\circ}C^{-1} ) $$ and bricks $$ (k = 0.9 W^{\circ}C^{-1} ) $$. The temperature of the hot gases on the inside surface of the chimney is $$140^{\circ}C$$, whereas the outside is exposed to the surrounding air, which is at $$ T= 10C^{\circ}C^{-1} $$. The dimensions of the chimney (in meters) are shown below. For the analysis, exploit the symmetry and consider 1/8 of the chimney cross-sectional area. Consider a mesh of eight elements as shown below. Determine the temperature and flux in the two materials.

Analyze the problem with $$2x2$$, $$4x4$$ and $$8x8$$ quadrilateral elements for 1/8 of the problem domain. $$A2x2$$ finite element mesh is shown in Figure 8.23. Symmetry implies insulated boundary conditions on edges AD and BC. Note that element boundaries have to coincide with the interface between the concrete and bricks.

Solution
Matlab Code 1: Heat2d.m

adfghasdf

Matlab Code 2:  Include_flags.m

qwrgasg

Matlab Code 3:  preprocessor.m

Matlab Code 4:  input_file_16ele.m

Matlab Code 5: mesh2d.m

Matlab Code 6:  plotmesh.m

Matlab Code 7:  setup_ID_LM.m

Matlab Code 8:  heat2Delem.m

Matlab Code 9:  gauss.m

Matlab Code 10:  BmatHeat2D.m

Matlab Code 11:  assembly.m

Matlab Code 12:  src_and_flux.m

Matlab Code 13:  solvedr.m Matlab Code 14:   postprocessor.m

Matlab Code 15:  postprocessor.m

Here is a link to the pseudo code that explains the processes behind the above matlab codes,

1. Graphs for 8.3
The above figure corresponds to Figure8.9 on page 200 of Fish and Belytschko.

The above figure corresponds to Figure8.10 on page 201 of Fish and Belytschko.

The above figure corresponds to Figure8.12 on page 202 of Fish and Belytschko.

3. Graphs and matlab code for problem 8.6
The first three graphs show the actual set up of the mesh for each given number of elements.

The temperature distribution throughout the chimney is shown below for each case. As we increase the number of elements in our mesh we see improved accuracy of temperature at places that are between nodes. The more elements and nodes in a mesh the better overall picture we will see.

These three graphs show the flux through the chimney's section. We can see that the flux is from the center to the outer edge. The flux does not go out of the side edges. This is to be expected because of the way in which we defined our problem. We specified that because of symmetry we have insulated boundary conditions along the vertical edges. Therefore flux should not pass through these boundaries.

Below you will find the relevant codes that were needed to solve this problem. The input file defining the boundary conditions for each, then the definition of the mesh for each case, and finally a preprocessing code. The preprocessing code is written such that the user can comment or uncomment the desired input code depending on which mesh they wish to use.

=Problem 6.10:Heat Problem with New Boundary Condition=

Given: The Governing Equation and Boundary Conditions
Governing equation was derived in meeting 32-4.
 * $$\displaystyle Q_{1}$$:Heat flow into w throgh $$\displaystyle \partial w$$
 * $$\displaystyle Q_{2}$$:Heat generated in w by heat source.
 * $$\displaystyle Q_{3}$$:Heat in w due to change in term u.


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

$$  \displaystyle Q_{1}=-\int_{\partial w}\mathbf{q}\mathbf{n}d(\partial w)=\int_{w}div\mathbf{q}dw $$     (6.10.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle Q_{2}=\int_{w}f(\mathbf{x},t)dw $$     (6.10.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle Q_{3}=\int_{w}\rho c\frac{\partial u}{\partial t}dw $$     (6.10.3)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle Q_{1}+Q_{2}=Q_{3} $$     (6.10.4)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \mathbf{q}=-\mathbf{K}grad(u) $$     (6.10.5)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

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

Boundry Conditions

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

$$  \displaystyle u(\mathbf{x},t)]_{\Gamma _{g}}=g(\mathbf{x},t) $$     (6.10.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle (q.n)]_{\Gamma _{h}}=h(\mathbf{x},t) $$     (6.10.8)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle (q.n)]_{\Gamma _{H}}=H(u-u_{\infty }) $$     (6.10.9)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

WRF

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

$$  \displaystyle \int_{\Omega }w[div(Kgrad(u))+f-\rho c\frac{\partial u}{\partial t}]d\Omega =0 $$     (6.10.10)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} K=  & \int_{\Omega }wdiv(Kgrad(u))d\Omega =\frac{\partial }{\partial x_{i}}(K_{ij}\frac{\partial u}{\partial x_{j}})\\ & \int_{\Omega }n_{i}wK_{ij}\frac{\partial u}{\partial x_{j}}d(\partial\Omega)-\int_{\Omega }\frac{\partial w}{\partial x_{i}}K_{ij}\frac{\partial u}{\partial x_{j}}d(\Omega) \end{align} $$     (6.10.11)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle K=-\int_{\Gamma _{g}}wnqd\Gamma _{g}-\int_{\Gamma _{h}}wnqd\Gamma _{h}-\int_{\Gamma _{H}}wnqd\Gamma _{H}+\int_{\Omega }div(w).K.div(u)d\Omega $$     (6.10.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

We deliberately chose w i.e. i vanishes on essential boundry. So the first term vanishes too. In the first term we have unknown flux that we can not deal with right now. But other fluxes are defined.We can write our continuous form of components as:


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

$$  \displaystyle \tilde{m}(w,u)=\int_{\Omega}w\rho c\frac{\partial u}{\partial t}d\Omega $$     (6.10.13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \tilde{k} (w,u)=\int_{\Omega}div(w)Kdiv(u)d\Omega+\int_{\Gamma_{H}}wHud\Gamma_{H} $$     (6.10.14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{f}(w,u)=\int_{\Omega}wfd\Omega-\int_{\Gamma_{h}}whd\Gamma_{h}+\int_{\Gamma_{H}}wHu_{\infty }d\Gamma_{H} $$     (6.10.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

!Note this time f depeneds on both w and u.

Discrete WF

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

$$  \displaystyle w^{app}=\sum_{I=1}^{n}N(x,y)_{I}w_{I} $$     (6.10.16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle u^{app}=\sum_{J=1}^{n}N(x,y)_{J}d_{J} $$     (6.10.17)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle div(u)=div(N(x,y)_{I})d_{I}=B(x,y)_{I}d_{I} $$     (6.10.18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{m}(w,u)+\tilde{k}(w,u)=\tilde{f}(w,u) $$     (6.10.19) Now lets substitute everything in the equation above.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \sum_{I=1}^{n}w^{T}\left [\sum_{J=1}^{n}\left \{ \int_{\Omega }\rho cN_{I}N_{J}d_{I}^{s}d\Omega+\int_{\Omega}B_{I}KB_{J}d_{I}d\Omega -\int_{\Omega}N_{I}fd\Omega+\int_{\Gamma_{h}}N_{I}hd\Gamma_{h}+\int_{\Omega}N_{I}N_{J}(\Gamma_{H})H(d(\Gamma_{H})-d_{\infty })\right \} \right ]=0 $$     (6.10.20)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{d}}=\begin{bmatrix} \bar{de}\\ df
 * style="width:95%" |
 * style="width:95%" |

\end{bmatrix} $$     (6.10.21)
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{w}}=\begin{bmatrix} 0\\ wf
 * style="width:95%" |
 * style="width:95%" |

\end{bmatrix} $$     (6.10.22)
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{K}}=\begin{bmatrix} K_{E} &K_{EF} \\ K_{EF}^{T} & K_{F} \end{bmatrix} $$     (6.10.23)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle\begin{align} & K_{E}=\int_{\Omega }B_{1}.K.B_{1}d\Omega&\\ & K_{EF}=\int_{\Omega }B_{1}.K.B_{J}d\Omega& j=2,n\\ & K_{FF}=\int_{\Omega }B_{I}.K.B_{J}d\Omega& i,j=2,n
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$     (6.10.24)
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{M}}=\begin{bmatrix} M_{E} &M_{EF} \\ M_{EF}^{T} & M_{F} \end{bmatrix} $$     (6.10.25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle\begin{align} & M_{E}=\int_{\Omega }N_{1}\rho cN_{1}d\Omega&\\ & M_{EF}=\int_{\Omega }N_{1}\rho cN_{J}d\Omega& j=2,n\\ & M_{FF}=\int_{\Omega }N_{I}\rho cN_{J}d\Omega& i,j=2,n
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$     (6.10.26)
 * <p style="text-align:right">
 * }

! When the convection coefficient is very large the temperature $$\displaystyle T_{\infty }$$ is immediately felt at on $$\displaystyle \Gamma_{H}$$ and the essential boundary condition is again enforced as a limiting case of natural boundary condition .We can use penalty method.


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

$$  \displaystyle \begin{matrix} d & Matrix \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }




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

$$  \displaystyle \begin{matrix} Temperature & Contours \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }



=Contributing Members & Referenced Lecture=