text:
User:Egm6341.s11.team4/HW6

=Problem 6.1 Find $$ P_{6} $$ and $$ P_{7} $$ in the HOTRE=

Problem Statement
Given on Lecture Slide 31-1

Given
$$ E=[P_{2}(t)g^{(1)}(t)+P_{4}(t)g^{(3)}(t)]_{-1}^{+1}-\int\limits_{-1}^{+1} P_{5}(t)g^{(5)}(t)dt $$

Objective
1. Perform steps 4ab to find $$(P_{6}, P_{7}) $$ and E

2.1 Find $$ t_{k}(x) $$

2.2 Find $$ d_{1}, d_{2} $$ and $$ d_{3} $$

Part 1.1:$$(P_{6},P_{7}) $$
We will start with the values we obtained for polynomials 4 and 5 in class

$$ P_{4}(t)=c_{1}(\frac{t^{4}}{4!})+c_{3}(\frac{t^{2}}{2!})+c_{5}$$

$$ P_{5}(t)=c_{1}(\frac{t^{5}}{5!})+c_{3}(\frac{t^{3}}{3!})+c_{5}$$

 * $$ c_{1}=-1 $$
 * $$ c_{3}=\frac{1}{6} $$
 * $$ c_{5}=\frac{-7}{360} $$

We will now perform step 4a as was done with 3a in class using integration by parts. This was performed in class on Lecture slides 30-5 and 30-6

$$ E=[P_2(t)g^{(1)}(t)+P_4(t)g^{(3)}(t)]_{-1}^{1}+ [P_6(t)g^{(5)}(t)]_{-1}^{1}-\int_{-1}^{1}P_6(t)g^(4)(t)dt$$

$$ P_6(t)=\int P_5(t)=c_1\frac{t^6}{6!}+c_3\frac{t^4}{4!}+c_5\frac{t^2}{2!}+c_7 $$

Now we must perform step 4b, once again using integration by parts.

$$ E=[P_2(t)g^{(1)}(t)+P_4(t)g^{(3)}(t)+P_{6}(t)g^{(5)}(t)]_{-1}^{1} +[P_7g^{(6)}(t)]_{-1}^{1} -int_{-1}^{1} P_6(t)g^{(5)}(t)dt$$

$$ P_7(t)=\int P_6(t)dt=c_{1}\frac{t^7}{7!}+c_3\frac{t^5}{5!}+c_5\frac{t^{3}}{3!}+c_{7}t+c_8 $$ If we now proceed such that

$$ P_5(0)=0 $$     $$ P_5(\pm 1)=0 $$

we can determine the higher powered coefficents.

The equation at t=0 becomes

$$ P_7(0)=c_{1}\frac{0^7}{7!}+c_3\frac{0^5}{5!}+c_5\frac{0^{3}}{3!}+c_{7}0+c_8 $$ $$ c_8=0 $$

and at $$ t=\pm 1 $$

$$ P_7(\pm 1)=c_{1}\frac{1^7}{7!}+c_3\frac{1^5}{5!}+c_5\frac{1^{3}}{3!}+c_{7}1+c_8 $$

and now through substitution of $$ c_8=0 $$ we obtain

$$ P_7(\pm 1)=0=c_{1}\frac{1^7}{7!}+c_3\frac{1^5}{5!}+c_5\frac{1^{3}}{3!}+c_{7}1 $$

Through substituting in the previously determined values for our other coefficients, we obtain

$$ c_7=\frac{31}{15120} $$

Now if we re-express our polynomials with our obtained coefficients.

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

 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ P_6(t)=\frac{-t^6}{6!}+\frac{t^4}{144}-\frac{7t^2}{720}+\frac{31}{151210} $$

 * }

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

 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ P_7(t)=\frac{-t^7}{7!}+\frac{t^5}{720}-\frac{7t^3}{2160}+\frac{31}{151210}t $$

 * }

And our E value is

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

 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$E=[P_2(t)g^{(1)}(t)+P_4(t)g^{(3)}(t)+P_6(t)g^{(5)}(t)]_{-1}^{1}-\int_{-1}^{1}P_7(t)g^{(7)}(t)dt $$

 * }

Part 2.1: Find $$ t_k(x) $$
Starting from the given equation.

$$ x(t) := \frac{x_k+x_{k+1}}{2}+t\frac{h}{2} $$

$$ t \in [-1, 1] $$

$$ t_k $$ can be found from algebraic manipulation.

$$ x(t_k) := \frac{x_k+x_{k+1}}{2}+t_k\frac{h}{2} $$

$$ t_k=\frac{2}{h}(x(t_k)-\frac{x_k+x_{k+1}}{2}) $$

$$ t_k=\frac{2}{h}(x-\frac{x_k+x_{k+1}}{2}) $$

As expected, this relation yields a value of -1 when evaluated at $$ x_k $$ and a value of 1 when evaluated at $$ x_{k+1} $$. This is due to the definition of h ($$ h:=x_{k+1}-x_k $$).

Part 2.2: Find $$ d_1,d_2, $$ and $$ d_3 $$
These can all be obtained from their definitions (which are all composed entirely of terms we've already developed.

$$ d_1=\overline{d}_{2.1}=\overline{d}_{2}=\frac{P_2(1)}{2^2}$$

$$ d_2=\overline{d}_{4}=\frac{P_4(1)}{2^4} $$

$$ d_3=\overline{d}_{6}=\frac{P_6(1)}{2^6} $$

Since we know our values of $$ P_2,P_4 $$ and $$ P_6 $$ we can evaluate these expressions

$$ P_2(1)=-\frac{1^2}{2!}+\frac{1}{6}=-\frac{1}{3} $$

$$ P_4(1)=-\frac{1^4}{4!}+\frac{1^2}{12}+\frac{-7}{360}=\frac{1}{45} $$

$$ P_6(1)=\frac{-1^6}{6!}+\frac{1^4}{144}-\frac{7*1^2}{720}+\frac{31}{151210}=-\frac{1}{945} $$

Now substituting these values into our expressions for our d values we obtain

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

 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ d_1=-\frac{1}{12} $$

 * }

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

 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ d_2=\frac{1}{720} $$

 * }

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

 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ d_3=-\frac{1}{30240} $$

 * }

This problem was solved by Brendan Mahon

Problem Statement
Refer to lecture notes [[media:nm1.s11.mtg31.djvu|31-3]] for more detailed description.

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

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

\begin{align} & E_{n}^{T}=I-{{T}_{0}}(n) \\ & =\sum\limits_{r=1}^{l}{{{h}^{2r}}{{{\bar{d}}}_{2r}}[{{f}^{2r-1}}(b)-{{f}^{2r-1}}(a)]-{{(\frac{h}{2})}^{2r}}\sum\limits_{k=0}^{n-1}{\int\limits_^{{{P}_{2l}}({{t}_{k}}(x)){{f}^{(2l)}}(x)dx}}} \\ \end{align}

$$ And :{| style="width:100%" border="0" $$\displaystyle
 * }
 * style="width:95%" |
 * style="width:95%" |

{{\bar{d}}_{2r}}=\frac{{{P}_{2r}}(1)}

$$
 * }.

Solution
In steps 2ab, we already got
 * {| style="width:100%" border="0"

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

\begin{align} & E=[{{P}_{2}}(t){{g}^{(1)}}(t)]_{-1}^{+1}-\int\limits_{-1}^{+1}{{{P}_{2}}(t){{g}^{(2)}}(t)dt} \\ & {{P}_{2}}(t)=-\frac{2!}+\frac{1}{6} \\ & {{P}_{3}}(t)=-\frac{3!}+\frac{1}{6}t \\ \end{align}

$$ Substitute $$ g(t)=f(x(t))$$, $${{g}^{(i)}}(t)={{(\frac{h}{2})}^{i}}{{f}^{(i)}}(x(t))$$, $$x=\frac{h}{2}t+\frac{{{x}_{k}}+{{x}_{k+1}}}{2} $$ and $$dx=\frac{h}{2}dt $$ into above equation, we have
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & E=[{{P}_{2}}(t){{g}^{(1)}}(t)]_{-1}^{+1}-\underset{-1}{\overset{+1}{\mathop \int }}\,{{P}_{2}}(t){{g}^{(2)}}(t)dt \\ & =[{{P}_{2}}(+1){{(\frac{h}{2})}^{1}}{{f}^{(1)}}({{x}_{k+1}})-{{P}_{2}}(-1){{(\frac{h}{2})}^{1}}{{f}^{(1)}}({{x}_{k}})]-\underset{\overset{\mathop \int }}\,{{P}_{2}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{(\frac{h}{2})}^{2}}{{f}^{(2)}}(x)\frac{2}{h}dx \\ & ={{(\frac{h}{2})}^{1}}{{P}_{2}}(1)[{{f}^{(1)}}({{x}_{k+1}})-{{f}^{(1)}}({{x}_{k}})]-{{(\frac{h}{2})}^{1}}\underset{\overset{\mathop \int }}\,{{P}_{2}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(2)}}(x)dx \\ \end{align}

$$ Hence,
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & E_{1}^{T}=\frac{h}{2}\sum\limits_{k=0}^{0}{E} \\ & ={{(\frac{h}{2})}^{2}}{{P}_{2}}(1)\sum\limits_{k=0}^{0}{[{{f}^{(1)}}({{x}_{k+1}})-{{f}^{(1)}}({{x}_{k}})]}-{{(\frac{h}{2})}^{2}}\sum\limits_{k=0}^{0}{\underset{\overset{\mathop \int }}\,{{P}_{2}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(2)}}(x)dx} \\ & ={{(\frac{h}{2})}^{2}}{{P}_{2}}(1)[{{f}^{(1)}}(b)-{{f}^{(1)}}(a)]-{{(\frac{h}{2})}^{2}}\sum\limits_{k=0}^{0}{\underset{\overset{\mathop \int }}\,{{P}_{2}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(2)}}(x)dx} \\ \end{align}

$$ And in step 3ab, we also obtained
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & E=[{{P}_{2}}(t){{g}^{(1)}}(t)+{{P}_{4}}(t){{g}^{(3)}}(t)]_{-1}^{+1}-\int\limits_{-1}^{+1}{{{P}_{4}}(t){{g}^{(4)}}(t)dt} \\ & {{P}_{4}}(t)=-\frac{4!}+\frac{1}{6}\frac{2!}-\frac{7}{360} \\ & {{P}_{5}}(t)=-\frac{5!}+\frac{1}{6}\frac{3!}-\frac{7}{360}t \\ \end{align}

