User:Egm6322.s12.team2/HW12

From Wikiversity
Jump to: navigation, search

R 12.1 - Integration Points and Gauss-Legendre Quadrature[edit]

Given[edit]

Legendre series

Find[edit]

Compute the 1st non-zero coefficient in the Legendre series with increasing number of integration points until convergence to within 1% accuracy. Use the error formula (1) p. 45b-7 to estimate the necessary number of integration points to achieve the desired accuracy.

Solution[edit]

Withn the class context, non-zero coefficients have generally been solved by applying the Fourier-Legendre series and orthogonality to find non-zero coefficients. For team 2, the previous solutions of 9.1 and 10.7 will be expanded to solve this problem. The general pattern is described below.

Step 1: Define 1st Non-Zero Coefficient

The coefficient of a Fourier-Legendre series requires the definition below.

A_j = \frac{2j+1}{2} \int_{-1}^1 P_j(x)T_0 \sinh{(1 - x^2)}d(x)

The coefficient is from the term for the Fourier-Legendre series of f(x).

f(x) = \sum_{n=0}^{\infty}a_n P_n(x)

Essentially, the series utilizes the Legendre polynomials and coefficients to define a function  f(x). The Legendre polynomials have the following orthogonality.

\int_{-1}^{1} P_n(x) P_m(x)dx = \frac{2}{2m+1} \delta_{mn}

The problem does not explicitly give a function, f(x), for finding the 1st non-zero coefficient. Therefore, the team's approach is to utilize the function from 9.1 and 10.7 below.

f(x) = T_0 \sinh(1-x^2)

The coefficients are then defined as follows.

A_j = \frac{2j+1}{2} \int_{-1}^1 P_j(x)T_0 \sinh{(1 - x^2)}d(x)

From 9.1 and 10.7, the 1st non-zero coefficient of the chosen f(x) is

A_0 = 0.7459T_0

Step 2: Apply Qauss-Quadrature to Integral of Coefficient

A_j = \frac{2j+1}{2} \int_{-1}^1 P_j(x)T_0\sinh{(1 - x^2)}d(x) \approx I_n(f) \sum_{j=1}^n w_jf(x_j)

The weighting function is

