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)