$$ Using the same technique for step2ab, we have
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & E=[{{P}_{2}}(t){{g}^{(1)}}(t)+{{P}_{4}}(t){{g}^{(3)}}(t)]_{-1}^{+1}-\underset{-1}{\overset{+1}{\mathop \int }}\,{{P}_{4}}(t){{g}^{(4)}}(t)dt \\ & ={{(\frac{h}{2})}^{1}}{{P}_{2}}(1)[{{f}^{(1)}}({{x}_{k+1}})-{{f}^{(1)}}({{x}_{k}})]+{{(\frac{h}{2})}^{3}}{{P}_{4}}(1)[{{f}^{(3)}}({{x}_{k+1}})-{{f}^{(3)}}({{x}_{k}})]-\underset{\overset{\mathop \int }}\,{{P}_{4}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{(\frac{h}{2})}^{4}}{{f}^{(4)}}(x)\frac{2}{h}dx \\ & ={{(\frac{h}{2})}^{1}}{{P}_{2}}(1)[{{f}^{(1)}}({{x}_{k+1}})-{{f}^{(1)}}({{x}_{k}})]+{{(\frac{h}{2})}^{3}}{{P}_{4}}(1)[{{f}^{(3)}}({{x}_{k+1}})-{{f}^{(3)}}({{x}_{k}})]-{{(\frac{h}{2})}^{3}}\underset{\overset{\mathop \int }}\,{{P}_{4}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(4)}}(x)dx \\ \end{align}

$$ Hence,
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & E_{2}^{T}=\frac{h}{2}\sum\limits_{k=0}^{1}{E} \\ & ={{(\frac{h}{2})}^{2}}{{P}_{2}}(1)\sum\limits_{k=0}^{1}{[{{f}^{(1)}}({{x}_{k+1}})-{{f}^{(1)}}({{x}_{k}})]}+{{(\frac{h}{2})}^{4}}{{P}_{4}}(1)\sum\limits_{k=0}^{1}{[{{f}^{(3)}}({{x}_{k+1}})-{{f}^{(3)}}({{x}_{k}})]}-{{(\frac{h}{2})}^{4}}\sum\limits_{k=0}^{1}{\underset{\overset{\mathop \int }}\,{{P}_{4}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(4)}}(x)dx} \\ & ={{(\frac{h}{2})}^{2}}{{P}_{2}}(1)[{{f}^{(1)}}(b)-{{f}^{(1)}}(a)]+{{(\frac{h}{2})}^{4}}{{P}_{4}}(1)[{{f}^{(3)}}(b)-{{f}^{(3)}}(a)]-{{(\frac{h}{2})}^{4}}\sum\limits_{k=0}^{1}{\underset{\overset{\mathop \int }}\,{{P}_{4}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(4)}}(x)dx} \\ & =\sum\limits_{r=1}^{2}{{{(\frac{h}{2})}^{2r}}{{P}_{2r}}(1)[{{f}^{(2r-1)}}(b)-{{f}^{(2r-1)}}(a)]-{{(\frac{h}{2})}^{4}}\sum\limits_{k=0}^{1}{\underset{\overset{\mathop \int }}\,{{P}_{4}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(4)}}(x)dx}} \\ \end{align}

$$ Using the same technique for step $$n$$ab, we can obtain
 * }
 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

\begin{align} & E_{n}^{T}=\frac{h}{2}\sum\limits_{k=0}^{n-1}{E} \\ & =\sum\limits_{r=1}^{l}{{{(\frac{h}{2})}^{2r}}{{P}_{2r}}(1)[{{f}^{(2r-1)}}(b)-{{f}^{(2r-1)}}(a)]-{{(\frac{h}{2})}^{2l}}\sum\limits_{k=0}^{1}{\underset{\overset{\mathop \int }}\,{{P}_{2l}}(\frac{x-{{x}_{k+1}}+x-{{x}_{k}}}{h}){{f}^{(2l)}}(x)dx}} \\ \end{align}

$$ This is also the higher order trapezoidal rule error theorem.
 * }

Problem was solved by Hailong Chen'

=Problem 6.3 Verify Bernoulli Number and Comparison with HOTRE=

Problem Statement
Refer to lecture notes [[media:nm1.s11.mtg31.djvu|31-3]] for more detailed description.

Given
Bernoulli numbers:

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

$$\displaystyle \frac{x}{{{e}^{x}}-1}=\sum\limits_{r=0}^{\infty }{\underbrace{\frac{(2r)!}}_{-{{{\bar{d}}}_{2r}}}{{x}^{2r}}} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

1).
Verify $${{\bar{d}}_{2}},{{\bar{d}}_{4}},{{\bar{d}}_{6}}$$

2).
Find $${{\bar{d}}_{8}},{{\bar{d}}_{10}}$$ using above given expression, compare with $${{\bar{d}}_{2r}}=\frac{{{P}_{2r}}(1)}$$, where $${{P}_{2r}}(1)$$ is the polynomial of HOTRE evaluated at $$x=1$$.

Solution
Using Taylor Series to expand $$\frac{1}{{{e}^{x}}-1}$$ at $${{x}_{0}}=0$$ we have
 * {| style="width:100%" border="0"

$$\displaystyle \frac{1}{{{e}^{x}}-1}=\frac{1}{x}-\frac{1}{2}+\frac{x}{12}-\frac{720}+\frac{30240}-\frac{1209600}+\frac{4700160}-\frac{691{{x}^{11}}}{130767436800}+\cdots $$
 * style="width:95%" |
 * style="width:95%" |
 * }

WolframAlpha

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

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

\begin{align} & \frac{x}{{{e}^{x}}-1}=1-\frac{x}{2}+\frac{12}-\frac{720}+\frac{30240}-\frac{1209600}+\frac{4700160}-\frac{691{{x}^{12}}}{130767436800}+\cdots \\ & ={{B}_{0}}+{{B}_{1}}x+\frac{2!}{{x}^{2}}+\frac{4!}{{x}^{4}}+\frac{6!}{{x}^{6}}+\frac{8!}{{x}^{8}}+\frac{10!}{{x}^{10}}+\frac{12!}{{x}^{12}}+\cdots \\ \end{align}

$$ Comparing above equation, we obtain Bernoulli numbers:
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & {{B}_{2}}=\frac{1}{6} \\ & {{B}_{4}}=-\frac{1}{30} \\ & {{B}_{6}}=\frac{1}{42} \\ & {{B}_{8}}=-\frac{1}{30} \\ & {{B}_{10}}=\frac{5}{66} \\ & \vdots \\ \end{align}

$$ Hence we have :{| style="width:100%" border="0" $$\displaystyle
 * }
 * style="width:95%" |
 * style="width:95%" |

\begin{align} & {{{\bar{d}}}_{2}}=-\frac{1}{12} \\ & {{{\bar{d}}}_{4}}=\frac{1}{720} \\ & {{{\bar{d}}}_{6}}=-\frac{1}{30240} \\ & {{{\bar{d}}}_{8}}=\frac{1}{1209600} \\ & {{{\bar{d}}}_{10}}=-\frac{1}{\text{4700160}} \\ \end{align}

$$ From higher order trapezoidal rule error theorem, we already know that
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & {{P}_{2}}=-\frac{2!}+\frac{1}{6} \\ & {{P}_{4}}=-\frac{4!}+\frac{1}{6}\frac{2!}-\frac{7}{360} \\ & {{P}_{6}}=-\frac{6!}+\frac{1}{6}\frac{4!}-\frac{7}{360}\frac{2!}+\frac{31}{15120} \\ & {{P}_{8}}=-\frac{8!}+\frac{1}{6}\frac{6!}-\frac{7}{360}\frac{4!}+\frac{31}{15120}\frac{2!}-\frac{5}{23811} \\ & {{P}_{10}}=-\frac{10!}+\frac{1}{6}\frac{8!}-\frac{7}{360}\frac{6!}+\frac{31}{15120}\frac{4!}-\frac{5}{23811}\frac{2!}+\frac{11}{62803} \\ \end{align}

$$ So, we can get
 * }
 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

\begin{align} & {{{\bar{d}}}_{2}}=\frac{{{P}_{2}}(1)}=\frac{-\frac{1}{2!}+\frac{1}{6}}{4}=-\frac{1}{12} \\ & {{{\bar{d}}}_{4}}=\frac{{{P}_{4}}(1)}=\frac{-\frac{1}{4!}+\frac{1}{6}\frac{1}{2!}-\frac{7}{360}}{16}=\frac{1}{720} \\ & {{{\bar{d}}}_{6}}=\frac{{{P}_{6}}(1)}=\frac{-\frac{1}{6!}+\frac{1}{6}\frac{1}{4!}-\frac{7}{360}\frac{1}{2!}+\frac{31}{15120}}{64}=-\frac{1}{30240} \\ & {{{\bar{d}}}_{8}}=\frac{{{P}_{8}}(1)}=\frac{-\frac{1}{8!}+\frac{1}{6}\frac{1}{6!}-\frac{7}{360}\frac{1}{4!}+\frac{31}{15120}\frac{1}{2!}-\frac{5}{23811}}{256}=\frac{1}{1209600} \\ & {{{\bar{d}}}_{10}}=\frac{{{P}_{10}}(1)}=\frac{-\frac{1}{10!}+\frac{1}{6}\frac{1}{8!}-\frac{7}{360}\frac{1}{6!}+\frac{31}{15120}\frac{1}{4!}-\frac{5}{23811}\frac{1}{2!}+\frac{11}{62803}}{1024}=-\frac{1}{\text{4700160}} \\ \end{align}

$$
 * }

Problem was solved by Hailong Chen

=Problem 6.4 Redo proof of HOTRE by trying to cancel terms with odd order of derivatives of g=

Given
From lecture note 30-2.

Let, $$ \ P_1(t) := -t $$

Then,

Objectives
1) Redo proof of Higher Order Trapezoidal Rule Error by trying to cancel terms with odd order of derivatives of g.

Solution
From the given equation.

Do integration by parts to eq(4.1)