w_j = \frac{-2}{(n+1) P_n'(x_j)P_{n+1}(x_j)}

The error formula is shown below from 45-7b.

E_n(f) = \frac{2^{2n+1}(n!)^4}{(2n+1)[(2n)!]^2} \frac{f^{(2n)}(\eta)}{(2n)!}, \eta \in [-1,1]

Step 3: Find Number of Integration Points

The following definition is used to find the number of integration points.

I(f) = I_n(f) + E_n(f)

This can be rearranged to apply the error limit of 1% or 0.01.

\begin{alignat}{2}
                                     \frac{I(f) - I_n(f)}{I(f)} & = \frac{E_n(f)}{0.7459} \\
                                     & = \frac{1}{0.7459}\frac{2^{2n+1}(n!)^4}{(2n+1)[(2n)!]^2} \frac{f^{(2n)}(\eta)}{(2n)!} <= 0.01
\end{alignat}

Note that f^{(2n)}(\eta) is the (2n)th order derivative of f(x)=\sinh({1-x^2)}. Normally, this higher order derivative would require a rigorous analysis with numerical methods and sophisticated programming. But that can be minimized by noting that a derivative of sinh(x) = cosh(x) and vice versa. So the following definition is utilized.

\begin{alignat}{2}
                                    f^{(2n)}(x) & = \frac{d^{2n} \sinh({(1-x^2)}}{dx^{(2n)}} \\
                                                   & = {(-2x)}^{(2n)}\sinh{(1-x^2)}
\end{alignat}

So the simplified error term is below.

\frac{E_n(f)}{0.7459} = \frac{1}{0.7459}\frac{2^{2n+1}(n!)^4}{(2n+1)[(2n)!]^2} \frac{{(-2x)}^{(2n)}\sinh{(1-x^2)}}{(2n)!}

The Matlab image and code are shown below.

Matlab code Matlab code http://en.wikiversity.org/wiki/File:PEA2_12_1_Matlab.PNG

http://en.wikiversity.org/wiki/File:PEA2_12_1_Error_Term.png

See images below. The data from the Excel shows that there are two integration points needed to get the error less than 0.01 or 1%. The error term derived was simply entered into MS exceo with \eta = 0.5 and calculated for n = 1 to 20. At n =2 E < 0.01, so n = 2 will suffice.

http://en.wikiversity.org/wiki/File:Team2_12_1_Excel_Data.PNG

Egm6322.s12.team2.steele.m2 (talk) 20:32, 28 March 2012 (UTC)

R 12.2 - Alternative formula for weights of GL quad[edit]

Given[edit]



Using the following fomula

\displaystyle (1-x^{2})P_{n}^{'}=(n+1)xP_{n}-(n+1)P_{n+1}

(12.2.1)

Find[edit]



Show that

\displaystyle (n+1)P_{n+1}(x_{j})=-(1-x_{j}^{2})P_{n}^{'}(x_{j})

(12.2.2)

and

\displaystyle w_{j}=\frac{2}{(1-x_{j}^{2})[P_{n}^{'}(x_{j})]^{2}}

(12.2.3)

Solution[edit]

Solved on my own



From the 12.2.1, we have

\displaystyle (1-x_{j}^{2})P_{n}^{'}(x_{j})=(n+1)x_{j}P_{n}(x_{j})-(n+1)P_{n+1}(x_{j})

(12.2.4)

Since \displaystyle P_{n}(x_{j})=0,(The proof is @ p.45-17) then

\displaystyle (n+1)P_{n+1}(x_{j})=-(1-x_{j}^{2})P_{n}^{'}(x_{j})

(12.2.5)

In terms of (5)p.45b-6

\displaystyle w_{j}=\frac{-2}{(n+1)P_{n+1}(x_{j})P_{n}^{'}(x_{j})}

(12.2.6)

then substitute 12.2.5 into 12.2.6, we obtain

\displaystyle w_{j}=\frac{2}{(1-x_{j}^{2})[P_{n}^{'}(x_{j})]^{2}}

(12.2.3)



R 12.3 - Proof of Thean of Smynna[edit]

Given[edit]




   y_n =2x_{n-1} + y_{n-1}

(12.3.1)



   x_n =x_{n-1} + y_{n-1}

(12.3.2)



   \lim_{n\to\infty}\frac{y_n}{x_n}=\sqrt2

(12.3.3)

x_0 = 1

y_0 = 1

Find[edit]



Plot the sequence \{y_n/x_n, \ n=0,1,...10\} versus the iteration number n

Solution[edit]



By using Eq(12.3.1) and Eq(12.3.2),

Finding {x}_{0} to {x}_{10} and {y}_{0} to {y}_{10}

 n=0, y_{0}=1, x_{0}=1, \frac{y_{0}}{x_{0}}=1
 n=1, y_{1}=3, x_{1}=2, \frac{y_{1}}{x_{1}}=1.500
y_{1}=2x_{0}+y_{0}= 2+1=3 and x_{1}=x_{0} +y_{0}=1+1=2
 n=2, y_{2}=7, x_{2}=5, \frac{y_{2}}{x_{2}}=1.400
y_{2}=2x_{1}+y_{1}= (2)(2)+3=7 and x_{2}=x_{1} +y_{1}=2+3=5
 n=3, y_{3}=17, x_{3}=12, \frac{y_{3}}{x_{3}}=1.4167
y_{3}=2x_{2}+y_{2}= (2)(5)+7=17 and x_{3}=x_{2} +y_{2}=2+3=5
 n=4, y_{4}=41, x_{4}=29, \frac{y_{4}}{x_{4}}=1.4138
y_{4}=2x_{3}+y_{3}= (2)(12)+17=41 and x_{4}=x_{3} +y_{3}=12+17=29
 n=5, y_{5}=99, x_{5}=70, \frac{y_{5}}{x_{5}}=1.4143
y_{5}=2x_{4}+y_{4}= (2)(29)+41=99 and x_{5}=x_{4} +y_{4}=29+41=70
 n=6, y_{6}=239, x_{6}=169, \frac{y_{6}}{x_{6}}=1.4142
y_{6}=2x_{5}+y_{5}= (2)(70)+99=239 and x_{6}=x_{5} +y_{5}=70+99=169
 n=7, y_{7}=577, x_{7}=408, \frac{y_{7}}{x_{7}}=1.4142
y_{7}=2x_{6}+y_{6}= (2)(169)+239=577 and x_{7}=x_{6} +y_{6}=169+239=408
 n=8, y_{8}=1393, x_{8}=985, \frac{y_{8}}{x_{8}}=1.4142
y_{8}=2x_{7}+y_{7}= (2)(408)+577=1393 and x_{8}=x_{7} +y_{7}=408+577=985
 n=9, y_{9}=3363, x_{9}=2378, \frac{y_{9}}{x_{9}}=1.4142
y_{9}=2x_{8}+y_{8}= (2)(985)+1393=3363 and x_{9}=x_{8} +y_{8}=985+1393=2378
 n=10, y_{10}=8119, x_{10}=5741, \frac{y_{10}}{x_{10}}=1.4142
y_{10}=2x_{9}+y_{9}= (2)(2378)+3363=8119 and x_{10}=x_{9} +y_{9}=2378+3363=5741

Plot the Graph
File:Ssskimmm1.png Egm6322.s12.sungsik (talk) 21:08, 27 March 2012 (UTC)

R 12.4 - Computing the convergents of \sqrt 2[edit]

Given[edit]


Continued fraction formula for \sqrt 2:


   \displaystyle 
      \sqrt 2 = 1 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2 + {2 + \frac{1}{2 + \frac{1}{\ddots}}}}}}}}

(12.4.1)



First four convergents:


   \displaystyle 
    \begin{align}
      C_0 &:= 1\\
      C_1 &:= 1 + \frac{1}{2} = \frac{3}{2}\\
      C_2 &:= 1 + \frac{1}{2 + \frac{1}{2}} = \frac{7}{5}\\
      C_3 &:= 1 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2}}} = \frac{17}{12}\\
    \end{align}

(12.4.2)


Find[edit]

Write a Matlab program to compute the convergents C_i for i = 1,2,...10 and plot these convergents along with the line \sqrt 2 to visualize the convergence \lim_{i \to\infty} C_i = \sqrt 2

Solution[edit]

Solved on my own


This problem mainly deals with Matlab implementation. From Eq. 12.4.1 it is clear that the quotient of the terms in the denominator is always 2 due to the relation \sqrt 2 = 1 + (\sqrt 2 - 1). It is seen in the Mtg 48e that the term on the right can be expressed as \frac{1}{\sqrt 2 + 1}; substituting in the first relation then produces the quotient of 2.

When implementing this code in Matlab, it is more convenient to 'build up' the residue term starting with the i_{max} subtraction. As all contributions to the evaluation after i_{max} are ignored, this value is simply 1/2. It also defines the final level of the denominator which must be computed. Thus, a for loop is constructed between [i_{max}:-1:2] to iteratively develop the residue term. This is seen in the Matlab code in lines 11 - 14. The computation itself is quite simple; each successive iteration simply adds the quotient of 2 onto the previously evaluated derivative.

The convergents are computed in Line 19 by adding the first quotient of 1 onto each residue term. The final values are then plotted against the exact value (to machine precision) of \sqrt 2 and are seen to converge within 5% of the exact value as soon as the second convergent. It can also be seen that the first 4 terms match those given by Eq. 12.4.2.

Plot of continued fraction approximations of

Plot of continued fraction approximations of

Matt Shields 11:40,24 March 2012 (UTC)

R 12.5 Legendre polynomials and continued fraction[edit]

Given[edit]

\frac{a_1}{x+}\frac{-\frac{1}{2}}{\frac{3}{2} x +} \frac{-\frac{3}{4}}{\frac{7}{4} x +} ... \frac{-1\frac{2(n-1)-1}{2(n-1)}}{\frac{4(n-1)-1}{2(n-1)} x +} ...

Find[edit]

(i) Verify the last term for the nth approximant as shown above.
(ii) The Legendre polynomials \{P_1, ..., P_6\} and compare the results to those in R10.2 on p. 44b-6.
(iii) What can you observe about the factor in P_n(x) when generated using the contnued fraction compared to the factor of P_x(x) as obtained in R10.2?
(iv) Would the difference in these factors change the orthogonality of the polynomials generated by the continued fraction (1) p. 48-24b?

Solution[edit]

Part (i)

The solution requires an analysis of the continued fraction given for this problem.

\frac{a_1}{x+}\frac{-\frac{1}{2}}{\frac{3}{2} x +} \frac{-\frac{3}{4}}{\frac{7}{4} x +} ... \frac{-1\frac{2(n-1)-1}{2(n-1)}}{\frac{4(n-1)-1}{2(n-1)} x +} ...

There is a general pattern of continued fractions.

c_1 = b_0 + \frac{a_1}{b_1} = 0 + \frac{a_1}{b_1}

c_2 = \frac{a_1}{b_1 + \frac{a_2}{b_2}}

c_3 = \frac{a_1}{b_1 + \frac{a_2}{b_2 + \frac{a_3}{b_3}}}

c_4 = \frac{a_1}{b_1 + \frac{a_2}{b_2 + \frac{a_3}{b_3 + \frac{a_4}{b_4}}}}

....

c_n = {\frac{a_1}{b_1} +} {\frac{a_2}{b_2} +} {\frac{a_3}{b_3} +}...{\frac{a_n}{b_n}}

In this case

a_n = {-1\frac{2(n-1)-1}{2(n-1)}}

b_n = \frac{4(n-1)-1}{2(n-1)}x

Part (ii)

By applying the formula from part(i), we get

c_1 = \frac{a_1}{x}

Therefore, the 1st Legendre polynomial is

P_1(x) = x

Next, we get the 2nd Legendre polynomial.

\begin{alignat}{3}
                                    C_2(x) & = \frac{a_1}{x + \frac{-1}{2*32x}} \\
                                           & = \frac{a_1}{x - \frac{1}{3x}} \\
                                           & = \frac{3a_1}{3x^2 - 1}
\end{alignat}

Therefore, the 2nd Legendre polynomial is

P_2(x) = 3x^2 - 1

The 3rd Legendre polynomial is just found by the 3rd convergent.

C_2(x) = \frac{a_1}{x + \frac{-0.5}{32x + \frac{-0.75}{{\frac{7}{4}} x}}}

The inverse to get the Legendre polynomial is (5x^3 - 3x). From this pattern, it is clear that for each Legendre polynomial the related factor is missing when compared to 10.6. To confirm this pattern, let us rely on these two sources.

http://en.wikipedia.org/wiki/Legendre_polynomials

http://digitalcommons.uconn.edu/cgi/viewcontent.cgi?article=1079&context=chem_educ

C.W. David of University of Connecticut above demonstrates that the Legendre polynomials are derived from continued fractions in such a way that the constant of the solution is not necessary since the derived polynomial fits the general Legendre differential equation as a solution.

P_4(x) = (35x^4 - 30x^2 + 3)

P_5(x) = (63x^5 - 70x^3 + 15x)

P_6(x) = (231x^6 - 315x^4 + 105x62 - 5)

Part (iii)

In part (ii) we found that the fractional term at the front of each Legendre component is absent in the solutions found by continued fractions. For example, the fraction (1/2) is missing from the 2nd and 3rd Legendre polynomials for continued fractions. The method of derivation based on continued fractions produce solutions based on the derivates of the Legendre differential equations and recursive insertion of solutions. The ratio of derivatives results in an integration to give the Legendre function.

Part (iv)

The orthogonal property is not dependent on the factor term. The orthogonality is based on this term below, so the solutions from continuing fractions with convergents are still orthonal. This solution is not dependent on the fractional term derived from Rodrigues formula. The absence of the fractional term does not affect orthongality, which is still valid for continued fractional solutions.

\int_{-1}^{1}P_n(x)P_m(x)dx = \frac{2}{2m+1}\delta_{mn}


Egm6322.s12.team2.steele.m2 (talk) 20:30, 28 March 2012 (UTC)

Contributing Members[edit]

Team Contribution Table
Problem Number Assigned To Solved By Typed By Proofread By
12.1 Manuel Steele Manuel Steele Manuel Steele Lang Xia
12.2 Lang Xia Lang Xia Lang Xia Manuel Steele
12.3 sungsik kim sungsik kim sungsik kim Matt Shields
12.4 Matt Shields Matt Shields Matt Shields Lang Xia
12.5 Manuel Steele Manuel Steele Manuel Steele sungsik kim



References for Report 12[edit]

King, A.C., J. Billingham and S.R. Otto. "Differential Equations: Linear, Nonlinear, Ordinary, Partial." New York, NY: Cambridge University Press, 2003.

http://en.wikipedia.org/wiki/Legendre_polynomials

http://digitalcommons.uconn.edu/cgi/viewcontent.cgi?article=1079&context=chem_educ

Vu-Quoc, L. Class Lecture: Principles of Engineering Analysis. University of Florida, Gainesville, FL, (Meeting 45h) Mtg 45h Spring 2012.

Vu-Quoc, L. Class Lecture: Principles of Engineering Analysis. University of Florida, Gainesville, FL, (Meeting 48a) Mtg 48a Spring 2012.

Vu-Quoc, L. Class Lecture: Principles of Engineering Analysis. University of Florida, Gainesville, FL, (Meeting 48b) Mtg 48b Spring 2012.

Vu-Quoc, L. Class Lecture: Principles of Engineering Analysis. University of Florida, Gainesville, FL, (Meeting 48c) Mtg 48c Spring 2012.

Vu-Quoc, L. Class Lecture: Principles of Engineering Analysis. University of Florida, Gainesville, FL, (Meeting 48d) Mtg 48d Spring 2012.

Vu-Quoc, L. Class Lecture: Principles of Engineering Analysis. University of Florida, Gainesville, FL, (Meeting 48e) Mtg 48e Spring 2012.

http://mathworld.wolfram.com/Fourier-LegendreSeries.html