To cancel terms with odd order of derivative of g, $$ \ P_2(1) = P_2(-1) = 0$$ in eq(4.2).

Do integration by parts again to eq(4.2).

Do integration by parts again to eq(4.4).

To cancel terms with odd order of derivative of g, $$ \ P_4(1) = P_4(-1) = 0$$ in eq(4.5).

Do integration by parts again to eq(4.6).

Do integration by parts again to eq(4.7).

Insert and combine eq(4.4), eq(4.5), eq(4.7) and eq(4.8) into eq(4.2).

After more integration by parts to eq(4.9).

Since $$ \ P_{2r+1}(1) = -P_{2r+1}(-1) $$ is odd function.

From the given equation.

Objective
Use the recurrence relation, equation [[media:Nm1.s11.mtg32.djvu|(1)p.32-2]]

to find $$\displaystyle \;\; (p_2,p_3),\;(p_4,p_5),\;(p_6,p_7),\;(p_8,p_9)$$.

Solution
From equations [[media:Nm1.s11.mtg32.djvu|(1)p.32-1]] and [[media:Nm1.s11.mtg32.djvu|(3)-(4)p.32-1]] the polynomials $$\displaystyle p_{2i}(t)$$ and $$\displaystyle p_{2i+1}(t)$$ take the following form:

Find $$(p_2,p_3)$$
Here, $$\displaystyle i=1$$. Equation (5.2) becomes

Using equations (5.3) and (5.4) to reconstruct $$\displaystyle p_2$$ and $$\displaystyle p_3$$

Find $$(p_4,p_5)$$
Now, $$\displaystyle i=2$$. Equation (5.2) becomes

Using equations (5.3) and (5.4) to reconstruct $$\displaystyle p_4$$ and $$\displaystyle p_5$$

Find $$(p_6,p_7)$$
Now, $$\displaystyle i=3$$. Equation (5.2) becomes

Using equations (5.3) and (5.4) to reconstruct $$\displaystyle p_6$$ and $$\displaystyle p_7$$

Find $$(p_8,p_9)$$
Now, $$\displaystyle i=4$$. Equation (5.2) becomes

Using equations (5.3) and (5.4) to reconstruct $$\displaystyle p_8$$ and $$\displaystyle p_9$$

Solved by William Kurth.

=Problem 6.6: Kessler's Trapezoidal Rule Error Matlab Code=

Refer to lecture slide [[media:nm1.s11.mtg32.djvu|32-3]] for the problem statement.

Given
Kessler's Matlab code given in the following link:

Kessler's Trapezoidal Rule Error Matlab Code

Objective
a.) Use (P2i,P2i+1), i = 1, 2 and 3, to understand Kessler's Matlab code line by line.

b.) Reproduce Kessler's results.

a.) Line by Line explanation of Kessler's Code
Below is Kessler's Matlab code with the addition of our comments.

b.) Reproduction of Kessler's Results
The code above was used to generate the following two tables which perfectly matches the results in Kessler's paper linked in the references section:

Objective
Compute circumference of : 1).bifolium; 2). ellipse using the following methods: a). composite trap b). Romberg c). clenshaw-curtis quad d). Fast CC winckel

1.bifolium
$$r(t)=2\sin (t)\cos {{(t)}^{2}},t\in [0,\pi ]$$ The circumference of bifolium can be computed by: $$C=\int\limits_{t=0}^{\pi }{\sqrt{{{r}^{2}}+{{(\frac{dr}{dt})}^{2}}}}dt$$

Four methods to compute the integral: a).Composite trap Use error estimate for Composite Trapezoidal Rule to find n, such that the error is to the order of $$\displaystyle O(10^{-8}) $$.

Use successive numerical Integration results as a stopping criterion, i.e,

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

$$ $$
 * $$\displaystyle
 * I_{n+2} - I_{n}| < 10^{-8}
 * I_{n+2} - I_{n}| < 10^{-8}
 * $$\displaystyle
 * }.
 * }.

 Matlab Code:  b).Romberg The first column of Romberg table is calculated by composite trap, and then use the following equation to continue. $$\displaystyle T_{k}(2n) = \frac{2^{2k}T_{k-1}(2n) - T_{k-1}(n)}{2^{2k}-1}. $$  Matlab Code:  c).Clenshaw-Curtis Using Clenshaw-Curtis matlab code to compute the result.  Matlab Code:  d).Fast CC Using fast Clenshaw-Curtis matlab code to compute the result.  Matlab Code:  The results from the different methods are compared here:

2.ellipse-A

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

$$\displaystyle r(\theta) = \frac{a(1-e^2)}{1-e \, cos(\theta)} $$ $$
 * $$\displaystyle
 * }.
 * }.

Here, $$\displaystyle a=10,b=3 $$ and $$\displaystyle e=\sqrt{1-{{\left( \frac{b}{a} \right)}^{2}}} $$

Arc Length for an elliptical curve with $$\displaystyle \theta$$ varying from 0 to $$\displaystyle 2\pi $$ is given by,

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

$$\displaystyle

I=\int_{0}^{2\pi }{\sqrt{{{r}^{2}}+{{\left( \frac{dr}{d\theta } \right)}^{2}}}}d\theta

$$ $$ Four methods to compute the integral: a).Composite trap Use error estimate for Composite Trapezoidal Rule to find n, such that the error is to the order of $$\displaystyle O(10^{-10}) $$.
 * $$\displaystyle
 * }.
 * }.

Use successive numerical Integration results as a stopping criterion, i.e,

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

$$ $$
 * $$\displaystyle
 * I_{n+2} - I_{n}| < 10^{-10}
 * I_{n+2} - I_{n}| < 10^{-10}
 * $$\displaystyle
 * }.
 * }.

b).Romberg The first column of Romberg table is calculated by composite trap, and then use the following equation to continue. $$\displaystyle T_{k}(2n) = \frac{2^{2k}T_{k-1}(2n) - T_{k-1}(n)}{2^{2k}-1}. $$

c).Clenshaw-Curtis Using Clenshaw-Curtis matlab code to compute the result. d).Fast CC Using fast Clenshaw-Curtis matlab code to compute the result.

The results from the different methods are compared here:

3.ellipse-B
The elliptic integral of the second kind is given as:

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

$$\displaystyle I=4aE(e)=4a\int_{0}^{\pi /2}{\sqrt{1-{{e}^{2}}{{\sin }^{2}}\theta }}\ d\theta $$ $$ Four methods to compute the integral: a).Composite trap Use error estimate for Composite Trapezoidal Rule to find n, such that the error is to the order of $$\displaystyle O(10^{-10}) $$.
 * $$\displaystyle
 * }.
 * }.

Use successive numerical Integration results as a stopping criterion, i.e,

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

$$ $$
 * $$\displaystyle
 * I_{n+2} - I_{n}| < 10^{-10}
 * I_{n+2} - I_{n}| < 10^{-10}
 * $$\displaystyle
 * }.
 * }.

b).Romberg The first column of Romberg table is calculated by composite trap, and then use the following equation to continue. $$\displaystyle T_{k}(2n) = \frac{2^{2k}T_{k-1}(2n) - T_{k-1}(n)}{2^{2k}-1}. $$

c).Clenshaw-Curtis Using Clenshaw-Curtis matlab code to compute the result. d).Fast CC Using fast Clenshaw-Curtis matlab code to compute the result.

The results from the different methods are compared here:

Observations
1.For all the problem, Romberg method always need less time than Composite trap method. 2.The second kind of ellips (ellip-B) need less time than ellip-A. 3.The fast CC method always performs better than CC method. 4.For the problem of ellip-B, the computation time is similar for the four methods. This problem was solved by shengfeng yang =Problem 6.8: Derive the Arc Length Equation Using Triangle OAB and the Law of Cosines=

Problem Statement
Given on Lecture slide 33-2

Given
Triangle OAB

Objective
Prove

 * $$ \overline{PQ}= \int_{\theta_{P}}^{\theta_{Q}} \left[ r^2+ \left( \frac{dr}{d\theta}\right)^2\right]^{\frac{1}{2}} d\theta$$

Solution
Using the Law of Cosines We can express the length as.

$$ \overline{AB}^{2}=\overline{OA}^{2}+\overline{OB}^{2}-2(\overline{OA})( \overline{OB})cos(d\theta) $$

Squaring both sides

$$ \overline{AB}=\sqrt{\overline{OA}^{2}+\overline{OB}^{2}-2(\overline{OA} )(\overline{OB})cos(d\theta)} $$

If we now make the following substitutions

$$ \overline{AB}=dl $$

$$ \overline{OA}=r(\theta) $$

$$ \overline{OB}=r(\theta+d\theta) $$

$$ \overline{AB}=dl $$

We obtain

$$ dl=\sqrt{(r(\theta))^{2}+(r(\theta+d\theta))^{2}-2(r(\theta))(r(\theta+d\theta))cos(d\theta)} $$

This equation can be better expressed if we regroup the first term into a quadratic. This will allow us to use some simplifying assumptions later on.

$$ dl=\sqrt{(r(\theta)+r(\theta+d\theta))^{2}+2(r(\theta))(r(\theta+d\theta))(1-cos(d\theta))} $$

We can now use the following given information to simplify the first component of this equation.

$$ r(\theta+d\theta)-r(\theta) \approx dr $$

$$ dl=\sqrt{(dr)^{2}+2(r(\theta))(r(\theta+d\theta))(1-cos(d\theta))} $$

We may now use the cosine double angle formula to simplify the second component.

 * $$ cos(2\theta)=1-2sin^{2}(\theta) $$

When applied to the second component we get

$$ (1-cos(d\theta))=1-(1-2sin^{2}(\frac{\theta}{2})=2sin^{2}(\frac{\theta}{2}) $$

If we continue to assume that the angle is small, this can be simplified to

$$ (1-cos(d\theta))\approx2(\frac{d\theta}{2})^{2} =\frac{d\theta^{2}}{2} $$

Through substitution, we now have

$$ dl \approx \sqrt{(dr)^{2}+2(r(\theta))(r(\theta+d\theta))\frac{d\theta^{2}}{2}} $$

$$ dl \approx \sqrt{dr^{2}+r(\theta)r(\theta+d\theta)d\theta^{2}} $$

$$ dl \approx d\theta \sqrt{(\frac{dr}{d\theta})^{2}+r(\theta)r(\theta+d\theta)} $$

In keeping with our small angle assumption we can say that the term $$ r(\theta)r(\theta+d\theta) $$ will be dominated by the $$r(\theta)$$ components, leading us to the following small angle approximation.

$$ r(\theta)r(\theta+d\theta) $$

Substituting this back into our equation, we obtain

$$ dl \approx d\theta \sqrt{(\frac{dr}{d\theta})^{2}+r(\theta)^{2}} $$

If we now integrate this to find the arc length $$\overline{PQ} $$ we come up with the following definite integral.

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

 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ \overline{PQ} = \int_{\theta_{P}}^{\theta_{Q}} \sqrt{r^{2}+(\frac{dr}{d\theta})^{2}} d\theta$$
 * }

Solved by Brendan Mahon

Given

 * {| style="width:70%" border="0" align="center"

\begin{Bmatrix} c_0 \\ c _1 \\ c _2 \\ c_3 \end{Bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -3 & -2 & 3 & -1 \\ 2 & 1 & -2 & 1 \end{bmatrix} \begin{Bmatrix} Z_i \\ Z_i^' \\ Z_{i+1} \\ Z_{i+1}^' \end{Bmatrix} =\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -3 & -2 & 3 & -1 \\ 2 & 1 & -2 & 1 \end{bmatrix} \begin{Bmatrix} \overline{d}_1 \\ \overline{d}_2 \\ \overline{d}_3 \\ \overline{d}_4 \\ \end{Bmatrix} $$     (9.1)
 * $$\displaystyle
 * $$\displaystyle
 *
 * }
 * }

Where:

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

\begin{Bmatrix} Z_i \\ Z_i^' \\ Z_{i+1} \\ Z_{i+1}^' \end{Bmatrix} = \begin{Bmatrix} \overline{d}_1 \\ \overline{d}_2 \\ \overline{d}_3 \\ \overline{d}_4 \\ \end{Bmatrix} = \begin{Bmatrix} d_1 \\ hd_2 \\ d_3 \\ hd_4 \end{Bmatrix} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

And equation (1) on lecture slide [[media:nm1.s11.mtg37.djvu|37-1]]:

 * {| style="width:70%" border="0" align="center"

$$     (9.2)
 * $$\displaystyle Z(s) = \sum_{i=0}^3 c_i s^i = \sum_{i=1}^4 \overline{N}_i(s)d_i
 * $$\displaystyle Z(s) = \sum_{i=0}^3 c_i s^i = \sum_{i=1}^4 \overline{N}_i(s)d_i
 *
 * }
 * }

Objective
a.) Identify the Basis Functions Ni(s) for i = 1,2,3,4

b.) Plot the Basis Functions Ni(s) for i = 1,2,3,4

a.) Basis Functions Ni(s) for i = 1,2,3,4
Expanding the Matrix equation given in (9.1) we get:

 * {| style="width:70%" border="0" align="center"

$$     (9.3)
 * $$\displaystyle c_0 = Z_i
 * $$\displaystyle c_0 = Z_i
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.4)
 * $$\displaystyle c_1 = Z_i^'
 * $$\displaystyle c_1 = Z_i^'
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.5)
 * $$\displaystyle c_2 = -3Z_i - 2Z_i^' + 3Z_{i+1} - Z_{i+1}^'
 * $$\displaystyle c_2 = -3Z_i - 2Z_i^' + 3Z_{i+1} - Z_{i+1}^'
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.6)
 * $$\displaystyle c_3 = 2Z_i + Z_i^' - 2Z_{i+1} + Z_{i+1}^'
 * $$\displaystyle c_3 = 2Z_i + Z_i^' - 2Z_{i+1} + Z_{i+1}^'
 *
 * }
 * }

By expanding equation (9.2) and then plugging the above values for Ci into the expanded equation we obtain the following:

 * {| style="width:70%" border="0" align="center"

$$     (9.7)
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4 = C_0s^0 + C_1s^1 + C_2s^2 + C_3s^3
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4 = C_0s^0 + C_1s^1 + C_2s^2 + C_3s^3
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.8)
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

= Z_i + (Z_i^')s^1 + (-3Z_i - 2Z_i^' + 3Z_{i+1} - Z_{i+1}^')s^2 + (2Z_i + Z_i^' - 2Z_{i+1} + Z_{i+1}^')s^3 $$
 * $$\displaystyle
 * $$\displaystyle

 * }
 * }

Rearranging the right side of equation (9.8) we get the following:

 * {| style="width:70%" border="0" align="center"

$$     (9.9)
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

= (1 -3s^2 + 2s^3)Z_i + (s - 2s^2 + s^3)z_i^' + (3s^2 - 2s^3)Z_{i+1} + (-s^2 + s^3)Z_{i+1}^' $$
 * $$\displaystyle
 * $$\displaystyle

 * }
 * }

By comparing the left and right sides of equation (9.9) we obtain the following basis functions:

 * {| style="width:70%" border="0" align="center"

$$     (9.10)
 * $$\displaystyle N_1(s) = 1 -3s^2 + 2s^3
 * $$\displaystyle N_1(s) = 1 -3s^2 + 2s^3
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.11)
 * $$\displaystyle N_2(s) = s - 2s^2 + s^3
 * $$\displaystyle N_2(s) = s - 2s^2 + s^3
 * <p style="text-align:right">
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.12)
 * $$\displaystyle N_3(s) = 3s^2 - 2s^3
 * $$\displaystyle N_3(s) = 3s^2 - 2s^3
 * <p style="text-align:right">
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.13)
 * $$\displaystyle N_4(s) = -s^2 + s^3
 * $$\displaystyle N_4(s) = -s^2 + s^3
 * <p style="text-align:right">
 * }
 * }

b.) Plots of the Basis Functions Ni(s) for i = 1,2,3,4
The following Matlab code was used to create the figure below.

>> s = linspace(0,1,100); >> N1 = 1 - 3.*s.^2 + 2.*s.^3; >> N2 = s - 2.*s.^2 + s.^3; >> N3 = 3.*s.^2 - 2.*s.^3; >> N4 = -s.^2 + s.^3; >> plot(s,N1,s,N2,s,N3,s,N4)

title: User:Egm6341.s11.team4/HW6


text:
User:Valekvespa

Valekvespa is a volunteer in sharing knowledge to other people. He knows that this is what Wikiversity purpose.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

title: User:Valekvespa


title: File:Fe1.s11.team1.hw6.fig6.1.1.jpg


text:
User:Egm6341.s11.team4.fields/HW6

HW6.6

HW6.9

title: User:Egm6341.s11.team4.fields/HW6


text: File:S.11.Team6.tupsakhare-Figure- 6.3.6-Pre-processed line between P1 and P2(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.7-Pre processing Points P1 to P5 with 4 lines(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.8-figure Points P1 to P5 with 4 lines(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.9-pnt and line names display(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.10-(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.11-(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.12-1-(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.12-2-(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.13(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.14(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.15(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.16(1).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.1-Save as(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.2-(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.3-(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.4(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.6(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.8-(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.9-(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.10-(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.11-(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.12-(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.13(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.14-1(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.14-2(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.15(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.16(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.17(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.18(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.19(2).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.1-(3).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.2-(3).png


title: File:Fe1.s11.team1.hw6.fig6.1.2.jpg


text: File:S.11.Team6.tupsakhare-Figure- 6.3.3-(3).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.4-(3).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.5-(3).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.1-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.2-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.3.-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.4-(4).png


title: File:Fe1.s11.team1.hw6.fig6.1.3.jpg


text: File:S.11.Team6.tupsakhare-Figure- 6.3.5-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.6-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.7-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.8-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.9-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.10-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.11-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.12-(4).png


text: File:S.11.Team6.tupsakhare-Figure- 6.3.13-(4).png


title: File:Fe1.s11.team1.hw6.fig6.1.4.jpg


text:
User:Norton~enwikiversity/Wikistution/Amendment XXI (Alcohol Consumption Age)

Amendment XXI (1933)

Section 1

The eighteenth article of amendment to the Constitution of the United States is hereby repealed.

Section 2

The transportation or importation into any State, Territory, or Possession of the United States for delivery or use therein of intoxicating liquors, in violation of the laws thereof, is hereby prohibited.

Section 3

This article shall be inoperative unless it shall have been ratified as an amendment to the Constitution by conventions in the several States, as provided in the Constitution, within seven years from the date of the submission hereof to the States by the Congress.

Section 4

The legal drinking age of the United States is 18 years of age. The citizen wishing to drink will be required to obtain authorization in the form of an alcohol consumption license. This license can be obtained by passing an alcohol awareness exam at a state run facility and completing a probationary period. Additionally, the alcohol consumption license can be suspended or revoked upon receiving alcohol related offenses, in accordance with the severity of the crime.


text:
User:Keith1

Mainly interested in pursuing a relevant, useful Wikistution.

title: User:Keith1


title: File:Fe1.s11.team3.HW6.9.fig1.jpg


title: File:Fe1.s11.team3.HW6.9.fig2.jpg


title: File:Fe1.s11.team3.HW6.9.fig3.jpg


title: File:Fe1.s11.team3.HW6.9.fig4.jpg


title: File:Fe1.s11.team3.HW6.9.fig5.jpg


title: File:Fe1.s11.team3.HW6.9.fig6.jpg


title: File:Fe1.s11.team3.HW6.9.fig7.jpg


title: File:Fe1.s11.team3.HW6.9.fig8.jpg


title: File:Fe1.s11.team3.HW6.9.fig9.jpg


title: File:Nm1.s11.team1.HW5.fig5.1.1.png


title: File:Nm1.s11.team1.HW5.fig5.1.2.png


text:
Quiz on Height of the Ottoman Empire, Part II

{Why did Ottoman army enter Morocco in period of Murad III? -For Jihad -For trade +For domination -For Spanish help
 * type=""}

{What happened in 1583? +Battle of Torches -Battle of Lepanto -Battle of Kanizsa -Treaty of Ferhat Pasha
 * type=""}

{What happened in North (End of 16th Century)? -Polish army captured Moskow -Persian army entered Moskow -Russian army destroyed Crimean Khanete +Russian army captured some Tatarian Khanetes.
 * type=""}

{When did Long War begin? -1581 -1587 +1593 -1599
 * type=""}

{What happened in Kanizsa? -Austrian - Hungarian army destroyed Ottoman army +Ottoman army destroyed Austrian - Hungarian army with cheating -Sultanete army reached Kanizsa and Austrian - Hungarian army was destroyed -Crimean soldiers reached Kanizsa and they were destroyed.
 * type=""}

{What happened after Battle of Kanizsa? -Austrian - Hungarian army entered Balkans -Austrian - Hungarian army captured all Hungary and Croatia -Ottoman army was destroyed in Slovenia +Ottoman army entered Austrian Area
 * type=""}

{What happened in Anatolia? -Persian Revolt began +Jelali Revolt began -Armenian Revolt began -Jelali Revolt end
 * type=""}

{Who killed Shehzade Mahmud? -Selim II -Kuyucu Murad Pasha -Lala Mehmed Pasha +Mehmed III
 * type=""}

{When did Zsitvatorok Peace was sign? -1604 -1605 +1606 -1607
 * type=""}

{(Research question) What did Ottomans in science? -They worked about Alchemy +They worked about Literature and architecture -They worked about chemistry -They worked about archeology
 * type=""}

{(Research question) Which is famous cartographer? -Lala Mehmed Pasha -Kuyucu Murad Pasha -Hasan Reşis +Piri Reis
 * type=""}


text:
User:Rgulami/Reflective Practical Analysis

Reflective Practical Analysis


text:
User:Eml5526.s11.team4.yang/hw6

=HW 6.1 problem: Solving G1DM1.0/D1 using QLEBF=

Problem Statement
Similar to HW5.{1,3,7},but using QLEBF with uniform discretization(equidistant element nodes),nel=2,4,6,8... 1). For nel=2, comp $$\tilde{K}=\sum\limits_{e=1}^{2}$$ ,with by (1)p30.5,display $${{\tilde{K}}^{e}}$$, e=1,2 2). comp $${{k}^{e}}$$ , $${{L}^{e}}$$, for e=1,2 3).comp $${{{\tilde{K}}}^{e}}={{L}^{e}}^{T}{{k}^{e}}{{L}^{e}}$$, for e=1,2. compare to 1). 4).plot all QLEBF for nel=3. 5).plot $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$, $$\displaystyle \left[ u_{nel}^{h}(0.5)-u(0.5) \right]$$ vs number of element

solution
1.For nel=2, comp $$\tilde{K}=\sum\limits_{e=1}^{2}$$ ,with by (1)p30.5,display $${{\tilde{K}}^{e}}$$, e=1,2 e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix} { 0} & {0} & {0} & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Then,sum on e=1,2:

$$\tilde{K}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {32.666666666666664}  &{-20.666666666666668} & {2.833333333333334}\\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

2.comp $${{k}^{e}}$$ , $${{L}^{e}}$$, for e=1,2 e=1:

$${{k}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668}  \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  \\ \end{matrix} \right]$$

$${{L}^{1}}=\left[ \begin{matrix} 1 & 0 & 0 & 0 & 0 \\   0 & 1 & 0 & 0 & 0  \\   0 & 0 & 1 & 0 & 0  \\ \end{matrix} \right]$$

e=2:

$${{k}^{2}}=\left[ \begin{matrix} {17.833333333333332} & {-20.666666666666668} & {2.833333333333334} \\   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

$${{L}^{2}}=\left[ \begin{matrix} 0 & 0 & 1 & 0 & 0 \\   0 & 0 & 0 & 1 & 0  \\   0 & 0 & 0 & 0 & 1  \\ \end{matrix} \right]$$

3.comp $${{{\tilde{K}}}^{e}}={{L}^{e}}^{T}{{k}^{e}}{{L}^{e}}$$, for e=1,2. compare to 1). e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix}  {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix}   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Compare with result on 1), we can see that the result is the same by using these two methods. 4.plot all QLEBF for nel=3 5.plot $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$,when number of element is 2,3,4,6,respectively. plot $$\displaystyle \left[ u_{nel}^{h}(0.5)-u(0.5) \right]$$ vs number of element

The error at x=0.5 can be listed in the following table: From the table, it can be found that the error will be lower than 10^(-6) when the number of elements is greater than n=14.

This problem was solved by shengfeng yang

=HW 6.2 problem: Solving G1DM1.0/D1b using QLEBF=

Problem Statement
Similar to HW5.{2,4,8},but using QLEBF with uniform discretization(equidistant element nodes),nel=2,4,6,8... 1). For nel=2, comp $$\tilde{K}=\sum\limits_{e=1}^{2}$$ ,with by (1)p30.5,display $${{\tilde{K}}^{e}}$$, e=1,2 2). comp $${{k}^{e}}$$ , $${{L}^{e}}$$, for e=1,2 3).comp $${{{\tilde{K}}}^{e}}={{L}^{e}}^{T}{{k}^{e}}{{L}^{e}}$$, for e=1,2. compare to 1). 4).plot all QLEBF for nel=3. 5).plot $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$, $$\displaystyle \left[ u_{nel}^{h}(0.5)-u(0.5) \right]$$ vs number of element

solution
1.For nel=2, comp $$\tilde{K}=\sum\limits_{e=1}^{2}$$ ,with by (1)p30.5,display $${{\tilde{K}}^{e}}$$, e=1,2 e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix} { 0} & {0} & {0} & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Then,sum on e=1,2:

$$\tilde{K}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {32.666666666666664}  &{-20.666666666666668} & {2.833333333333334}\\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

2.comp $${{k}^{e}}$$ , $${{L}^{e}}$$, for e=1,2 e=1:

$${{k}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668}  \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  \\ \end{matrix} \right]$$

$${{L}^{1}}=\left[ \begin{matrix} 1 & 0 & 0 & 0 & 0 \\   0 & 1 & 0 & 0 & 0  \\   0 & 0 & 1 & 0 & 0  \\ \end{matrix} \right]$$

e=2:

$${{k}^{2}}=\left[ \begin{matrix} {17.833333333333332} & {-20.666666666666668} & {2.833333333333334} \\   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

$${{L}^{2}}=\left[ \begin{matrix} 0 & 0 & 1 & 0 & 0 \\   0 & 0 & 0 & 1 & 0  \\   0 & 0 & 0 & 0 & 1  \\ \end{matrix} \right]$$

3.comp $${{{\tilde{K}}}^{e}}={{L}^{e}}^{T}{{k}^{e}}{{L}^{e}}$$, for e=1,2. compare to 1). e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix}  {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix}   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Compare with result on 1), we can see that the result is the same by using these two methods. 4.plot all QLEBF for nel=3 5.plot $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$,when number of element is 2,3,4,6,10,respectively.

plot $$\displaystyle \left[ u_{nel}^{h}(0.5)-u(0.5) \right]$$ vs number of element

The error at x=0.5 can be listed in the following table: From the table, it can be found that the error will be lower than 10^(-6) when the number of elements is greater than n=8.

This problem was solved by shengfeng yang

title: User:Eml5526.s11.team4.yang/hw6


title: File:Hw6.1-me=2.jpg


title: File:Hw6.1-me=3.jpg


title: File:Hw6.1-me=4.jpg


title: File:Hw6.1-me=6.jpg


title: File:Hw6.1-me=10.jpg


title: File:Hw6.1-err.jpg


title: File:Hw6.1-me=3-basis.jpg


title: File:Hw6.2-me=2.jpg


title: File:Hw6.2-me=3.jpg


title: File:Hw6.2-me=4.jpg


title: File:Hw6.2-me=6.jpg


text:
Talk:Cosmological simulations

Links
The following link was removed from the page, as it is broken: Gadget-2 initial conditions for the Antennae model --mikeu talk 18:08, 3 April 2011 (UTC)
 * Link to archive.org added: Data Files and Public Movies. --mikeu talk 18:47, 21 January 2018 (UTC)

How to enroll
Students at WSU-Lake will get course credit for contributing to this resource. Those who are not enrolled but wish to edit this resource are welcome to do so. If you are not a WSU-Lake student and  wish to see your user name in the spreadsheet shown below, leave a message at below.

To do

 * 1) Find a place for non WSU-Lake students to sign up.
 * 2) Create a table like the one at Wright State University Lake Campus/2016-1/Phy1060
 * 3) Google search to verify that nobody important has taken the name EduV.

title: Talk:Cosmological simulations


title: File:Hw6.2-err.jpg


text:
User:Eml5526.s11.team01/Hwk6

 HOMEWORK 6

=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 $$
 * <p style="text-align:right;">$$\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:

=Problem 2:Weak Form Using Quadratic Lagrangian elemental Basis Function With Uniform Discretization.=

Given
General 1-D Model 1 Data set 1. The Partial Differential equation for this data set is,

Find
1)For $$\displaystyle n_{el}=2$$ Compute $$\underline{\tilde{K}}$$ Using Equation 1 On page 30-5. Display $$\underline{\tilde{K}},e=1,2$$

2)Compute $$K^{e}$$, $$L^{e}$$ for $$e=1,2$$

3)Compute $$\underline{\tilde{K}}^{e}=\underline{\tilde{L}}^{e}\underline{K}^{e}\underline{L}^{e}, for e=1,2$$ & Compare to 1

4)Plot all QLEBF for $$n_{el}=3$$.

5)Plot $$\displaystyle \tilde{u_{n}}^{h}$$ vs $$\displaystyle u$$ and [$$\tilde{u_{n}}^{h}\left ( 0.5\right )-u\left (0.5\right )$$] vs $$\tilde{n}$$

Solution
Lets Obtain the Weak form first of all as follow

Performing Integration by parts technique we get

Solving this we get,

As we do not know $$ \frac{du}{dx}\left ( 1\right )$$ we choose $$w\left ( 1 \right )=0$$ and moreover $$u^{h}\left(1 \right )=4$$. To use quadratic Lagrangian elemental basis function we divide the dominium in elements, on each element we apply the weak form by using quadratic Lagrangian Polynomial. For that we need 3 nodes per element. Each element with global nodes i,i+1,i+2 will have 3 local nodes 1,2,3. Lets call $$x_{1},x_{2},x_{3}$$the co-ordinates at node 1,2,3 respectively.

To use natural co-ordinate $$\epsilon$$ instead of Co-ordinate $$x$$ we use the equation 4 on page 30-3 for that we know the natural co-ordinates of the nodes are -1,0,1 & local nodes 1,2,3 respectively. Using Lagrange Polynomial to map between two Co-ordinate system we have:

substituting the above values in the $$\displaystyle (6.2.7)$$ we get,

In this way we have found a shape function, and the later equation can be written as follow,

In the same way we can interpolate the function $$u$$ & $$w$$, therefore,

The Derivative of function $$u$$ with respect to $$\epsilon$$ is given by, In the same way, If we place the Node 2 at the midpoint of the element we can write the relation between the two co-ordinates as Taking the derivative on both sides we get, The derivative of $$u$$ & $$w$$ with respect to $$x$$ are calculated by using chain rule which is as below, Replacing all the terms in weak form as a function of $$\epsilon$$ we obtain This can be further simplified as

From the later Equation we can get the Stiffness Matrix of one element as follows,

and the force vector of the element will be,

From the later Equation we can get the Stiffness Matrix of one element as follows,

There will be additional following factor,

From the later Equation we can get the Stiffness Matrix of one element as follows,

Now only count for degree of freedom at x=1.

1) For $$\displaystyle n_{el}=2$$ Compute $$\underline{\tilde{K}}$$ Using Equation 1 On page 30-5. Display $$\underline{\tilde{K}},e=1,2$$.
Now according to the Lagrange Interpolation we have following equation,

Also,

Replacing in the weak form above we get, The last term is taken only for the last element. Further the stiffness matrix is

Performing the integration for the code shown below for one element we obtain, Hence the global stiffness matrix will be,

2) Compute $$K^{e}$$, $$L^{e}$$ for $$e=1,2$$.
$$K^{e}$$ is already calculated above. For $$L^{e}$$ As the element has 3 degrees of freedom & the total number of degrees of freedom is 5 The L matrix will be of size $$3\times5$$. Therefore,

3) Compute $$\underline{\tilde{K}}^{e}=\underline{\tilde{L}}^{e}\underline{K}^{e}\underline{L}^{e}, for e=1,2$$ & Compare to 1.
Now, Solving the above matrix we get,

For element 2 we have,

Solving the above matrix we get,

Hence Combining equations $$\displaystyle (6.2.17)$$ & $$\displaystyle (6.2.18)$$ we get global stiffness matrix as follows,

Eventually Comparing $$\displaystyle (6.2.15)$$ with $$\displaystyle (6.2.19)$$ we found that the Stiffness matrices are identical.

4) Plot all QLEBF for $$n_{el}=3$$.
The Quadratic Lagrange Elemental basis Function is plotted using the following MATLAB code for $$\displaystyle n_{el}=3$$.

Appendix
 Matlab Code:

The plot is as follow.

5) Plot $$\displaystyle \tilde{u_{n}}^{h}$$ vs $$\displaystyle u$$ and [$$\tilde{u_{n}}^{h}\left ( 0.5\right )-u\left (0.5\right )$$] vs $$\tilde{n}$$.
These Two plots are plotted as below Using the mentioned MATLAB code.

Appendix
 Matlab Code:

===Comment: From the above plot we can observe that the solution converges at nel=16. The Convergence of $$10^{-6}$$ was difficult to achieve therefore the convergence criteria was modified to $$E=\left(u_{0.5}^{i}- u_{0.5}^{i-1}\right )$$.===

=Problem 3: Calculix=

Given
Given the tutorial found in this tutorial.

Find
Reproduce the actions of the author in Calculix

Solution
The tutorial was followed and reproduced sucessfully. The actions and results were summerized in the following

=Problem 4: FB pp.148-149 Problems 6.1,6.2,6.3=

Given
6.1)

Given vector $$ q_x = -y^2 ~,~ q_y = -2xy $$ on domain shown in Figure 6.2.

6.2)

Given a vector field $$ q_x = 3x^2 y +y^3 ~,~ q_y = 3x + y^3 $$ on the domain shown in Figure 6.9. The curved boundary of the domain is a parabola.

6.3) Given $$ \oint_{T} \mathbf{n} d\Gamma = 0 $$

Find
6.1)

Verify the divergence theorem.

6.2)

Verify the divergence theorem.

6.3)

Prove using divergence theorem

FB 6.1)

We have to show that:

Where $$ \Omega $$ is the domain, $$ \Gamma $$ is the boundary, n is the normal vector and the notation $$ q_{i,i} $$  means $$ \frac{ \partial q_i}{ \partial x_i} $$

The integration on the boundary is

FB 6.2)
To calculate the integral on the surface we see that part of the surface is the parabola and the other part is in the x axis; therefore the integral will be:

The first integral is on the parabola and the second part is on the x axis.

On the curved surface the, normal vector on a segment of surface ds is shown in the figure. From that figure we observe that:

Also, it is clear that on the x axis we have y=0. Replacing in the integral we obtain:

Taking the derivative on the equation of the parabola we obtain:

Replacing on the integral we have:

Changing the limits of integration in the first two terms and organizing we obtain:

The integral on the volume will be:

Performing first the integration with respect to y we have:

FB 6.3)
The theorem says that for any tensor $$ T_{i...j} $$ (of any order) inside a volume $$ \Omega $$ surrounded by a surface $$ \Gamma $$ with a normal vector n

In this case the tensor T is the constant 1; therefore, changing the integral of surface to an integral of volume we have

But as 1 is a constant, its derivatives with respect to any coordinate are zero and we have

=Problem 6.5: Solve G2DM10/D1 using 2DLIBF=

Given: 2DLIBF from p. 29-2
PDE

$$ K = I $$ (identity matrix)

$$ f = 0, \frac{\partial u}{\partial t} = 0 $$

$$ \frac{\partial}{\partial x_i} \left(k_{ij} \frac{\partial u}{\partial x_j} \right) = 0 $$

Essential B.C.

$$ g = 2 on  \partial \Omega $$

Natural B.C.

none

Find
Solve G2DM10/D1 for $$ u^h $$ until accuracy $$ 10^{-6} $$ at center $$ (x,y) = (0,0) $$

Solution
Using the weak form to be

and

2DLIBF

where

$$ I = i + \left( j-1 \right) m $$

$$ n = m = 2, 4, 6, 8... $$

Thus

Therefore by substituting equation 6.5.3 into equation 6.5.1 and integrating, the stiffness matrix can be solved for. This is done by using Matlab.

The solution yields

=Problem 6=

Solution
=Problem 7=

Solution
= Problem 6.8 =

Verify the Table of {(wi, xi), i=1,2,3,4,5} against the NIST Handbook and that in the textbook
The Table of $$\left\{ \left( {{w}_{i}},{{x}_{i}} \right),\ i=1,2,3,4,5 \right\}$$ is shown as follows:

Complare the positions of points and the weight with respect to each of them in the Table above with NIST Handbook and Fish and Belytschko's book on page 89.

Solution
Calculation of the value of each term in the above table gives the following result:

The nodes and weights for the 5-point Gauss–Legendre formula in the NIST Handbook is the following:

The position of Gauss points and corresponding weights in Fish and Belytschko's book on page 89 can be seen as follows:

Conclusion:

It can be seen from the above several tables that the table of $$\left\{ \left( {{w}_{i}},{{x}_{i}} \right),\ i=1,2,3,4,5 \right\}$$ given in Lecture 36 gives us quite precise results of the position of Gauss points and corresponding weights compared to the other two authoritative sources.

Use Gauss quadrature to obtain exact values for the given integrals and verify by analytical integration.
(c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB code.

Solution
(a)

Select the number of Gauss point, n:

$$2n-1\ge 2$$, select n=2.

Then,

By analytical integration:

So, the result performed by analytical integration is identical to that obtained by Gauss quadrature.

(b)

Select the number of Gauss point, n:

$$2n-1\ge 4$$, select n=3.

By analytical integration:

So, the result performed by analytical integration is identical to that obtained by Gauss quadrature.

(c)

Comments:

The results for (a) and (b) using MATLAB code are 25.3333 and 1.7333  respectively, which are the same with those calculated by analytical integration.

Use three-point Gauss quadrature to evaluate the given integrals and compare to the analytical integration.
(c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB code.

Solution
(a)

n=3, then we have the following according to the Table of {(wi, xi)}:

So,

Integrate using WolframAlpha, we get the same result.

(b)

n=3, then we have the following according to the Table of {(wi, xi)}:

So,

Integrate using WolframAlpha, we get the result of I=1, which is different from what we get using three-point Gauss quadrature.

(c)

Comments:

The result for (a) using MATLAB code is 0, which is the same with that calculated by analytical integration.

And the result for (b) using MATLAB code is 1.53, which is different from that calculated by analytical integration.

Employ one-point and three-point Gauss quadrature and evaluate the accuracy of the given integral compared to the result by two-point quadrature.
The integral $$\int\limits_{-1}^{1}{\left( 3{{\xi }^{3}}+2 \right)d\xi }$$ 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.

Check your calculations against MATLAB code.

Solution
a.

n=1, then we have the following according to the Table of {(wi, xi)}:

So,

b.

n=3, then we have the following according to the Table of {(wi, xi)}:

So,

Comments:

The results for (a) and (b) using MATLAB code are all equal to 4 , which are the same with those calculated by analytical integration.

Conclusion:

It can be seen from the above Gauss quadrature calculation and the verification using MATLAB code that the accuracy is not affected when the number of Gauss points is either increased or decreased in this problem.

Matlab Code
=Problem 9=

Solve the 2D heat conduction problem with 16 elements using the finite element code.
Consider the heat conduction problem depicted in Figure 6.9.1. The coordinates are given in meters. The conductivity is isotropic, with $$\mathbf{D}=k\left[ \begin{matrix} 1 & 0 \\   0 & 1  \\ \end{matrix} \right]$$, and $$k=5\text{W}{}^{\text{o}}{{\text{C}}^{\text{-1}}}$$. The temperature $$T=0$$ is prescribed along edges AB and AD. The heat fluxes $$\bar{q}=0$$ and $$\bar{q}=20\text{W}{{\text{m}}^{\text{-1}}}$$ are prescribed on edges BC and CD, respectively. A constant heat source $$s=6\text{W}{{\text{m}}^{\text{-2}}}$$ is applied over the plate. Consider this problem modeled with 16 quadrilateral finite elements as shown in Figure 6.9.2.

Solve this problem using the finite element code given in Section 12.5.

Solution

The MATLAB code are downloaded at FB Student Companion Site.

Show the temperature distribution for the coarsest (34-element) and the finest (502-elements) meshes.
In this example, we consider a manufactured solution of the form

defined over the domain of a square plate with a hole. For heat equation with isotropic conductivity and k=1, the corresponding source term is given by

The essential boundary conditions on $${{\Gamma }_{T}}$$ are

The natural boundary conditions on $${{\Gamma }_{q}}$$ are $$\left( \overline{q}=-k{{\mathbf{n}}^{T}}\nabla T \right)$$

Solution
Due to the hole in the midle of the domain, we find dificult to generate the mesh by hand, therefore we use a comercial program to solve the problem. First we start with a coarse mesh withwith 244 elements (Figure 6.9.5). The results are shown in Figure 6.5.6.

With 244 elements the results are show in the following figure.

Figure 6.9.7 shows the results across the plate at y=0.5.

Then the mesh was refine to 1116 elements (Figure 6.9.8), and again we show the results at y = 0.5 (Figure 6.9.9). When this results were compared with the analytical solution (Figure 6.9.11) an important difference was found. This diffence may be due to some error in our procedure, but also can be due to the size of the hole. As in the book the diemeter of the hole is not given we chose a diameter of 0.5 which apparently is lower than that used in the book.

The following figure shows the results with 1116 elements.

The analitical solution is shown below.

Determine the temperature and flux in the two materials of a chimney with 2x2, 4x4 and 8x8 quadrilateral elements for 1/8 of the problem domain.
Consider a chimney constructed of two isotropic materials: dense concrete ($$k=2.0W{}^{\circ }{{C}^{-1}}$$) and bricks ($$k=0.9W{}^{\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=10{}^{\circ }C$$. 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 2×2, 4×4 and 8×8 quadrilateral elements for 1/8 of the problem domain. A2×2 finite element mesh is shown in Figure below. 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
Figure 6.9.13 shows the initial mesh used to solve the problem. With this mesh, we obtained the temperature shown in Figure 6.9.14 and the heat flux shown in Figure 6.9.15. Then, we refine the mesh as shown in Figure 6.9.16 and obtained the temperature and heat flux shown in Figure 6.9.16 and 6.9.17 respectively. In the heat flux results we observed a high flux in the right border in the contact between the two materials. Finally, an extrafine mesh was used (Figure 6.9.19) to obtain the temperature shown in Figure 6.9.20 and the heat flux shown in Figure 6.9.21. It can be observed that with extrafine mesh the maximum values of heat flux obtained are higher than those obtained with coarse mesh and the regions of maximum heat flux are narrow. This behaviour was not observed on the results of temperature where the results are similar with the three mesh.

=Problem 10: Newtons Law of Colling=

Given
The partial Differential Equation is

1.Develop Weak form similar to (3)-(6) page 34-3.
The weak form will be as below,

Lets calculate the derivative of the following factor, Therefore, Replacing in the integral we obtain, Applying the Gauss Theorem on the first term we obtain, We can recognize that the Heat Flux is Given by, Eventually We have,

Organizing We get, Simplifying we get, For the Case when Conductivity if the Identity Matrix, we have following equations, Also We can Write as,

2) Develop Semidiscrete Equations (ODE's) similar to (3) page 23-3. Give Detailed Description for all the quantities.
Performing the Discretization on the Dominium we have,

Where N is the vector of the Interpolation (Lagrange Polynomial) and $$\displaystyle u_{e},w_{e}$$ are the Nodal Values of u & w Respectively. The Derivative With respect to $$\displaystyle x_{1}$$, $$\displaystyle x_{2}$$ are Replacing on the Semi Discrete Weak form we obtain, From the Last Equation we can Obtain the Conductivity Matrix, the Heat source vector and the Capacitance matrix. The Conductivity Matrix is Given by

Further the Capacitance Matrix is Given by,

And the Heat Source is

3.Solve G2DM2.0/D1 Using 2D LIBF (similar to HW 6.6) till $$10^{-6}$$ accuracy at center. Plot Solution in 3D w/Contour.
To solve G2DM2.0/D1 Following Routine was developed. The following figure shows the result for 16 points.

To solve the problem using LIBF we will devide the x-axis and Y-axis in 2 parts initially. In that way the Approximated solution by LIBF is given by:

We Have developed a MATLAB code to plot this routine. The MATLAB code is as follows

Appendix
 Matlab Code:

The Plot for the Routine Developed along with the intended 3-D contour is as below.

Comment: The Contour is plotted for $$\displaystyle n_{el}=4$$. The above MATLAB code is developed for larger values of $$n_{el}$$.
=Contributing Members=

= References =

title: User:Eml5526.s11.team01/Hwk6


text:
User:Egm6341.s11.team4.fields/HW6.6

=Problem 6.6: Kessler's Trapezoidal Rule Error Matlab Code=

Refer to lecture slide [[media:nm1.s11.mtg32.djvu|32-3]] for the problem statement.

Given
Kessler's Matlab code given in the following link:

Kessler's Trapezoidal Rule Error Matlab Code

Objective
a.) Use (P2i,P2i+1), i = 1, 2 and 3, to understand Kessler's Matlab code line by line.

b.) Reproduce Kessler's results.

a.) Line by Line explanation of Kessler's Code
Below is Kessler's Matlab code with the addition of our comments.

b.) Reproduction of Kessler's Results
The code above was used to generate the following two tables which perfectly matches the results in Kessler's paper linked in the references section:

title: User:Egm6341.s11.team4.fields/HW6.6


text:
User:Eml5526.s11.team6.kurth/HW6

=Problem 6.2: Solving Data Set G1DM1.0/D1 Using Quadratic Lagrangian Element Basis Functions and Uniform Discretization=

Boundary Conditions
The Essential Boundary condition as applied on the boundary $$\Gamma_g = [1]$$, is

The Natural Boundary condition as applied on the boundary $$\Gamma_h = [0]$$, is

Descritization
Here we want to use Uniform Descritization such that Element nodes are equidistant with number of elements as $$ nel = 2, 4, 6, 8, ..$$

Find
1.) For $$n_{el} = 2$$, compute $$\tilde{\underline{K}}=\sum_{e=1}^2\tilde{\underline{K}}^e$$, with $$\tilde{\underline{K}}^e$$ by [[media: Fe1.s11.mtg30.djvu| (1)p.30-5]], display ,$$\tilde{\underline{K}}^e,\;e=1,2.$$

2.) Comp. $$\underline{k}^e,\underline{L}^e$$, for $$\displaystyle e = 1,2.$$

3.) Comp. $$\tilde{\underline{K}}^e={\underline{L}^e}^T\underline{k}^e\underline{L}^e$$, for $$e = 1,2.$$ Compare to 1).

4.) Plot all QLEBF for $$n_{el}=3.$$

5.) Plot $$u_{\tilde{n}}^h$$ vs. $$\displaystyle u$$, $$[u_{\tilde{n}}^h(0.5)-u(0.5)]$$ vs. $$\tilde{n}$$

Solution
When using quadratic Lagrange element basis functions for this case we will employ 1-D quadratic elements, i.e. elements consisting of 3 nodes each. Figure 1 below shows the element structure for the case when nel=2.

Element Basis Functions and Their Derivatives
The quadratic Lagrange element basis functions are

Note that these functions take these values only within their corresponding element and are zero everywhere else.

Since the discretization used is uniform these basis functions may be simplified knowing the locations of the nodes. The simplified basis functions are

Then it follows that the derivatives of these functions are

Constructing the Global Stiffness Matrix
From equation [[media:Fe1.s11.mtg30.djvu|(1)p.30-5]]

Note: This matrix denotes the contribution due to element 1 to the global stiffness matrix. Since the basis functions $$b_4(x)$$ and $$b_5(x)$$ are zero within this element, all the terms in the fourth and fifth rows and columns are equivalently zero.

Similarly for the contribution from element 2:

Note: This matrix denotes the contribution due to element 2 to the global stiffness matrix. Since the basis functions $$b_1(x)$$ and $$b_2(x)$$ are zero within this element, all the terms in the first and second rows and columns are equivalently zero.

Now, constructing the global stiffness matrix using equation [[media:Fe1.s11.mtg30.djvu|(2)p.30-4]]

Computing the Element Stiffness Matrices
For the purpose of comparison we will now compute the global stiffness matrix by computing the element stiffness matrices and then using scatter matrices to scatter the coefficients into the proper place in the global stiffness matrix.

From equation [[media:Fe1.s11.mtg31.djvu|(1)p.31-3]]

A scatter matrix is a matrix that relates the local element degrees of freedom to the global degrees of freedom as seen in equation [[media:Fe1.s11.mtg31.djvu|(3)p.31-1]]

Due to the similarity of the structure of the element stiffness matrix and the element stiffness contribution matrix previously calculated, it can be seen that the element stiffness matrices are equivalent to the non-zero sub-matrices of the stiffness contribution matrices. i.e.

And by inspection the scatter matrices are found to be

Alternate Computation of the Global Stiffness Matrix
The element contribution to the global stiffness matrix is then computed via equation [[media:Fe1.s11.mtg31.djvu|(2)p.31-1]]

Now, as before,

Note:    $$\displaystyle (Eq. 6.2.21) = (Eq. 6.2.11)$$

QLEBF Used for nel=3
The following MATLAB code was used to generate the QLEBF used for nel=3, shown in figure 2 below

Numerical Solution
This problem was then solved using the QLEBF and increasing the number of elements, $$ nel = 2,4,6...$$ until the error converged below 10-6. The following MATLAB code was used to develop figures 3 and 4 below, showing a comparison between the numerical and exact solution and the error convergence, respectively.

Solved by William Kurth.

=Problem 6.9: Reproduce F&B Examples 8.3, 8.4 and Solve F&B Problem 8.6=

Objectives
We were asked to use a finite element analysis software or routine to:

1.) Reproduce example 8.3, FB, p. 200

2.) Reproduce example 8.4, FB, p. 205

3.) Use the same codes to solve problem 8.6, FB, p. 212

Solution
The MATLAB finite element codes developed by F&B, corresponding to chapter 8, were used to reproduce and solve the following problems, these codes are available at | the FB companion site and are also shown below.

The basic procedure for using these codes is as follows:

 * define the problem parameters, such as material properties, boundary conditions, mesh and element structure in an appropriate input file (i.e. input_file_64ele.m)
 * specify desired plots created in this input file
 * define the problem geometry, nodal coordinates, and element connectivity matrix either in this input file or in the mesh2d.m subroutine
 * specify in the preprocessor.m file which input file to load
 * if any new variables are introduced, add them to the include_flags.m file to establish them as global variables
 * run the heat2d.m file to solve the problem and produce the plots specified in the input file

The function assembly.m assembles passed element stiffness and force matrices into the current global stiffness and force matrices. Called iteratively for each element.

BmatHeat2D.m computes the basis function derivative matrix.

gauss.m gets Gauss points and weight factors.

get_flux.m calculates the flux vectors.

heat2d.m calls the necessary functions and .m files to carry out the computations.

heat2Delem.m calculates the element stiffness and force matrices as well as the elemental source term.

include_flags.m initializes global variables

input_file_1ele.m is the input file configured to compute example 8.3 with 1 element.

input_file_16ele.m is the input file containing parameters to solve example 8.3 with 16 elements.

input_file_64ele.m is the input file containing parameters to solve example 8.3 with 64 elements.

mesh2d.m constructs the nodes and mesh for example 8.3.

NmatHeat2D.m constructs the matrix of basis functions.

plotmesh.m is used to plot the mesh/element structure.

postprocess.m displays the calculated results.

preprocessor.m specifies which input file to load.

setup_ID_LM.m sorts the nodes which are constrained on the essential boundary condition from those which are not.

solvedr.m partitions the K and F matrices and solves for the degrees of freedom.

src_and_flux.m computes the nodal flux vectors.

Example 8.3
Consider a plate with geometry shown in figure 1 below. The conductivity is isotropic with $$D = k\bigl[\begin{smallmatrix} 1&0\\0&1\end{smallmatrix}\bigr]$$ with k = 5 W°C-1. The temperature is prescribed at T=0 at the bottom and right edges. Heat fluxes $$ \bar{q} = 0$$ and $$ \bar{q} = 20 Wm^{-1}$$ on the right and top boundaries, respectively. A constant heat source of $$\displaystyle s = 6 Wm^{-2}$$ is applied over the plate. Solve the problem using 16 quadrilateral finite elements.

The accompaniment codes come already configured to solve this problem, one must only simply specify in the preprocessor.m file that the input_file_16ele.m is the input file to use (comment out the others and save the file). Running heat2d.m displays these results:

Example 8.4
Consider a plate with geometry shown in figure 4 below, with prescribed temperature at x = ± b and prescribed flux at y = ± b.

The essential boundary conditions are given by

The natural boundary conditions on $$\Gamma_q$$ are $$(\bar{q}=-k\mathbf{n}^T\nabla T)$$

Lastly the source term is defined as

The MATLAB code used for example 8.3 was adapted to solve this problem by rewriting the input file. The results are shown below.

Problem 8.6
Consider a chimney constructed of two isotropic materials: dense concrete (k = 2.0 W°C-1 ) and bricks (k = 0.9 W°C-1 ). The temperature of the hot gasses on the inside surface of the chimney is 140 °C, whereas the outside is exposed to the surrounding air, which is at T = 10 °C. The dimensions of the chimney (in meters) are shown below. For the analysis, exploit the symmetry and consider 1/8 of the chimney crosssectional area. Consider a mesh of eight elements as shown below. Determine the temperature and flux in the two materials. Analyze the problem with 2 x 2, 4 x 4 and 8 x 8 quadrilateral elements for 1/8 of the problem domain. A 2 x 2 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.

The previous MATLAB code was once again adapted to solve the given problem. Below are the results for the cases of nel = 4, 16, and 64.

Solved by William Kurth.

title: User:Eml5526.s11.team6.kurth/HW6


text:
User:Egm6341.s11.team4.fields/HW6.9

=Problem 6.9: Identify the Basis Functions Ni(s) for i = 1,2,3,4 and Plot Them=

Refer to lecture slide [[media:nm1.s11.mtg37.djvu|37-3]] for the problem statement.

Given

 * {| style="width:70%" border="0" align="center"

\begin{Bmatrix} c_0 \\ c _1 \\ c _2 \\ c_3 \end{Bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -3 & -2 & 3 & -1 \\ 2 & 1 & -2 & 1 \end{bmatrix} \begin{Bmatrix} Z_i \\ Z_i^' \\ Z_{i+1} \\ Z_{i+1}^' \end{Bmatrix} =\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -3 & -2 & 3 & -1 \\ 2 & 1 & -2 & 1 \end{bmatrix} \begin{Bmatrix} \overline{d}_1 \\ \overline{d}_2 \\ \overline{d}_3 \\ \overline{d}_4 \\ \end{Bmatrix} $$     (9.1)
 * $$\displaystyle
 * $$\displaystyle
 *
 * }
 * }

Where:

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

\begin{Bmatrix} Z_i \\ Z_i^' \\ Z_{i+1} \\ Z_{i+1}^' \end{Bmatrix} = \begin{Bmatrix} \overline{d}_1 \\ \overline{d}_2 \\ \overline{d}_3 \\ \overline{d}_4 \\ \end{Bmatrix} = \begin{Bmatrix} d_1 \\ hd_2 \\ d_3 \\ hd_4 \end{Bmatrix} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

And equation (1) on lecture slide [[media:nm1.s11.mtg37.djvu|37-1]]:

 * {| style="width:70%" border="0" align="center"

$$     (9.2)
 * $$\displaystyle Z(s) = \sum_{i=0}^3 c_i s^i = \sum_{i=1}^4 \overline{N}_i(s)d_i
 * $$\displaystyle Z(s) = \sum_{i=0}^3 c_i s^i = \sum_{i=1}^4 \overline{N}_i(s)d_i
 *
 * }
 * }

Objective
a.) Identify the Basis Functions Ni(s) for i = 1,2,3,4

b.) Plot the Basis Functions Ni(s) for i = 1,2,3,4

a.) Basis Functions Ni(s) for i = 1,2,3,4
Expanding the Matrix equation given in (9.1) we get:

 * {| style="width:70%" border="0" align="center"

$$     (9.3)
 * $$\displaystyle c_0 = Z_i
 * $$\displaystyle c_0 = Z_i
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.4)
 * $$\displaystyle c_1 = Z_i^'
 * $$\displaystyle c_1 = Z_i^'
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.5)
 * $$\displaystyle c_2 = -3Z_i - 2Z_i^' + 3Z_{i+1} - Z_{i+1}^'
 * $$\displaystyle c_2 = -3Z_i - 2Z_i^' + 3Z_{i+1} - Z_{i+1}^'
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.6)
 * $$\displaystyle c_3 = 2Z_i + Z_i^' - 2Z_{i+1} + Z_{i+1}^'
 * $$\displaystyle c_3 = 2Z_i + Z_i^' - 2Z_{i+1} + Z_{i+1}^'
 *
 * }
 * }

By expanding equation (9.2) and then plugging the above values for Ci into the expanded equation we obtain the following:

 * {| style="width:70%" border="0" align="center"

$$     (9.7)
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4 = C_0s^0 + C_1s^1 + C_2s^2 + C_3s^3
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4 = C_0s^0 + C_1s^1 + C_2s^2 + C_3s^3
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.8)
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

= Z_i + (Z_i^')s^1 + (-3Z_i - 2Z_i^' + 3Z_{i+1} - Z_{i+1}^')s^2 + (2Z_i + Z_i^' - 2Z_{i+1} + Z_{i+1}^')s^3 $$
 * $$\displaystyle
 * $$\displaystyle

 * }
 * }

Rearranging the right side of equation (9.8) we get the following:

 * {| style="width:70%" border="0" align="center"

$$     (9.9)
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 * $$\displaystyle N_1(s)d_1 + N_2(s)d_2 + N_3(s)d_3 + N_4(s)d_4
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

= (1 -3s^2 + 2s^3)Z_i + (s - 2s^2 + s^3)z_i^' + (3s^2 - 2s^3)Z_{i+1} + (-s^2 + s^3)Z_{i+1}^' $$
 * $$\displaystyle
 * $$\displaystyle

 * }
 * }

By comparing the left and right sides of equation (9.9) we obtain the following basis functions:

 * {| style="width:70%" border="0" align="center"

$$     (9.10)
 * $$\displaystyle N_1(s) = 1 -3s^2 + 2s^3
 * $$\displaystyle N_1(s) = 1 -3s^2 + 2s^3
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.11)
 * $$\displaystyle N_2(s) = s - 2s^2 + s^3
 * $$\displaystyle N_2(s) = s - 2s^2 + s^3
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.12)
 * $$\displaystyle N_3(s) = 3s^2 - 2s^3
 * $$\displaystyle N_3(s) = 3s^2 - 2s^3
 *
 * }
 * }

 * {| style="width:70%" border="0" align="center"

$$     (9.13)
 * $$\displaystyle N_4(s) = -s^2 + s^3
 * $$\displaystyle N_4(s) = -s^2 + s^3
 *
 * }
 * }

b.) Plots of the Basis Functions Ni(s) for i = 1,2,3,4
The following Matlab code was used to create the figure below.

>> s = linspace(0,1,100); >> N1 = 1 - 3.*s.^2 + 2.*s.^3; >> N2 = s - 2.*s.^2 + s.^3; >> N3 = 3.*s.^2 - 2.*s.^3; >> N4 = -s.^2 + s.^3; >> plot(s,N1,s,N2,s,N3,s,N4)

title: User:Egm6341.s11.team4.fields/HW6.9


text: File:Fe1.s11.team6.kurth.hw8.prob2.fig1.jpg


text: File:S.11.Team6.tupsakhare-Figure- 6.3.5-(2).png


text:
MediaWiki:Wikimedia-copyrightwarning

By saving changes, you are agreeing to the Terms of Use and the Privacy Policy, and you irrevocably agree to release your contribution under the [//creativecommons.org/licenses/by-sa/4.0/ CC-BY-SA 4.0 License] and the [//www.gnu.org/copyleft/fdl.html GFDL]. You agree that a hyperlink or URL is sufficient attribution under the Creative Commons license.

title: MediaWiki:Wikimedia-copyrightwarning