text:
User:Eml5526.s11.team6.gravois

Team Homework Report Page Homework 1 Homework 2 Homework 3 Homework 4 = Homework 1 =

Find
Derive the following equation by performing a balance of forces:

$$ \frac{\partial}{\partial x}\left [ A(x)E(x)\frac{\partial u}{\partial x} \right ]+f(x,t)=A(x) \rho(x)\cdot \frac{\partial ^{2}u}{\partial t^{2}}\ $$

Solution
Consider the FBD of an infinitesimally small section of length dx of the 1-D elastic bar with varying material properties at a distance x from the origin. Where, A(x) - Varying cross section of the elastic bar, E(x) - Modulus of elasticity, rho (x) - Mass density, u - Displacement constrained only to x- axis (1-D bar), t - Time.

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

Where, $$\; {A}'=A(x) (dx)\cdot \rho (x)\cdot\frac{\partial^{2}u }{\partial t^{2}} =mass\cdot acceleration=Inertia\, Force$$
 * $$\displaystyle (Eq. 1.1) $$
 * }
 * }

As per Newton's law of motion, the involved forces should contribute to inertia force.

Force equilibrium in the x-direction yields:

$$\sum F_{x}=0 $$

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

$$\sum F_{x}=-N(x,t)+N(x+dx,t)+f(x,t) dx-m(x)\cdot \frac{\partial ^{2}u}{\partial t^{2}}=0$$
 * $$\displaystyle (Eq. 1.2) $$
 * }
 * }

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

Where, $$\;\; m(x)=A(x)(dx)\rho (x)$$
 * $$\displaystyle (Eq. 1.3) $$
 * }
 * }

Including Taylor series expansion for $$N(x+dx,t)$$ and neglecting the higher order terms (h.o.t) yields:

$$ \cancel {-N(x,t)}+[\cancel {N(x,t)}+\frac {\partial }{\partial x}(N(x,t))+(h.o.t)]+f(x,t)=\frac{m(x)}{dx}\cdot\frac{\partial ^{2} u}{\partial t^{2}}$$

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

$$\frac{\partial }{\partial x}(N(x,t))+f(x,t)=\frac{m(x)}{dx}\cdot\frac{\partial ^{2} u}{\partial t^{2}}$$
 * $$\displaystyle (Eq. 1.4) $$
 * }
 * }

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

$$Also,\; N(x,t)=A(x) \cdot \sigma (x)=A(x) \cdot [E(x)\cdot \epsilon (x,t)] $$
 * $$\displaystyle (Eq. 1.5) $$
 * }
 * }

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

$$\;\;\;\;\;\;\;\;\;\;\;\Rightarrow N(x,t)= A(x)\cdot E(x) \cdot\frac{\partial u}{\partial x}(x,t)$$
 * $$\displaystyle (Eq. 1.6) $$
 * }
 * }

Thus proving:

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

 * $$\displaystyle (Eq. 1.7) $$
 * }
 * }

Find
Discuss the case in which the bar has a rectangular cross section shown below.

Solution
We have from reference of class notes for mtg 2, eq(3) pg.2-3,

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

$$m=\rho (x) \cdot A(x) $$
 * $$\displaystyle (Eq. 2.1) $$
 * }
 * }

Also with reference to class notes for mtg 6, pg 6-1, point 2), for A(x)=bˑh(h),

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

$$m=\rho \left ( x+\frac{dx}{2} \right )\cdot \frac{1}{2}\left [ h(x)+h(x+dx) \right ]b $$
 * $$\displaystyle (Eq. 2.2) $$
 * }
 * }

Substituting Eq 2.1 and Eq 2.2 into Eq 1.7 for case of A(x)=bˑh(x),
 * {| style="width:100%" border="0"

$$\frac{\partial}{\partial x}\left [ \left ( b\cdot h(x)\cdot E(x)\cdot \frac{\partial u}{\partial x}\right ) \right ]+f(x,t)=\frac{1}{2}\left [ h(x)+h(x+dx) \right ]\rho (x+dx)\cdot \frac{\partial^{2} u}{\partial t^{2}}$$
 * $$\displaystyle (Eq. 2.3) $$
 * }
 * }

Note:

1) $$\left [ \frac{1}{2}\left [ h(x)+h(x+dx) \right ]\rho (x+dx) \right ]$$ denotes mass per unit length at a distance of (x+dx/2), which is the middle of the elementary section considered in Figure 1.2.

2) $$\frac {1}{2}\left [h(x)+h(x+dx) \right ]$$ denotes average height of elementary section of figure 1.2.

3) Since we are considering an infinitesimally small section (dx) of the entire bar, we can write,

$$h(x)\approx h(x+dx)$$

$$\Rightarrow \frac{1}{2}\left [ h(x)+h(x+dx) \right ]\approx h(x)$$

4) Also, b is the width of the cross section and is constant and ρ(x+dx) is approximately equal to ρ(x) for a small section as in 3) above.

Eq 2.3 converts to the following:

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

 * $$\displaystyle (Eq. 2.4) $$
 * }
 * }

Eq 2.4 is the modified form of Eq 1.7 for the specific case of A(x)=bˑh(x).

Eml5526.s11.team6.gravois 04:37, 25 January 2011 (UTC)

= Homework 2 =

Find
1) Construct Γ5x5(F) and observe the proportionality of Γ.

2) Find det Γ(F)

3) Conclude F is orthogonal basis i.e. $$\Gamma _{ij}=\left \langle b_{i},b_{j} \right \rangle=\delta _{ij}$$

Given
$$F=\left \{ 1,cos(\omega x),cos(2\omega x),sin(\omega x),sin(2\omega x) \right \} , \Omega =[0,T], i=1,2$$

Solution
1) $$ \mathbf{\Gamma} _{5x5}=\int_{\Omega }{}b_{i}\cdot b_{j}\;dx=\int_{0}^{T}b_{i}\cdot b_{j}\;dx $$

Constructing a matrix of bi ˑ bj yields:

$$ \Rightarrow \begin{bmatrix} 1\cdot 1 & cos(\omega x)\cdot 1 & cos(2\omega x)\cdot 1 & sin(\omega x)\cdot 1 & sin(2\omega x)\cdot 1\\ 1\cdot cos(\omega x) & cos(\omega x)\cdot cos(\omega x) & cos(2\omega x)\cdot cos(\omega x) & sin(\omega x)\cdot cos(\omega x) & sin(2\omega x)\cdot cos(\omega x)\\ 1\cdot cos(2\omega x) & cos(\omega x)\cdot cos(2\omega x) & cos(2\omega x)\cdot cos(2\omega x) & sin(\omega x)\cdot cos(2\omega x) & sin(2\omega x)\cdot cos(2\omega x)\\ 1\cdot sin(\omega x) & cos(\omega x)\cdot sin(\omega x) & cos(2\omega x)\cdot sin(\omega x) & sin(\omega x)\cdot sin(\omega x) & sin(2\omega x)\cdot sin(\omega x)\\ 1\cdot sin(2\omega x) & cos(\omega x)\cdot sin(2\omega x) & cos(2\omega x)\cdot sin(2\omega x) & sin(\omega x)\cdot sin(2\omega x) & sin(2\omega x)\cdot sin(2\omega x) \end{bmatrix} $$

Multiplying out the above matrix and integrating the result using a table of integrals yielded the following matrix:

$$ \Rightarrow \begin{bmatrix} x & sin(\omega x) & sin(\omega x)cos(\omega x) & -cos(\omega x) & -\frac{1}{2}cos(2\omega x)\\ sin(\omega x) & \frac{1}{2}(x+sin(\omega x)cos(\omega x)) & \frac{1}{6}(3sin(\omega x)+sin(3\omega x)) & -\frac{1}{2}cos^2(\omega x) & -\frac{2}{3}cos^3(\omega x)\\ sin(\omega x)cos(\omega x) & \frac{1}{6}(3sin(\omega x)+sin(3\omega x)) & \frac{1}{8}(4x+sin(4\omega x)) & \frac{1}{6}(3cos(\omega x)-cos(3\omega x)) & -\frac{1}{8}cos(4\omega x))\\ -cos(\omega x)) & -\frac{1}{2}cos^2(\omega x) & \frac{1}{6}(3cos(\omega x)-cos(3\omega x)) & \frac{1}{2}(x-sin(\omega x)cos(\omega x)) & \frac{2}{3}sin^3(\omega x)\\ -\frac{1}{2}cos(2\omega x) & -\frac{2}{3}cos^3(\omega x) & -\frac{1}{8}cos(4\omega x)) & \frac{2}{3}sin^3(\omega x) & \frac{1}{8}(4x-sin(4\omega x)) \end{bmatrix} $$

By definition of Angular Frequency

$$2\pi =\omega T\rightarrow \omega=\frac{2\pi}{T}$$

Solving the above matrix over the limits of the bar from 0 to T yielded:

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

 * }
 * }

Г is a symmetric matrix about the diagonal.

2) The determinate of the above Г matrix was computed using MATLAB as shown.

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

 * }
 * }

3) An orthogonal matrix is a matrix such that the transpose of the matrix multiplied by the original matrix should yield the identity matrix I, a matrix with only 1's on the diagonal and zeros everywhere else.

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

 * }
 * }

Find
1) Construct Γ5x5(F) and observe the proportionality of Γ.

2) Find det Γ(F)

3) Conclude F is orthogonal basis i.e. $$\Gamma _{ij}=\left \langle b_{i},b_{j} \right \rangle=\delta _{ij}$$

Given
$$F=\left \{ 1,x,x^2,x^3,x^4 \right \},

\Omega =[0,1]$$

Solution
1) $$ \mathbf{\Gamma} _{5x5}=\int_{\Omega }{}b_{i}\cdot b_{j}\;dx=\int_{0}^{1}b_{i}\cdot b_{j}\;dx $$

Constructing a matrix of bi ˑ bj yields:

$$ \Rightarrow \begin{bmatrix} 1\cdot 1 & x\cdot 1 & x^2\cdot 1 & x^3\cdot 1 & x^4\cdot 1\\ 1\cdot x & x\cdot x & x^2\cdot x & x^3\cdot x & x^4\cdot x\\ 1\cdot x^2 & x\cdot x^2 & x^2\cdot x^2 & x^3\cdot x^2 & x^4\cdot x^2\\ 1\cdot x^3 & x\cdot x^3 & x^2\cdot x^3 & x^3\cdot x^3 & x^4\cdot x^3\\ 1\cdot x^4 & x\cdot x^4 & x^2\cdot x^4 & x^3\cdot x^4 & x^4\cdot x^4 \end{bmatrix} $$

Multiplying out the above matrix yields:

$$ \Rightarrow \begin{bmatrix} 1 & x & x^2 & x^3 & x^4\\ x & x^2 & x^3 & x^4 & x^5\\ x^2 & x^3 & x^4 & x^5 & x^6\\ x^3 & x^4 & x^5 & x^6 & x^7\\ x^4 & x^5 & x^6 & x^7 & x^8 \end{bmatrix} $$

Integrating the above matrix with respect to x yields:

$$ \Rightarrow \begin{bmatrix} x & \frac{x^2}{2} & \frac{x^2}{2} & \frac{x^3}{3} & \frac{x^4}{4}\\ \frac{x^2}{2} & \frac{x^3}{3} & \frac{x^4}{4} & \frac{x^5}{5} & \frac{x^6}{6}\\ \frac{x^3}{3} & \frac{x^4}{4} & \frac{x^5}{5} & \frac{x^6}{6} & \frac{x^7}{7}\\ \frac{x^4}{4} & \frac{x^5}{5} & \frac{x^6}{6} & \frac{x^7}{7} & \frac{x^8}{8}\\ \frac{x^5}{5} & \frac{x^6}{6} & \frac{x^7}{7} & \frac{x^8}{8} & \frac{x^9}{9} \end{bmatrix}_{0}^{1} $$

Evaluating the integrated matrix with the limits of [0,1] yielded:

$$ \Rightarrow \begin{bmatrix} 1-0 & \frac{1^2}{2}-0 & \frac{1^3}{3}-0 & \frac{1^4}{4}-0 & \frac{1^5}{5}-0\\ \frac{1^2}{2}-0 & \frac{1^3}{3}-0 & \frac{1^4}{4}-0 & \frac{1^5}{5}-0& \frac{1^6}{6}-0\\ \frac{1^3}{3}-0 & \frac{1^4}{4}-0 & \frac{1^5}{5}-0 & \frac{1^6}{6}-0& \frac{1^7}{7}-0\\ \frac{1^4}{4}-0 & \frac{1^5}{5}-0 & \frac{1^6}{6}-0 & \frac{1^7}{7}-0& \frac{1^8}{8}-0\\ \frac{1^5}{5}-0 & \frac{1^6}{6}-0 & \frac{1^7}{7}-0 & \frac{1^8}{8}-0& \frac{1^9}{9}-0 \end{bmatrix} $$

The family of linearly independent basis functions matrix is thus:

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

 * }
 * }

Г is a symmetric matrix about the diagonal.

2) The determinate of the above Г matrix was computed using MATLAB as shown.

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

 * }
 * }

3) An orthogonal matrix is a matrix such that the transpose of the matrix multiplied by the original matrix should yield the identity matrix I, a matrix with only 1's on the diagonal and zeros everywhere else.

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

 * }
 * }

Eml5526.s11.team6.gravois 01:01, 1 February 2011 (UTC)

=Homework 3=

Given

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

 * $$\alpha =b_i \left [ a_{2}'b_{j}'+a_{2}b_{j}''\right ]$$
 * $$\displaystyle (Eq. 4.1) $$
 * }
 * {| style="width:100%" border="0"
 * }
 * {| style="width:100%" border="0"

Where,
 * $$\beta =b_j \left [ a_{2}'b_{i}'+a_{2}b_{i}''\right ]$$
 * $$\displaystyle (Eq. 4.2) $$
 * }
 * }
 * }

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

$$\overline{K}_{ij}=\int_{\Omega }b_i (x)\left [ \frac{d}{dx}\left (a_2\frac{d}{dx}b_j (x) \right ) \right ]dx$$
 * $$\displaystyle (Eq. 4.3) $$
 * }
 * }

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

$$\overline{K}_{ji}=\int_{\Omega }b_j (x)\left [ \frac{d}{dx}\left (a_2\frac{d}{dx}b_i (x) \right ) \right ]dx$$
 * $$\displaystyle (Eq. 4.4) $$
 * }
 * }

α represents the integrand in the stiffness matrix (Eq 4.3) and β represents the integrand in stiffness matrix (Eq 4.4).

Find
Show that the above $$\alpha \neq \beta $$ using the following example:

$$\left\{\begin{matrix} b_i(x)=cos(ix)\\ b_j(x)=cos(jx) \end{matrix}\right.$$ where $$i \neq j$$

Solution
Substituting bi and bj into the α and β equations yielded:

$$\alpha = cos(ix)\left [ a_2'sin(jx)'+a_2cos(jx)'' \right ]$$ $$\beta = cos(jx)\left [ a_2'sin(ix)'+a_2cos(ix)'' \right ]$$

Carrying out the derivations yielded the following:

$$\alpha = cos(ix)\left [ -ja_2'sin(jx)-j^{2}a_2cos(jx) \right ]$$ $$\beta = cos(jx)\left [ -ia_2'sin(ix)-j^{2}a_2cos(ix) \right ]$$

Thus, $$\alpha \neq \beta$$

Since $$\alpha \neq \beta$$, it proves that the stiffness matrix $$\overline{K}$$ is not symmetric.

Given
The three-bar structure below is subjected to the given load of 103N at point B. Th Young's modulus is E=1011Pa, the cross sectional area of the bar BC is 2x10-2m2 and BD is 10-2m2. Point D is free to move in the x-direction and all distances are given in meters.

Find
a. Construct the global stiffness matrix and load matrix. b. Partition the matrices and solve for the unknown displacements at Point B and Point D. c. Find the stresses in the three bars. d. Find the reactions at nodes C, D, F.

a. Construct the global stiffness matrix and load matrix.
The stiffness matrix for each bar element is given below:

$$\mathbf{K}^{(n)}=\frac{AE}{l}\begin{bmatrix} cos(\phi)^2 & cos(\phi)sin(\phi) & -cos(\phi)^2 & -cos(\phi)sin(\phi)\\ cos(\phi)sin(\phi) & sin(\phi)^2 & -cos(\phi)sin(\phi) & -sin(\phi)^2\\ -cos(\phi)^2 & -cos(\phi)sin(\phi) & cos(\phi)^2 & cos(\phi)sin(\phi)\\ -cos(\phi)sin(\phi) & -sin(\phi)^2 & cos(\phi)sin(\phi) & sin(\phi)^2 \end{bmatrix}$$

For bar (1):

$$cos(\phi)=cos(90)=0$$ $$sin(\phi)=sin(90)=1$$

Substituting the above results into the general element stiffness matrix yielded the following:

$$\mathbf{K}^{(1)}=\frac{2AE}{l}\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & -1\\ 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 1 \end{bmatrix}$$

For bar (2):

$$cos(\phi)=cos(45)=\frac{1}{\sqrt{2}}$$ $$sin(\phi)=sin(45)=\frac{1}{\sqrt{2}}$$

Substituting the above results into the general element stiffness matrix yielded the following:

$$\mathbf{K}^{(2)}=\frac{AE}{\sqrt{2}l}\begin{bmatrix} \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2}\\ \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2}\\ -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & \frac{1}{2}\\ -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{bmatrix}$$

For bar (3):

$$cos(\phi)=cos(135)=-\frac{1}{\sqrt{2}}$$ $$sin(\phi)=sin(135)=\frac{1}{\sqrt{2}}$$

Substituting the above results into the general element stiffness matrix yielded the following:

$$\mathbf{K}^{(3)}=\frac{AE}{\sqrt{2}l}\begin{bmatrix} \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2}\\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}\\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}\\ \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \end{bmatrix}$$

Inserting the above element stiffness matrices into the global stiffness matrix yielded the following:

$$\mathbf{K}=\frac{AE}{l}\begin{bmatrix} 0+\frac{1}{2\sqrt{2}}+\frac{1}{2\sqrt{2}} & 0+\frac{1}{2\sqrt{2}}+\frac{1}{2\sqrt{2}} & 0 & 0 & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}\\ 0+\frac{1}{2\sqrt{2}}-\frac{1}{2\sqrt{2}} & 2+\frac{1}{2\sqrt{2}}+\frac{1}{2\sqrt{2}} & 0 & -2 & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -2 & 0 & 2 & 0 & 0 & 0 & 0\\ \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & 0\\ \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & 0\\ -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}}\\ -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} \end{bmatrix}$$

After combining terms in the global stiffness matrix, the result is as follows:

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

 * }
 * }

The general load matrix is given below:

$$\mathbf{F}=\begin{bmatrix} F_{1x} & F_{1y} & F_{2x} & F_{2y} & F_{3x} & F_{3y} & F_{4x} & F_{4y}\\ \end{bmatrix}^{T}$$

Substituting in the force acting at Point B yielded the load matrix below:

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

 * }
 * }

b. Partition the matrices and solve for the unknown displacements at Point B and Point D.
The Global systems of equations become:

F = K d

$$\begin{bmatrix} 1000 \\ 0 \\ r_{2x} \\ r_{2y} \\ r_{3x} \\ r_{3y} \\ 0 \\ r_4y\\ \end{bmatrix}=1x10^9\begin{bmatrix} \frac{1}{\sqrt{2}} & 0 & 0 & 0 & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}\\ 0 & 2+\frac{1}{\sqrt{2}} & 0 & -2 & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -2 & 0 & 2 & 0 & 0 & 0 & 0\\ \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & 0\\ \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & 0\\ -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}}\\ -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} \end{bmatrix}\cdot\begin{bmatrix} u_{1x}\\ u_{1y}\\ 0\\ 0\\ 0\\ 0\\ u_{4x}\\ 0

\end{bmatrix}$$

Partitioning to solve for the displacement of Point B yields the following matrices:

$$\begin{bmatrix} 1000\\0 \\0 \end{bmatrix}=10^9 \begin{bmatrix} \frac{1}{\sqrt2} & 0 & -\frac{1}{2\sqrt2}\\ 0 & 2+\frac{1}{\sqrt2} & -\frac{1}{2\sqrt2}\\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \end{bmatrix}\cdot \begin{bmatrix} u_{1x}\\u_{1y} \\u_{4x} \end{bmatrix}$$

$$\begin{bmatrix} 1000\\0 \\0 \end{bmatrix}= \begin{bmatrix} 7.071x10^8 & 0 & -3.536x10^8\\ 0 & 2.7071x10^9 & -3.536x10^8\\ -3.536x10^8 & -3.536x10^8 & 3.536x10^8 \end{bmatrix}\cdot \begin{bmatrix} u_{1x}\\u_{1y} \\u_{4x} \end{bmatrix}$$

Solving for the displacements yields:

$$\mathbf{d}=\begin{bmatrix} 3.328x10^{-6} \\ 5x10^{-7} \\ 3.828x10^{-6} \end{bmatrix}m$$

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

 * }
 * }

c. Find the stresses in the three bars.
The stress in each member is obtained using the following equation:

$$\sigma^{(1)}=\frac{E^{(n)}}{l^{(n)}}\begin{bmatrix} -cos(\phi^{(n)})& -sin(\phi^{(n)})& cos(\phi^{(n)})& sin(\phi^{(n)}) \end{bmatrix}\mathbf{d}^{(n)}$$

For Bar 1:

$$\sigma^{(1)}=\frac{10^9}{1}\begin{bmatrix} 0 & -1& 0& 1 \end{bmatrix}\mathbf\cdot \begin{bmatrix} 3.328x10^{-6}\\ 5x10^{-7}\\ 0\\ 0 \end{bmatrix}$$

For Bar 2:

$$\sigma^{(2)}=\frac{10^9}{1}\begin{bmatrix} -\frac{1}{\sqrt2} & -\frac{1}{\sqrt2}& \frac{1}{\sqrt2}& \frac{1}{\sqrt2} \end{bmatrix}\mathbf\cdot \begin{bmatrix} 3.328x10^{-6}\\ 5x10^{-7}\\ 3.828x10^{-6}\\ 0 \end{bmatrix}$$

For Bar 3:

$$\sigma^{(2)}=\frac{10^9}{1}\begin{bmatrix} \frac{1}{\sqrt2} & -\frac{1}{\sqrt2}& -\frac{1}{\sqrt2}& \frac{1}{\sqrt2} \end{bmatrix}\mathbf\cdot \begin{bmatrix} 3.328x10^{-6}\\ 5x10^{-7}\\ 0\\ 0 \end{bmatrix}$$

d. Find the reactions at nodes C, D, F.
$$\mathbf{r}=\begin{bmatrix} r_{2x}\\ r_{2y}\\ r_{3x}\\ r_{3y}\\ r_{4y}\\ \end{bmatrix}=\mathbf{Kd}=10^{9}\begin{bmatrix} 0 &0 & 0\\ 0 & -2 & 0\\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2}\\ -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0\\ \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 \end{bmatrix}\cdot \begin{bmatrix} 3.328x10^{-6}\\ 0.5x10^{-6}\\ 3.828x10^{-6} \end{bmatrix}= \begin{bmatrix} 0\\ -1000\\ 0\\ -1000\\ 1000 \end{bmatrix}N $$

= Homework 4 =

Given
Consider the elastic bar with a variable distributed spring p(x) along its length shown below in Figure 4.6.1. The distributed distributed s[pring imposes an axial force on the bar in proportion to the displacement. The bar is of length l, cross-sectional area A(x), Young's modulus E(x) with body force b(x) and boundary conditions shown in the following figure.

The above problem statement is given in A First Course in Finite Elements FB p74 Problem 3.9.

Find
a. Construct the strong form. b. Construct the weak form.

Solution
The bar must satisfies the following conditions: 1. The bar must be in equilibrium. 2. The bar must satisfy Hooke's Law σ(x)=E(x)ε(x). 3. The displacement field must be compatible. 4. It must satisfy the strain-displacement equation.

The following Figure 4.6.2 shows the internal force balance of a segment of the bar.

a. Strong Form
Performing a force balance on the section in Figure 4.6.2 yields the following differential equation 4.6.1.

Dividing the above equation through by dx to simplify terms yielded the following equation 4.6.2.

Now take the limit as Δx goes to 0. The result is presented below in equation 4.6.3.

By Hooke's Law:

Substituting Eq.4.6.6 into Eq.4.6.5 and then plugging that result into Eq.4.6.4 yielded the following equation.

Substituting Eq.4.6.7 into Eq.4.6.3 yields the final strong for differential equation as follows:

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

 * }
 * }

b. Weak Form
For the weak form, we multiply the strong form solution by the weighting function w and integrate the result over the length of the bar 0 to l as shown in the following equation:

Integrating the first term of Eq.4.6.8 by using integration by parts yielded the following equation:

Knowing that: $$A(x)E(x)\frac{du}{dx}=\bar{t}$$ and accounting for the b(x)-p(x) term yielded the following:

The final weak form becomes the following:

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

 * }
 * }

Figure
Following figure 4.7.8 shows the geometry and the mesh generated by CGX for a sphere volume segment

Geometry generation code
Following figure 4.7.10 shows the code to generate the geometry and the mesh by CGX for the cylinder in Figure 4.7.9.

Explanation of the commands used is given in the succeeding subsection.

Explanation of the commands used
Please refer to the 'Disc' section above for explanation of all the following CGX commands used to build the geometry and mesh of the cylinder Following commands are used to generate the geometry and mesh of the cylinder in figure 4.7.3 -

' PNT  For entering a point - (refer to the 'Disc' section for explanation)

' LINE  For entering a Line - (refer to the 'Disc' section for explanation)

' GSUR  For entering a surface - (refer to the 'Disc' section for explanation)

' ELTY  To assign a specific element type to a set of entities and generating the mesh - (refer to the 'Disc' section for explanation)

=Homework 5=

Given
The strong form of GM1.0b/d1b on page 2 of Mtg 9 shown below

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 $$,$$ \forall x \in \left[ 0,1 \right]$$, and $$\Gamma_g = [0]$$, is $$u(0)=4$$, $$\Gamma_h = [1]$$, is $$ -\frac{\partial u(x=1)}{\partial x}=6 $$

Solution
Using weighting function $$ w(x) $$ such as $$ w(0) = 0 $$ and the technique of integration by parts, we have derived the following weak form.

Fourier basis function
The trial solution then becomes the following for n=1 to convergence:

The essential boundary condition,$$u(x=0)=4$$ must be satisfied, and the basis functions must be selected so $$ w(n=0) = 0 $$. Solving the above basis function for u(x=1)=4 yielded the following: $$ {d}_{1}=4, i=1,2,\cdots $$, $$ j=1,2,\cdots $$

The remaining displacement coefficients must be solved for $$d_{j}, \;\;j=1,2,\cdots $$, by constructing the following stiffness and force matrices using the switching Fourier basis function.

Matlab was used to generate the following stiffness and force matrices shown below.

Since Kd=F, d can be solved for using Matlab to perform the matrix algebra. The resulting d matrix is shown below.

The trial solution of the governing equation then becomes the following.

Now we plot for Exact solution u(x) and Approximate solution uh(x) and also Convergence of error vs. number of degrees of freedom (n)

As seen from the Figure 5.1.1 below, it is evident that the approximate solution uh(x) very closely imitates the exact solution u(x) Also it can be seen from Figure 5.1.2. that the error converges to almost zero value when more than 8 basis functions are used

Find
Plot three figures similar to those of the linear Lagrangian element basis functions..

Solution
The generic equation for the quadratic Lagrangian element basis functions (QLEBF) is presented below in Equation 5.6.1;

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

$$ L_{j,n}=\prod_{k=1}^{n}\frac{(x-x_{k})}{(x_{i}-x_{k})} $$
 * $$\displaystyle (Eq.5.6.1) $$
 * }
 * }

Equation 5.6.1 can be expanded to the following equation to fit each of our cases as shown in Equations 5.6.2-4 below with x1=1, x2=2, and x3=3.

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

$$ L_{1,3}=\frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})} $$
 * <p style="text-align:right;">$$\displaystyle (Eq.5.6.2) $$
 * }
 * }

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

$$ L_{2,3}=\frac{(x-x_{3})(x-x_{1})}{(x_{2}-x_{3})(x_{2}-x_{1})} $$
 * <p style="text-align:right;">$$\displaystyle (Eq.5.6.3) $$
 * }
 * }

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

$$ L_{3,3}=\frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})} $$
 * <p style="text-align:right;">$$\displaystyle (Eq.5.6.4) $$
 * }
 * }

Since the quadratic Lagrangian element basis functions have three nodes, three QLEBFs were used to generate the following plots 5.6.1-4 in MATLAB. The code used to generate these plots is presented below.

Comments
Each of the plotted quadratic Lagrangian element basis functions has three nodes as opposed to the two nodes of the linear basis functions due to the n=3 condition. It is observed from the plots that the sum of the basis functions at each node must be equal to 1 regardless of the number of nodes. This also holds true for the linear Lagrangian element basis functions as well.

=Homework 6=

Given
The G2Dm1.0/D1 PDE

The G2Dm1.0/D1 PDE can be reduced to where $$\left[ \right]$$ is an identity matrix.

The essential boundary condition is $$g = 2\ \ \ \ on\ \ \ \ \partial \Omega $$.

Find
1) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) and verify results with any FE codes with detailed documentation.

2) Plot resulting solution in 3d with contours.

Solution
1) For Eq.6.5.1, we can write the weak form

And then we choose Lagrange Interpolation Basis Function for each element with n=m=2

 * $$N_{I}(x,y)=L_{i,m}(x)\cdot L_{j,n}(y),\;\; where \;\;I=i+(j-1)m$$

After substituting Eq.6.5.4 into Eq.6.5.3 in which Lagrange Interpolation Basis Function are both trial function and weight function, we integrate over the whole domain to get stiffness matrix and achieve final results using the following Matlab codes

The resulting stiffness and force matrices for the n=m=2 case is as follows:

$$K_{EF}=\begin{bmatrix} 0.6667 & -0.1667 & 0 & -0.1667 & 0 & 0 & 0 & 0 & -0.3333\\ -0.1667 & 1.3333 & -0.1667 & -0.3333 & -0.3333 & 0 & 0 & 0 & -0.3333\\ 0 & -0.1667 & 0.6667 & 0 & -0.1667 & 0 & 0 & 0 & -0.3333\\ -0.1667 & -0.3333 & 0 & 1.3333 & 0 & -0.1667 & -0.3333 & 0 & -0.3333\\ 0 & -0.3333 & -0.1667 & 0 & 1.3333 & 0 & -0.3333 & -0.16667 & -0.3333\\ 0 & 0 & 0 & -0.1667 & 0 & 0.6667 & -0.16667 & 0 & -0.3333\\ 0& 0 & 0 & -0.3333 & -0.3333 & -0.1667 & 1.3333 & -0.16667 & -0.3333\\ 0 & 0 & 0 & 0 & -0.1667 & 0 & -0.1667 & 0.6667 & -0.3333\\ -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & 2.6667 \end{bmatrix}$$

$$F=\begin{bmatrix} 0.1110x10^{-15}& -0.2220x10^{-15}& 0.1110x10^{-15}& 0x10^{-15}& -0.1110x10^{-15}& 0.1110x10^{-15}& 0.1110x10^{-15}& 0x10^{-15}& \end{bmatrix}^{T}$$

Knowing that Kd=F and solving for the d matrix yielded:

 * $$d=

\begin{bmatrix} 2& 2& 2& 2& 2& 2& 2& 2& 2& \end{bmatrix}^{T}$$

The result from the MATLAB code is as follows:

The following is the matlab code used to generate the following plots as well as calculate the above matrices.

2)

Based on the above results we see that regardless of how many elements we use, we get a uniform temperature distribution. The exact solution is easily obtained, because this is static problem with uniform boundary conditions. The 3D contour plots are flat surfaces with no contour due to the uniform temperature distribution all the way around the perimeter of the plate. If the temperature varied along one of the edges or at a node, then the plots would have some contour and an uneven temperature distribution.

=My Team Members= Nakul Deshpande Shreyas Joglekar William Kurth Yongtao Li Anand Tupsakhare Chris Vork

title: User:Eml5526.s11.team6.gravois


text:
User:EGM6341.s11.TEAM1.HW2

Find
1) Do integration by parts on last term of (2.1) to reveal three more terms in Taylor series, i.e.

2) Use Intergral Mean Value Theorem to express remainder in terms of $$ f^{(5)}( \xi\ ) $$ for $$ \xi\ \in [x_0,x]  $$

3) Assume (1.2) and (1.3) are correct. Perform integration by parts once more to verify (1.2) and (1.3) for $$ (n+1) \ $$ expansion with $$ R_{(n+2)}(x) \ $$

4) Use Intergral Mean Value Theorem on (1.3) to show (1.4)

Solution
Required

Integrating by parts

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

\begin{align} \int U'V = UV - \int UV'\\ \int UV' = UV - \int U'V\\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 *  (1.1)
 * }
 * }

Part 1) 1st term Derive the first term down below from

Taking integration by parts from the integrating part of the LHS,

Back to the main equation including the factorial term,

2nd term Derive the first term down below from

Note: Integrating by parts

Taking integration by parts from the integrating part of the LHS,

Back to the main equation including the factorial term,

3rd term Derive the first term down below from

Taking integration by parts from the integrating part of the LHS,

Back to the main equation including the factorial term,

In conclusion,

Part.2)

Consider the remainder as a function that we derived in HW1, prob3

Applying the remainder function to below equation,

Integrating the integration part of the RHS of the equation yields,

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

$$\displaystyle \begin{align} \therefore \frac{1}{4!} \int_{x_0}^{x} (x-t)^4 \cdot f^{(5)}(t) dx = \frac{(x-x_0)^5}{5!} f^{(5)}(\xi) \quad _{for} \quad w(x) \geqq 0 \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = | (1.17)
 * }
 * }

Part.3) Derive the first term down below from

With same method as used in part a), use integrating by parts again,

Back to the main equation including the factorial term,

Part.4) Consider the remainder as a function that we derived in HW1, prob3

Applying below IMVT to equation (1.22)

Integrating the integration part of the RHS of the equation yields,

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

$$\displaystyle \begin{align} \therefore \frac{1}{(n+1)!}\int_{x_0}^{x} (x-t)^{n+1} \cdot f^{(n+2)}(t) dx = \frac{(x-x_0)^{n+2}}{(n+2)!}f^{(n+2)}(\xi) \quad _{for} \quad w(x) \geqq 0 \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = | (1.25)
 * }
 * }

Author EGM6341.s11.TEAM1.Yoon 11:13, 26 January 2011 (UTC) - Primary Author

Find
1) Construct Taylor series around (2.2)

2) Plot these series for each n.

3) Find/Estimate max $$ \left | R_{n+1}(x=\frac{3 \pi\ }{4}) \right \vert $$

Solution
1) Construct Taylor series

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

$$\displaystyle \begin{align} P_0(x)=\sin(\frac{3\pi}{8})=0.4730 \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_1(x)=f(\frac{3\pi}{8}) + \frac{(x-x_0)^{(1)}}{1!} f^{(1)}(x_0) $$

$$\begin{align}       =0.9239 + (x-\frac{3\pi}{8})(\cos(\frac{3\pi}{8})) \end{align}$$
 * }

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

$$\displaystyle \begin{align} P_1(x)=0.4730 + 0.3827x \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_2(x) = P_1(x) + \frac{(x-x_0)^{(2)}}{2!} f^{(2)}(x_0) $$

$$\displaystyle        = (0.3827x + 0.4730) + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8}))$$
 * }

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

$$\displaystyle \begin{align} P_2(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_3(x) = P_2(x) + \frac{(x-x_0)^{3}}{3!} f^{(3)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8}))  $$
 * }

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

$$\displaystyle \begin{align} P_3(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_4(x) = P_3(x) + \frac{(x-x_0)^{4}}{4!} f^{(4)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (\sin(\frac{3\pi}{8})) $$
 * }

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

$$\displaystyle \begin{align} P_4(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (0.9239) \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_5(x) = P_4(x) + \frac{(x-x_0)^{5}}{5!} f^{(5)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (\cos(\frac{3\pi}{8})) $$
 * }

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

$$\displaystyle \begin{align} P_5(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (0.9239) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (0.3827) \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_6(x) = P_5(x) + \frac{(x-x_0)^{6}}{6!} f^{(6)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (\cos(\frac{3\pi}{8})) $$ $$\displaystyle + \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-\sin(\frac{3\pi}{8})) $$
 * }

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

$$\displaystyle \begin{align} P_6(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (0.9239) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (0.3827) \\+ \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-0.9239)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_7(x) = P_6(x) + \frac{(x-x_0)^{7}}{7!} f^{(7)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (\cos(\frac{3\pi}{8})) $$ $$\displaystyle + \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-\cos(\frac{3\pi}{8})) $$
 * }

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

$$\displaystyle \begin{align} P_7(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (0.9239) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (0.3827) \\+ \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-0.9239) \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-0.3827)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_8(x) = P_7(x) + \frac{(x-x_0)^{8}}{8!} f^{(8)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (\cos(\frac{3\pi}{8})) $$ $$\displaystyle + \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{8}}{8!} (\sin(\frac{3\pi}{8})) $$
 * }

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

$$\displaystyle \begin{align} P_8(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (0.9239) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (0.3827) \\+ \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-0.9239) \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-0.3827) + \frac{(x-\frac{3\pi}{8})^{8}}{8!} (0.9239)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_9(x) = P_8(x) + \frac{(x-x_0)^{9}}{9!} f^{(9)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (\cos(\frac{3\pi}{8})) $$ $$\displaystyle + \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{8}}{8!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{9}}{9!} (\cos(\frac{3\pi}{8})) $$
 * }

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

$$\displaystyle \begin{align} P_9(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (0.9239) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (0.3827) \\+ \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-0.9239) \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-0.3827) + \frac{(x-\frac{3\pi}{8})^{8}}{8!} (0.9239) + \frac{(x-\frac{3\pi}{8})^{9}}{9!} (0.3827)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

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

 * $$\displaystyle P_{10}(x) = P_9(x) + \frac{(x-x_0)^{10}}{10!} f^{(10)}(x_0) $$

$$\displaystyle = 0.3827x + 0.4730 + \frac{(x-\frac{3\pi}{8})^{2}}{2!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (\cos(\frac{3\pi}{8})) $$ $$\displaystyle + \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{8}}{8!} (\sin(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{9}}{9!} (\cos(\frac{3\pi}{8})) + \frac{(x-\frac{3\pi}{8})^{10}}{10!} (-\sin(\frac{3\pi}{8})) $$
 * }

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

$$\displaystyle \begin{align} P_{10}(x) = 0.4730 + 0.3827x + \frac{(x-\frac{3\pi}{8})^{2}}{2!}(-0.9239) + \frac{(x-\frac{3\pi}{8})^{3}}{3!}(-0.3827) + \frac{(x-\frac{3\pi}{8})^{4}}{4!} (0.9239) + \frac{(x-\frac{3\pi}{8})^{5}}{5!} (0.3827) \\+ \frac{(x-\frac{3\pi}{8})^{6}}{6!} (-0.9239) \frac{(x-\frac{3\pi}{8})^{7}}{7!} (-0.3827) + \frac{(x-\frac{3\pi}{8})^{8}}{8!} (0.9239) + \frac{(x-\frac{3\pi}{9})^{9}}{9!} (0.3827) + \frac{(x-\frac{3\pi}{8})^{10}}{10!} (-0.9239)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

2) Plot Taylor series

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

 *  Matlab code:
 * }

3)Find/Estimate Max Remainder $$ \left | R_{n+1}(x=\frac{3 \pi\ }{4}) \right \vert $$

$$\begin{align}R_{n+1}(x) \end{align}$$ in eqn.(2.3) is defined as following;

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

 * $$\begin{align} max \bigg| R_1(x) \bigg|=\frac{x-x_0}{1!} = \frac{\frac{3\pi}{4}-\frac{3\pi}{8}}{1!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_1(x) \bigg| = 1.1781$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_2(x) \bigg|=\frac{{x-x_0}^2}{2!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^2}{2!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_2(x) \bigg| = 0.6940$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_3(x) \bigg|=\frac{{x-x_0}^3}{3!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^3}{3!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_3(x) \bigg| = 0.2725$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_4(x) \bigg|=\frac{{x-x_0}^4}{4!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^4}{4!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_4(x) \bigg| = 0.0803$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_5(x) \bigg|=\frac{{x-x_0}^5}{5!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^5}{5!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_5(x) \bigg| = 0.0189$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_6(x) \bigg|=\frac{{x-x_0}^6}{6!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^6}{6!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_6(x) \bigg| = 0.0037$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_7(x) \bigg|=\frac{{x-x_0}^7}{7!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^7}{7!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_7(x) \bigg| = 6.2494e-004$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_8(x) \bigg|=\frac{{x-x_0}^8}{8!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^8}{8!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_8(x) \bigg| = 9.2030e-005$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_9(x) \bigg|=\frac{{x-x_0}^9}{9!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^9}{9!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_9(x) \bigg| = 1.2047e-005$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_{10}(x) \bigg|=\frac{{x-x_0}^{10}}{10!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^{10}}{10!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_{10}(x) \bigg| = 1.4192e-006$$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

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

 * $$\begin{align} max \bigg| R_{11}(x) \bigg|=\frac{{x-x_0}^{11}}{11!} = \frac{(\frac{3\pi}{4}-\frac{3\pi}{8})^{11}}{11!}\end{align}$$
 * }

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

$$\displaystyle max \bigg| R_{11}(x) \bigg| = 1.5200e-007 $$
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style="width:35%; padding:10px; border:2px solid #8888aa" |
 * style= |
 * }
 * }

Author --EGM6341.s11.TEAM1.Sujin 00:03, 30 January 2011 (UTC)

Find
1) Expand $$ e^x \ $$ in Taylor Series up to order n, with remainder $$ R_{n+1}[e^x_ix] $$

2) Find Taylor Series of $$ f(x) \ $$ up to order n, with remainder $$ R_{n+1}[f_{i}x] \ $$ equal to:

Solution
1) Expand Taylor series of exponential function $$\begin{align}, e^x \end{align}$$
 * Taylor series expansion of $$\displaystyle f(x) $$ is defined as following, refer eqn.(2.4);

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

$$\begin{align} f_n(x) =\underbrace{P_n(x)}_{\text{Taylor Series}} +\underbrace{R_{n+1}(x)}_{\text{Remainder}} \end{align}$$ where,
 * style="width:95%" |
 * style="width:95%" |
 * }

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

$$\begin{align} P_n(x)= f(x_0) + \frac{{(x-x_0)}^{(n)}}{1!} f^{(1)}(x_0) + ... + \frac{(x-x_0)^{(n)}}{n!} f^{(n)}(x_0) \end{align}$$
 * style="width:95%" |
 * style="width:95%" |

(3.4)
 *
 * }
 * ,and
 * {| style="width:100%" border="0"

$$\begin{align} R_{n+1}(x)=\frac{1}{n!}\int_{x_{0}}^{x}(x-t)^{n}f^{(n+1)}(t)dt=\frac{(x-x_{0})^{n+1}}{(n+1)!}f^{(n+1)}(\xi), \quad \xi\in[x_{0},x] \end{align}$$
 * style="width:95%" |
 * style="width:95%" |

(3.5)
 *
 * }

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

By using eqn.(2.4) we can expand $$\begin{align} e^x \end{align}$$ in Taylor series around $$\begin{align} x_0=0 \end{align}$$ ;
 * style="width:95%" |
 * style="width:95%" |

Derivative of $$\begin{align}e^x\end{align}$$ is its own; $$\begin{align} \frac{d^{n}}{dx^{n}}(e^{x})=e^{x} \end{align}$$ So, by using above property of $$\begin{align}e^x\end{align}$$ and eqn.(3.4) $$\begin{align} P_n(x) \end{align}$$ can obtaine $$\begin{align} P_n(x)=e^{0}+\frac{(x-0)}{1!}e^{0}+\frac{(x-0)^{2}}{2!}e^{0}+\cdots+\frac{(x-0)^{n}}{n!}e^{0} \end{align}$$

$$\begin{align} =P_n(x) = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \cdots + \frac{x^n}{n!} = \sum_{i=0}^{n}\frac{x^i}{i!} \end{align}$$

(3.6)
 *
 * }

 * From eqn.(3.5), we can also find $$\begin{align} R_{n+1}(x) \end{align}$$.
 * {| style="width:100%" border="0"

$$\begin{align} R_{n+1}(x)=\frac{(x-0)^{n+1}}{(n+1)!}e^{\xi}=\frac{x^{n+1}}{(n+1)!}e^{\xi} \end{align}$$
 * style="width:95%" |
 * style="width:95%" |

(3.7)
 *
 * }
 * Plugging eqn.(3.6) and eqn.(3.7) into eqn(2.4), we can meet Taylor series expansion of $$\displaystyle e^x $$

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

$$\displaystyle \begin{align} e^{x}=P_{n}(x)+R_{n+1}(x)=\sum_{j=0}^{n}\frac{x^{j}}{j!}+\underbrace{\frac{x^{n+1}}{(n+1)!}e^{\xi}}_{R_{n+1}[e^x;x]} {}, \xi\epsilon [0,x] \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |

(3.8)
 *
 * }

2) Expand Taylor series of f(x) $$\begin{align}, f(x)=\frac{e^x-1}{x} \end{align}$$

 * Taylor expansion of $$\begin{align}, f(x)=\frac{e^x-1}{x} \end{align}$$ can be expressed as following;
 * {| style="width:100%" border="0"

$$\begin{align} f(x) = \frac{e^x-1}{x} = \frac{1}{x}(e^{x}-1)=\frac{1}{x}(\cancel{1}+\underbrace{\frac{x}{1!}+\frac{x^2}{2!}+\cdots+\frac{x^n}{n!}}_{P_n(x)}+\underbrace{\frac{x^{n+1}}{(n+1)!}e^{\xi}}_{R_{n+1}} { } \cancel{-1}) \end{align}$$
 * style="width:95%" |
 * style="width:95%" |

$$\begin{align} = \cancel{\frac{1}{x}} \sum_{i=1}^{n}\frac{x^ {\cancelto{i-1}{i}} }{i!}+\frac{x^ {\cancelto{n}{n+1}} }{\cancel{x}(n+1)!} e^{\xi} = \sum_{i=1}^{n}\frac{x^ {i-1} }{i!}+\frac{x^{n}}{(n+1)!} e^{\xi} \end{align}$$
 * }

 * Thus,

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

$$\displaystyle \begin{align} f(x) = \frac{e^x-1}{x}=P_{n}(x)+R_{n+1}(x)=\sum_{j=0}^{n}\frac{x^{j-1}}{j!}+\underbrace{\frac{x^{n}}{(n+1)!}e^{\xi}}_{R_{n+1}[f(x);x]} {}, \xi\epsilon [0,x] \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |

(3.9)
 *
 * }

 * Furthermore, we can confirm equation [[media:Egm6341.s11.mtg7.pdf|(2)p.7-1]] using eqn.(3.8) and eqn.(3.9);
 * {| style="width:100%" border="0"

From eqn.(3.8), we can find $$\begin{align}R_n[e^x;x]\end{align}$$; $$\begin{align} R_n[e^x;x]=\frac{(x-0)^n}{n!}e^{\xi} { },{ } \xi \epsilon [0,x] \end{align}$$ Also, we obtained $$\begin{align} R_{n+1}[f(x)=\frac{e^x-1}{x};x] \end{align}$$ in eqn.(3.9); $$\begin{align} R_{n+1}[f;x] = \frac{x^{n}}{(n+1)!}e^{\xi} { },{ } \xi \epsilon [0,x] \end{align}$$ $$\begin{align} f(x)-f_n(x)=R_{n+1}(x)=\frac{x^{n}}{(n+1)!}e^{\xi} = \underbrace{\frac{(x)^n}{n!}e^{\xi}}_{R_n[e^x;x]}\frac{1}{n+1} \end{align}$$ $$\begin{align} \Rightarrow R_{n+1}[f;x]=\frac{R_n[e^x;x]}{n+1} \end{align}$$
 * style="width:95%" |
 * style="width:95%" |
 * }

Author --EGM6341.s11.TEAM1.Sujin 00:03, 30 January 2011 (UTC)

Find
1) Plot f(x)

2) Integrate f (find $$ I_n \ $$) for n=2,4,8...until error of order $$ 10^6 \ $$,using the following methods:

2a) Taylor Series Expansion $$ f_n \ $$

2b) Computational Trapezoidal Rule

2c) Composite Simpson Rule

2d) Gauss Legendre Quadrature

Solution
Part 1

 Matlab code:

Part 2A (Taylor series)

From Homework 1.2 ,

To integrate,

Using the error definition given in Meeting 7-1 :

From Homework 1.2 ,

By the Integral Mean Value Theorem,

Finding the maximum and minimum of $$f(x)$$ on the interval and applying the Intermediate Value Theorem,

So the error will be bounded depending on $$n$$ in the following table

Using equation 4.5, the value for the integral for each respective $$n$$ are:

 Matlab code:

Part 2B (Trapezoidal Rule)

For the following parts, Wolfram Alpha was used to determine the "actual" value of $$I(x)$$ to be

The trapezoidal rule states:

Exercising this for a variety of $$n$$ and replacing the value of $$f(0)$$ whenever encountered with $$f(0)=1$$ (from HW1.1 ), the following table is generated:

 Matlab code:

Part 2C (Simpson's Rule)

Simpson's rule states:

Exercising this for a variety of $$n$$ and replacing the value of $$f(0)$$ whenever encountered with $$f(0)=1$$ (from HW1.1 ), the following table is generated:

 Matlab code:

Part 2D (Gauss-Legendre Quadrature)

The Gauss-Legendre Quadrature states:

where the values for $$x_i$$ and $$w_i$$ are the roots of the Legendre polynomial of order $$n$$ and their associated weights, respectively.

Exercising this for a variety of $$n$$ and replacing the value of $$f(0)$$ whenever encountered with $$f(0)=1$$ (from HW1.1 ), the following table is generated:

 Matlab code:

Author Egm6341.s11.team1.arm 02:43, 31 January 2011 (UTC)

Find
Verify given table against NIST Handbook.

Solution
The table checks with the values presented in the | NIST Handbook of Mathmatical Functions (NIST values are in decimal format).

Author Egm6341.s11.team1.arm 19:59, 2 February 2011 (UTC)

Given
Where $$ \in \Rho\ _n= $$ set of polynomials of degrees less than or equal to n

Find
Verify Eq6.3 thru Eq6.7 using Eq6.1 and Eq6.2

Solution
When n = 0:

Matlab generated Polynomials n=0 through 4

When n = 1:

Matlab generated Polynomials n=0 through 4

When n = 2:

Matlab generated Polynomials n=0 through 4

When n = 3:

Matlab generated Polynomials n=0 through 4

When n = 4:

Matlab generated Polynomials n=0 through 4

Author EGM6341.s11.TEAM1.Yoon 11:13, 26 January 2011 (UTC) - Primary Author

Find
1) Verify the family of Legendre Polynomials as shown in $$Eq. 7.1$$ is orthogonal for $$ n=0,1,...,5 \ $$

2) Construct Gram Matrix: $$ \underline{ \Gamma\ } $$ as shown in $$Eq. 7.2$$

3) Show Gram Matrix:$$ \underline{ \Gamma\ } $$ is diagonal and compute the determinate.

for $$ P_0, P_1,P_3 \ $$, do by hand.

for $$ P_4, P_5,P_6 \ $$, use WolframAlpha (http://www.wolframalpha.com/) and link in the solution.

HINT: If family of Legendre Polynomials is orthogonal, the Gram Matrix:$$ \underline{ \Gamma\ } $$ should be diagonal.

Solution
1) Here we verified for the general case of the family of Legendre Polynomials as shown in $$Eq. [7.1]$$ is orthogonal and in Part 2 and 3 we varified for the Particular cases of $$ n=0,1,...,5 \ $$

We start off by writing out the Legendre polynomials as follows

Next we want to multiply the above $$Eq. [7.2]$$ by $$P_m(x)$$ and integrating over the closed interval $$[-1,1].$$

We are going to want to use integration by parts where

We want to show that for each time we use integration by parts that the boundary term of each integration will be zero over the interval $$[-1,1].$$

Now let's interchange $$n$$ and $$m$$.

If we take $$Eq. 7.7$$ and subtract it from $$Eq. 7.6$$ we will obtain an expression which tells us something about the orthogonality of the Legendre Polynomials.

Therefore,for the general result of $$m\ne n$$

For $$ m = n \,\!$$

where $$N_n(x)$$ is a normalization constant. We can obtain an expression for the normalization constant by considering the Rodriguez formula for the Legendre Polynomial.

Again, as we did before in $$Eq. [7.4],$$ we are going to integrate by parts, but this time we are going to repeat the integration $$ n $$ times.

To take the nth derivative of the product in the integrand above in $$ Eq. [7.14]$$ we use the general Leibnitz Rule for the nth derivative of a product.

After applying the above Leibnitz Rule to $$ Eq. [7.14]$$ we come to the following evaluation of the 2n-fold derivative of the product $$(x^2-1)^n.$$

allowing us to express the normalization constant $$N_n(x)$$ as the following:

Now we can rewrite the integrand in $$ Eq. [7.15].$$

Plugging $$ Eq. [7.16]$$ into $$Eq. [7.15]$$ gives us:

Note that:

Recalling $$ Eq. [7.15],$$ we can now express the $$2^{nd}$$ term of $$Eq. [7.20]$$ as:

Using $$ Eq. [7.21]$$ we can rewrite $$ Eq. [7.20]$$ as:

Rearranging shows:

$$Eq. [7.23]$$ shows that the value of $$(2n+1)N_n(x)$$ is a constant and independent of $$n.$$ This is due to the fact that both the R.H.S. and L.H.S. of $$Eq. [7.23]$$ have the same dependence on $$n$$, therefore showing that overall, the R.H.S. with respect to the L.H.S. of $$Eq. [7.23]$$ is independent of $$n$$.This allows us to evaluate $$(2n+1)N_n(x)$$ at $$n=0$$:

Dividing $$Eq. [7.24]$$ through by $$2n+1$$ gives us the sought after result:

2) Construct Gram Matrix: $$ \underline{ \Gamma\ } $$ as shown in $$Eq. [7.2]$$

Let's calculate the inner products for the particular cases $$n=0,1,...,5.$$

We will integrate the rest of the inner products via Wolframalpha to not only obtain the remaining terms of the Gram Matrix, but as well test our predictions in Part 1 for the general solution of the inner products of orthogonal functions.

Note that the prediction for the normalization constant of orthogonal functions in Part 1 of this problem in $$Eq. [7.25]$$ was proven correct with $$Eq.[7.33], Eq.[7.39], Eq.[7.40], Eq.[7.44-45].$$ Now that we have obtained all the components of the Gram Matrix, we can plug in the values and is that it is indeed a diagonal matrix.

3) Show Gram Matrix:$$ \underline{ \Gamma\ } $$ is diagonal and compute the determinate.

We can already see from above in Part 2 that the Gram Matrix is a diagonal matrix. The determinant is easy computed for diagonal matrices because it's nothing more than the product of the pivot terms (i.e. diagonal terms).

Author Egm6341.s11.team-1.langpm 21:46, 31 January 2011 (UTC)

Find
Derive Eq8.2 from Eq8.1

Solution
Use Newton-Cotes Method. Integration of

and

Set $$x0 = a$$, $$x1 = b$$, $$x$$ in $$[a,b]$$, we get

Where

So,

which equals 1 for $$x=x_0$$ and 0 for $$x=x_1$$, and

which equals 0 for $$x=x_0$$ and 1 for $$x=x_1$$.

Author EGM6341.s11.team1.Chiu 13:50, 2 February 2011 (UTC)

Find
Use $$Eq. [9.1]$$ to generate Legendre Polynomial $$ P_5(x) \ $$.

Use MATLAB command 'roots" to compare the roots of $$ P_5(x) \ $$ to check values in table below.

Plot the roots on $$\left[ -1,1\right]$$ using MATLAB "plot" command. Plot dots "." with coordinates $$ \left ( x_i,y_i \right ), i=1,...,5 $$

Use roots of $$ P_5(x) = x_i \ $$ and $$ y_i=0 \ $$, markersize=14.

Repeat above steps for $$ P_{10}(x) \ $$. Observe the location of the roots near end points $$\left[-1\right]$$ and $$\left[1\right].$$

Solution
Simplifying the coefficients we come to the finally expression for the $$5th$$ order Legendre Polynomial:

Author Egm6341.s11.team-1.langpm 20:37, 31 January 2011 (UTC)

Find
Use Eq10-1 and Eq10-2 to find an expression for $$ \left \{ c_i \right \} $$ in terms of $$ \left ( x_i,f(x_i), i=0,1,2 \right )$$

Solution
Substitute for $$ i=0,1,2 \ $$, into Eqs 10.1 and 10.2 to obtain the following three equations with three unknowns:

Solve for $$ c_0, c_1, \ $$ and $$ c_2 \ $$ to obtain expressions of $$ \left \{ c_i \right \} $$ in terms of $$ \left ( x_i,f(x_i), i=0,1,2 \right )$$

Author

Find
Use Eq11-1 to derive Simple Simpson's Rule

Solution
$$p_2 \ $$ must be integrated from $$x_0 \ $$ to $$x_2 \ $$. To simplify, $$x_0 \ $$ and $$x_2 \ $$ can be set to -1 and 1, respectively, without loss of generality.

Author Egm6341.s11.team1.arm 03:21, 1 February 2011 (UTC)

Find
Consider $$ n=1,2,4,8,16 \ $$.

Construct $$ f_n(x) \ $$ as shown in Eq12.2.

Plot $$ f,f_n,n=1,2,4,8,16 \ $$

Compute $$ I_n= \int_{a}^{b} f_n(x)\, dx $$ for $$ n=1,2,4,8 \ $$  and compare to $$ I \ $$ (use Wolfram Alpha with more digits)

Plot $$ l_0,l_1,l_2 \ $$ for $$ n=5 \ $$

How would $$ l_3,l_4 \ $$ look?

Solution
1) Construct $$ f_n(x) \ $$ as shown in Eq12.2. (a) $$\displaystyle n=1 $$

 * $$\displaystyle x_0=-1,\quad x_1=1 $$

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

\begin{align} f_1(x) &= P_1(x)\\ &=\sum_{i=0}^{1} l_{i,1}(x)f(x_i)\\ &=l_{0,1}(x)f(x_0)+l_{1,1}f(x_1) \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }

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

\begin{align} l_{0,1}(x)&=\frac{(x-x_1)}{(x_0-x_1)}=\frac{x-1}{-2}, \quad f(x_0)=f(-1)=1-e^{-1} \\ l_{1,1}(x)&=\frac{(x-x_0)}{(x_1-x_0)}=\frac{x+1}{2}, \quad f(x_1)=f(1)=e-1 \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }

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

$$\displaystyle \begin{align} f_1(x)&=p_1(x)\\ &=(\frac{x-1}{-2})(1-e^{-1}) + (\frac{x+1}{2})(e-1)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = | (12.3)
 * }
 * }

(b) $$\displaystyle n=2 $$
 * $$\displaystyle x_0=-1,\quad x_1=0,\quad x_2=1 $$
 * {| style="width:100%" border="0" align="left"

\begin{align} f_2(x) &= P_2(x)\\ &=\sum_{i=0}^{2} l_{i,2}(x)f(x_i)\\ &=l_{0,2}(x)f(x_0)+l_{1,2}f(x_1)+l_{2,2}f(x_2)\\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * Where
 * $$\displaystyle

\begin{align} l_{0,2}(x)&=\frac{(x-x_1)}{(x_0-x_1)}\frac{(x-x_2)}{(x_0-x_2)}=\frac{x(x-1)}{2}, &&f(x_0)=f(-1)=1-e^{-1} \\ l_{1,2}(x)&=\frac{(x-x_0)}{(x_1-x_0)}\frac{(x-x_2)}{(x_1-x_2)}=\frac{(x+1)(x-1)}{-1}, && f(x_1)=f(0)=1 \\ l_{2,2}(x)&=\frac{(x-x_0)}{(x_2-x_0)}\frac{(x-x_1)}{(x_2-x_1)}=\frac{(x+1)x}{2}, && f(x_1)=f(1)=e+1 \end{align} $$
 * Thus,

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

$$\displaystyle \begin{align} f_2(x)&=P_2(x)\\ &=( \frac{x(x-1)}{2} )(1-e^{-1}) + (\frac{(x+1)(x-1)}{-1})(1) + (\frac{(x+1)x}{2})(e+1)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = | (12.4)
 * }
 * }

(c)$$\displaystyle n=4 $$
 * $$\displaystyle

x_0=-1,\quad x_1=-\frac{1}{2}, \quad x_2=0, \quad x_3=\frac{1}{2}, \quad x_4=1 $$

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

\begin{align} f_4(x)&=P_4(x)\\ &=\sum_{i=0}^{4}l_{i,4}(x)f(x_{i})\\ &=l_{0,4}(x)f(x_{0})+l_{1,4}(x)f(x_{1})+l_{2,4}(x)f(x_{2})+l_{3,4}(x)f(x_{3})+l_{4,4}(x)f(x_{4}) \end{align}$$
 * $$\displaystyle
 * $$\displaystyle
 * }

 * Where,
 * $$\displaystyle

\begin{align} l_{0,4}&=\frac{(x+\frac{1}{2})(x-0)(x-\frac{1}{2})(x-1)}{(-1+\frac{1}{2})(-1-0)(-1-\frac{1}{2})(-1-1)}, &&f(x_{0})=f(-1)=1-e^{-1} \\ l_{1,4}&=\frac{(x+1)(x-0)(x-\frac{1}{2})(x-1)}{(-\frac{1}{2}+1)(-\frac{1}{2}-0)(-\frac{1}{2}-\frac{1}{2})(-\frac{1}{2}-1)}, &&f(x_{1})=f(-\frac{1}{2})=\frac{e^{-\frac{1}{2}}-1}{-\frac{1}{2}} \\ l_{2,4}&=\frac{(x+1)(x+\frac{1}{2})(x-\frac{1}{2})(x-1)}{(0+1)(0+\frac{1}{2})(0-\frac{1}{2})(0-1)}, &&f(x_{2})=f(0)=1 \\ l_{3,4}&=\frac{(x+1)(x+\frac{1}{2})(x-0)(x-1)}{(\frac{1}{2}+1)(\frac{1}{2}+\frac{1}{2})(\frac{1}{2}-0)(-\frac{1}{2}-1)}, &&f(x_{3})=f(\frac{1}{2})=\frac{e^{\frac{1}{2}}-1}{\frac{1}{2}} \\ l_{4,4}&=\frac{(x+1)(x+\frac{1}{2})(x-0)(x-\frac{1}{2})}{(-1+1)(-1+\frac{1}{2})(-1-0)(-1-\frac{1}{2})}, &&f(x_{4})=f(1)=e^{1}-1 \\ \end{align} $$

 * Thus,

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

$$\displaystyle \begin{align} f_4(x) =& P_4(x)\\ =&(\frac{(x+\frac{1}{2})(x-0)(x-\frac{1}{2})(x-1)}{(-1+\frac{1}{2})(-1-0)(-1-\frac{1}{2})(-1-1)}) (1-e^{-1}) +(\frac{(x+1)(x-0)(x-\frac{1}{2})(x-1)}{(-\frac{1}{2}+1)(-\frac{1}{2}-0)(-\frac{1}{2}-\frac{1}{2})(-\frac{1}{2}-1)})(\frac{e^{-\frac{1}{2}}-1}{-\frac{1}{2}})\\ &+(\frac{(x+1)(x+\frac{1}{2})(x-\frac{1}{2})(x-1)}{(0+1)(0+\frac{1}{2})(0-\frac{1}{2})(0-1)})(1) +(\frac{(x+1)(x+\frac{1}{2})(x-0)(x-1)}{(\frac{1}{2}+1)(\frac{1}{2}+\frac{1}{2})(\frac{1}{2}-0)(-\frac{1}{2}-1)})(\frac{e^{\frac{1}{2}}-1}{\frac{1}{2}})\\ &+(\frac{(x+1)(x+\frac{1}{2})(x-0)(x-\frac{1}{2})}{(-1+1)(-1+\frac{1}{2})(-1-0)(-1-\frac{1}{2})})(e^{1}-1)\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = | (12.5)
 * }
 * }

(d) $$\displaystyle n=8 $$
 * $$\displaystyle

x_0=-1,\quad x_1=-\frac{3}{4},\quad x_2=-\frac{1}{2},\quad x_3=-\frac{1}{4},\quad x_4=0,\quad x_5=\frac{1}{4},\quad x_6=\frac{1}{2},\quad x_7=\frac{3}{4}, \quad x_8=1 $$
 * {| style="width:100%" border="0" align="left"

\begin{align} f_8(x)=P_8(x)=\sum_{i=0}^{8}l_{i,8}(x)f(x_{i})&=l_{0,8}(x)f(x_{0})+l_{1,8}(x)f(x_{1})+l_{2,8}(x)f(x_{2})+l_{3,8}(x)f(x_{3})\\ &+l_{4,8}(x)f(x_{4})+l_{5,8}(x)f(x_{5})+l_{6,8}(x)f(x_{6}) + l_{7,8}(x)f(x_{7})+l_{8,8}(x)f(x_{8}) \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }

 * where,
 * $$\displaystyle

\begin{align} l_{0,8}&=\frac{(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-1+\frac{3}{4})(-1+\frac{1}{2})(-1+\frac{1}{4})(-1-0)(-1-\frac{1}{4})(-1-\frac{1}{2})(-1-\frac{3}{4})(-1-1)}, \quad f(x_{0})=f(-1)=1-e^{-1} \\ l_{1,8}&=\frac{(x+1)(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-\frac{3}{4}+1)(-\frac{3}{4}+\frac{1}{2})(-\frac{3}{4}+\frac{1}{4})(-\frac{3}{4}-0)(-\frac{3}{4}-\frac{1}{4})(-\frac{3}{4}-\frac{1}{2})(-\frac{3}{4}-\frac{3}{4})(-\frac{3}{4}-1)}, \quad f(x_{1})=f(-\frac{3}{4})==\frac{e^{-\frac{3}{4}}-1}{-\frac{3}{4}} \\ l_{2,8}&=\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-\frac{1}{2}+1)(-\frac{1}{2}+\frac{3}{4})(-\frac{1}{2}+\frac{1}{4})(-\frac{1}{2}-0)(-\frac{1}{2}-\frac{1}{4})(-\frac{1}{2}-\frac{1}{2})(-\frac{1}{2}-\frac{3}{4})(-\frac{1}{2}-1)}, \quad f(x_{2})=f(-\frac{1}{2})=\frac{e^{-\frac{1}{2}}-1}{-\frac{1}{2}} \\ l_{3,8}&=\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-\frac{1}{4}+1)(-\frac{1}{4}+\frac{3}{4})(-\frac{1}{4}+\frac{1}{2})(-\frac{1}{4}-0)(-\frac{1}{4}-\frac{1}{4})(-\frac{1}{4}-\frac{1}{2})(-\frac{1}{4}-\frac{3}{4})(-\frac{1}{4}-1)}, \quad f(x_{2})=f(-\frac{1}{4})=\frac{e^{-\frac{1}{4}}-1}{-\frac{1}{4}} \\ l_{4,8}&=\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{3}{4})(x-1)}{(0+1)(0+\frac{3}{4})(0+\frac{1}{2})(0+\frac{1}{4})(0-\frac{1}{4})(0-\frac{1}{2})(0-\frac{3}{4})(0-1)}, \quad f(x_{4})=f(0)=1 \\ l_{5,8}&=\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{2})(x-\frac{3}{4})(x-1)}{(\frac{1}{4}+1)(\frac{1}{4}+\frac{3}{4})(\frac{1}{4}+\frac{1}{2})(\frac{1}{4}+\frac{1}{4})(\frac{1}{4}-0)(\frac{1}{4}-\frac{1}{2})(\frac{1}{4}-\frac{3}{4})(\frac{1}{4}-1)}, \quad f(x_{5})=f(\frac{1}{4})=\frac{e^{\frac{1}{4}}-1}{\frac{1}{4}} \\ l_{6,8}&=\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(\frac{1}{2}+1)(\frac{1}{2}+\frac{3}{4})(\frac{1}{2}+\frac{1}{2})(\frac{1}{2}+\frac{1}{4})(\frac{1}{2}-0)(\frac{1}{2}-\frac{1}{4})(\frac{1}{2}-\frac{3}{4})(\frac{1}{2}-1)}, \quad f(x_{6})=f(\frac{1}{2})=\frac{e^{\frac{1}{2}}-1}{\frac{1}{2}} \\ l_{7,8}&=\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{2})(x-\frac{1}{4})(x-1)}{(\frac{3}{4}+1)(\frac{3}{4}+\frac{3}{4})(\frac{3}{4}+\frac{1}{2})(\frac{3}{4}+\frac{1}{4})(\frac{3}{4}-0)(\frac{3}{4}-\frac{1}{4})(\frac{3}{4}-\frac{1}{2})(\frac{3}{4}-1)}, \quad f(x_{7})=f(\frac{3}{4})=\frac{e^{\frac{3}{4}}-1}{\frac{3}{4}} \\ l_{8,8}&=\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})}{(1+1)(1+\frac{3}{4})(1+\frac{1}{2})(1+\frac{1}{4})(1-0)(1-\frac{1}{4})(1-\frac{1}{2})(1-\frac{3}{4})}, \quad f(x_{8})=f(1)=\frac{e-1}{1} \\ \end{align} $$

 * Thus,

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

$$\displaystyle \begin{align} f_8(x)=&(\frac{(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-1+\frac{3}{4})(-1+\frac{1}{2})(-1+\frac{1}{4})(-1-0)(-1-\frac{1}{4})(-1-\frac{1}{2})(-1-\frac{3}{4})(-1-1)}) (1-e^{-1}) \\ &+(\frac{(x+1)(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-\frac{3}{4}+1)(-\frac{3}{4}+\frac{1}{2})(-\frac{3}{4}+\frac{1}{4})(-\frac{3}{4}-0)(-\frac{3}{4}-\frac{1}{4})(-\frac{3}{4}-\frac{1}{2})(-\frac{3}{4}-\frac{3}{4})(-\frac{3}{4}-1)}) (\frac{e^{-\frac{3}{4}}-1}{-\frac{3}{4}})\\ &+(\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-\frac{1}{2}+1)(-\frac{1}{2}+\frac{3}{4})(-\frac{1}{2}+\frac{1}{4})(-\frac{1}{2}-0)(-\frac{1}{2}-\frac{1}{4})(-\frac{1}{2}-\frac{1}{2})(-\frac{1}{2}-\frac{3}{4})(-\frac{1}{2}-1)}) (\frac{e^{-\frac{1}{2}}-1}{-\frac{1}{2}} )\\ &+(\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x-0)(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(-\frac{1}{4}+1)(-\frac{1}{4}+\frac{3}{4})(-\frac{1}{4}+\frac{1}{2})(-\frac{1}{4}-0)(-\frac{1}{4}-\frac{1}{4})(-\frac{1}{4}-\frac{1}{2})(-\frac{1}{4}-\frac{3}{4})(-\frac{1}{4}-1)}) (\frac{e^{-\frac{1}{4}}-1}{-\frac{1}{4}})\\ &+(\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-\frac{1}{4})(x-\frac{1}{2})(x-\frac{3}{4})(x-1)}{(0+1)(0+\frac{3}{4})(0+\frac{1}{2})(0+\frac{1}{4})(0-\frac{1}{4})(0-\frac{1}{2})(0-\frac{3}{4})(0-1)}) (1)\\ &+(\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{2})(x-\frac{3}{4})(x-1)}{(\frac{1}{4}+1)(\frac{1}{4}+\frac{3}{4})(\frac{1}{4}+\frac{1}{2})(\frac{1}{4}+\frac{1}{4})(\frac{1}{4}-0)(\frac{1}{4}-\frac{1}{2})(\frac{1}{4}-\frac{3}{4})(\frac{1}{4}-1)}) (\frac{e^{\frac{1}{4}}-1}{\frac{1}{4}})\\ &+(\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{4})(x-\frac{3}{4})(x-1)}{(\frac{1}{2}+1)(\frac{1}{2}+\frac{3}{4})(\frac{1}{2}+\frac{1}{2})(\frac{1}{2}+\frac{1}{4})(\frac{1}{2}-0)(\frac{1}{2}-\frac{1}{4})(\frac{1}{2}-\frac{3}{4})(\frac{1}{2}-1)}) (\frac{e^{\frac{1}{2}}-1}{\frac{1}{2}})\\ &+(\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{2})(x-\frac{1}{4})(x-1)}{(\frac{3}{4}+1)(\frac{3}{4}+\frac{3}{4})(\frac{3}{4}+\frac{1}{2})(\frac{3}{4}+\frac{1}{4})(\frac{3}{4}-0)(\frac{3}{4}-\frac{1}{4})(\frac{3}{4}-\frac{1}{2})(\frac{3}{4}-1)}) (\frac{e^{\frac{3}{4}}-1}{\frac{3}{4}})\\ &+(\frac{(x+1)(x+\frac{3}{4})(x+\frac{1}{2})(x+\frac{1}{4})(x-0)(x-\frac{1}{2})(x-\frac{1}{4})(x-\frac{3}{4})}{(1+1)(1+\frac{3}{4})(1+\frac{1}{2})(1+\frac{1}{4})(1-0)(1-\frac{1}{4})(1-\frac{1}{2})(1-\frac{3}{4})}) (\frac{e-1}{1})\\ \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = | (12.6)
 * }
 * }

(e) $$\displaystyle n=16 $$
 * $$\displaystyle

\begin{align} x_0&=-1,\quad x_1=-\frac{7}{8},\quad x_2=-\frac{3}{4},\quad x_3=-\frac{5}{8}, \cdots ,\quad x_{15}=\frac{7}{8},\quad x_{16}=1 \end{align} $$
 * {| style="width:100%" border="0" align="left"

\begin{align} f_{16}(x)=P_{16}(x)=\sum_{i=0}^{16}l_{i,16}(x)f(x_{i})&=l_{0,16}(x)f(x_{0}) + l_{1,16}(x)f(x_{1}) + l_{2,16}(x)f(x_{2})+l_{3,16}(x)f(x_{3})\\ &+ l_{4,16}(x)f(x_{4}) + l_{5,16}(x)f(x_{5}) + l_{6,16}(x)f(x_{6}) + l_{7,16}(x)f(x_{7}) + l_{8,16}(x)f(x_{8})+l_{9,16}(x)f(x_{9})\\ &+l_{10,16}(x)f(x_{10}) + l_{11,16}(x)f(x_{11}) + l_{12,16}(x)f(x_{12}) + l_{13,16}(x)f(x_{13})+l_{14,16}(x)f(x_{14})\\ &+l_{15,16}(x)f(x_{15}) + l_{16,16}(x)f(x_{16}) \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }

 * $$\displaystyle

\begin{align} l_{0,16}&=\prod_{j=0}^{16}\frac{x-x_{j}}{x_{0}-x_{j}}, \quad f(x_{0})=f(-1)=\frac{e^{-1}-1}{-1} \\ l_{1,16}&=\prod_{j=0}^{16}\frac{x-x_{j}}{x_{1}-x_{j}}, \quad f(x_{1})=f(-\frac{7}{8})=\frac{e^{-\frac{7}{8}}-1}{-\frac{7}{8}}\\ & \vdots \\ l_{16,16}&=\prod_{j=0}^{16}\frac{x-x_{j}}{x_{16}-x_{j}}, \quad f(x_{16})=f(1)=\frac{e-1}{1} \end{align} $$

Thus,

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

$$\displaystyle \begin{align} f_{16}(x)=&(\prod_{j=0}^{16}\frac{x-x_{j}}{x_{0}-x_{j}})(\frac{e^{-1}-1}{-1}) + (\prod_{j=0}^{16}\frac{x-x_{j}}{x_{1}-x_{j}}) (\frac{e^{-\frac{7}{8}}-1}{-\frac{7}{8}}) + \cdots + (\prod_{j=0}^{16}\frac{x-x_{j}}{x_{16}-x_{j}}) (\frac{e-1}{1}) \end{align} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = | (12.7)
 * }
 * }

2) Compute $$ I_n= \int_{a}^{b} f_n(x)\, dx $$ for $$ n=1,2,4,8 $$

I = 2.1145wolframAlpha
 * {|class = "wikitable collapsible"

!colspan="3"| 3) Plot $$ l_0,l_1,l_2 \ $$ for $$ n=5 \ $$
 * n
 * In
 * En
 * 1
 * 2.35040238
 * -0.23590238
 * 2
 * 2.11680663
 * -0.00230663
 * 4
 * 2.11451758
 * -0.00001758
 * 8
 * 2.11450788
 * -0.00000788
 * }
 * 8
 * 2.11450788
 * -0.00000788
 * }
 * }

Matlab generated LIET plot for n={1,2,4,8,16}

Matlab generated Lagrange-Interpolation for n={1,2,4,8,16}

When n=1 When n=2 When n=4 When n=8 When n=16

Co-Author --EGM6341.s11.TEAM1.Sujin 20:03, 2 February 2011 (UTC)

Find
Show that Eq13.2 can be derived from Eq13.1.

Show that Eq13.4 can be derived from Eq13.3.

Solution
Trapezoidal Rule

The trapezoidal rule states:

When $$b-a \ $$ is not sufficiently small, the trapezoidal rule is not useful. Then, we break it into a sum of integrals over small subintervals.

Set $$n \geq 1$$, $$h=\frac{b-a}{n}$$, and $$x_j=a+jh \ $$ for $$j=0,1, \ldots ,n$$ with $$x_j-1 \leq \eta_j \leq x_j$$

We get:

The first terms in the sum can be combined to give the composite trapezoidal rule:

Simpson's Rule

Simpson's rule states:

Set $$n \geq 2$$ and even, define $$h=\frac{b-a}{n}$$, $$x_j = a+jh \ $$ for $$j=0,1,\ldots,n$$ with $$x_{2j-2} \leq \eta_j \leq x_{2j}$$,

Hence, we get the composite Simpson's rule by using the same method as the above,

Author EGM6341.s11.team1.Chiu 14:51, 2 February 2011 (UTC)

Find
Give details to obtain Eq14.1

Solution
Knowing that:

From Lecture 2-2, Eq1:

From Atkinson, pg270, Eq (5.3.1)

Also from Atkinson, pg270, Eq (5.3.2)

By combining terms in Eq14.3, Eq14.4 and Eq14.5 and given the definition of Eq14.2, it can be shown that:

AuthorEGM6341.s11.TEAM1.WILKS 06:39, 30 January 2011 (UTC)EGM6341.s11.TEAM1.WILKS

Signatures
Solved problem 1 -- EGM6341.s11.TEAM1.Yoon 11:13, 26 January 2011 (UTC) - Primary Author

Solved problem 2 -- EGM6341.s11.TEAM1.Sujin 00:05, 30 January 2011 (UTC)

Solved problem 3 -- EGM6341.s11.TEAM1.Sujin 00:05, 30 January 2011 (UTC)

Solved problem 4 -- Egm6341.s11.team1.arm 02:43, 31 January 2011 (UTC)

Solved problem 5 -- Egm6341.s11.team1.arm 19:59, 2 February 2011 (UTC)

Solved problem 6 -- EGM6341.s11.TEAM1.Yoon 11:13, 26 January 2011 (UTC) - Primary Author

Solved problem 7 -- Egm6341.s11.team-1.langpm 21:46, 31 January 2011 (UTC)

Solved problem 8 -- EGM6341.s11.team1.Chiu 13:50, 2 February 2011 (UTC)

Solved problem 9 -- Egm6341.s11.team-1.langpm 21:46, 31 January 2011 (UTC)

Solved problem 10 -- EGM6341.s11.TEAM1.WILKS 18:46, 29 January 2011 (UTC)EGM6341.s11.TEAM1.WILKS

Solved problem 11 -- Egm6341.s11.team1.arm 03:21, 1 February 2011 (UTC)

Solved problem 12 -- EGM6341.s11.TEAM1.Yoon 11:13, 31 January 2011 (UTC) - Co-Author, EGM6341.s11.TEAM1.Sujin 00:05, 30 January 2011 (UTC) - Co-Author

Solved problem 13 -- EGM6341.s11.team1.Chiu 14:51, 2 February 2011 (UTC)

Solved problem 14 -- EGM6341.s11.TEAM1.WILKS 06:40, 30 January 2011 (UTC)EGM6341.s11.TEAM1.WILKS

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

title: User:EGM6341.s11.TEAM1.HW2


text:
Fluid Mechanics for MAP/Analytical solutions of internal and external flows

>back to Chapters of Fluid Mechanics for MAP

>back to Chapters of Microfluid Mechanics

Concept of Fully Developed Flow

Consider the flow in a channel between two plates having a height of $$D$$ and an infinite depth in $$x_{3}$$ direction. Starting from the entrance, the boundary layers develop due to the no-slip condition on the wall. At a finite distance, the boundary layers merge and the inviscid core (field with no velocity gradient in $$x_{2}$$ direction) vanishes. The flow becomes fully viscous. The velocity field in $$x_{1}$$ direction adjusts slightly further until $$x_{1}=L_e$$ and it no longer changes with $$x_{1}$$ direction. This state of the flow is called fully-developed. At that state:

$$\frac{\partial U_{i}}{\partial x_{1}} = 0 \rightarrow \frac{\partial \tau_{ij}}{\partial x_{1}} = 0$$

Because of the two-dimensional nature of the flow, no gradient of the velocity quantities in $$x_{3}$$ direction is expected starting from the entrance.

$$U_{3}= 0\, \ \frac{\partial U_{i}}{\partial x_{3}}= 0\ , \ \frac{\partial \tau_{ij}}{\partial x_{3}}= 0\ , \ \frac{\partial P}{\partial x_{3}}= 0,$$ $$U_{2}= U_{3}= 0$$

Hence,

$$U_{1} = U_{1}(x_{2})\ ,\ \tau_{ij}=\tau_{ij}(x_2).$$

The entrance lengths for laminar pipe and channel flows are, respectively :

$$	\frac{Le_{pipe}}{D} = \left[(0.619)^{1.6} + (0.0567Re)^{1.6}\right]^{\frac{1}{1.6}} $$

$$ \frac{Le_{channel}}{D} = \left[(0.631)^{1.6} + (0.0442Re)^{1.6}\right]^{\frac{1}{1.6}} $$

Fully Developed Laminar Flow Between Infinite Parallel Plates

Consider the fully developed laminar flow between two infinite plates.

Consider the continuity equation and momentum equation in $$x_{1}$$ direction for an incompressible steady flow between two infinite plates as shown.

Continuity Equation

$$ \underbrace{\frac{\partial \rho}{\partial t}}_{= 0} + \frac{\partial \rho U_{i}}{\partial x_{i}} = 0\ \stackrel{\rho = cst.}{\longrightarrow}\ \frac{\partial U_{i}}{\partial x_{i}} = 0 $$

Since $$\frac{\partial U_{1}}{\partial x_{1}} = 0$$, $$\frac{\partial U_{3}}{\partial x_{3}} = 0 \rightarrow \frac{\partial U_{2}}{\partial x_{2}} = 0$$ because it is a fully-developed and two dimensional flow. Hence, $$\ U_2$$ reads $$\ U_2=constant $$ As $$\ U_2$$ is zero on the walls, it should be zero in the whole fully developed region, i.e. $$\ U_2(x_2)=0 $$

Momentum Equation in j-direction

$$ \underbrace{\rho\frac{\partial U_{j}}{\partial t}}_{= 0\ Steady} + \rho U_{i}\frac{\partial U_{j}}{\partial x_{i}} = -\frac{\partial P}{\partial x_{j}} + \rho g_{j} + \frac{\partial \tau_{ij}}{\partial x_{i}} $$

in $$x_{1}$$ direction ,$$g_{1} = 0$$

$$\underbrace{\rho U_{i}\frac{\partial U_{1}}{\partial x_{i}}}_{A} = -\frac{\partial P}{\partial x_{1}} + \underbrace{\frac{\partial \tau_{i1}}{\partial x_{i}}}_{B}$$

Consider term A:

$$\rho U_{i}\frac{\partial U_{1}}{\partial x_{i}} = \rho \left[U_{1}\underbrace{\frac{\partial U_{1}}{\partial x_{1}}}_{= 0\ Fully-developed} + \underbrace{U_{2}}_{= 0}\frac{\partial U_{1}}{\partial x_{2}} + \underbrace{U_{3}}_{=0\ fully-developed\ 2D }\underbrace{\frac{\partial U_{1}}{\partial x_{3}}}_{= 0\ 2D}\right]=0$$

Consider term B:

$$\frac{\partial \tau_{i1}}{\partial x_{i}} = \underbrace{\frac{\partial \tau_{11}}{\partial x_{1}}}_{= 0\ fully-dev.} + \frac{\partial \tau_{21}}{\partial x_{2}} + \underbrace{\frac{\partial \tau_{31}}{\partial x_{3}}}_{= 0\ 2D} = \frac{\partial \tau_{\tau_{21}}}{\partial x_{2}} = \frac{d\tau_{21}}{dx_{2}}$$

hence $$\tau_{21} = \tau_{21}(x_{2})$$.

Thus the momentum equation in $$x_{1}$$ direction reads:

$$0 = -\frac{\partial P}{\partial x_{1}} + \frac{d\tau_{21}}{dx_{2}}$$

This equation should be valid for all $$x_{1}$$ and $$x_{2}$$. This requires that $$\frac{\partial P}{\partial x_{1}} = \frac{d \tau_{21}}{d x_{2}}$$ = constant.

Remember $$\tau_{21}$$ is the stress in $$x_{1}$$ direction on a face normal to $$x_{2}$$ direction.

$$ \tau_{ij} = \mu \left(\frac{\partial U_{i}}{\partial x_{j}} + \frac{\partial U_{j}}{\partial x_{i}}\right) + \frac{2}{3}\delta_{ij}\mu\frac{\partial U_{k}}{\partial x_{k}} $$

thus,

$$\tau_{21} = \mu\left(\underbrace{\frac{\partial U_{2}}{\partial x_{1}}}_{= 0} + \frac{\partial U_{1}}{\partial x_{2}}\right) + 0$$

$$\tau_{21} = \mu\frac{\partial U_{1}}{\partial x_{2}}\ \text{since,}\ U_{1} = U_{1}(x_{2})\ \tau_{21} = \mu\frac{dU_{1}}{dx_{2}}$$

Thus the momentum equation reads:

$$ \frac{\partial P}{\partial x_{1}} = \mu\frac{d^{2}U_{1}}{dx_{2}^{2}} $$

This equation can be obtained also by using the Reynold's transport equations for a differential volume.

The momentum equation in $$x_{1}$$ direction,

$$ F_{S1} + \underbrace{F_{B1}}_{=0} = \frac{\partial}{\partial t} \underbrace{\int_{CV}\rho U_{1}dV}_{=0} + \underbrace{\int_{CS}U_{1}\rho U_{i}n_{i}dA}_{=0\ fully-dev.} $$

The flux term becomes zero since for fully-developed flow incoming flux is equal to the outgoing flux. Thus,

$$\displaystyle F_{S1} = 0$$

That is:

$$ \left(p - \frac{\partial P}{\partial x_{1}}\frac{dx_{1}}{2}\right)dx_{2}dx_{3} - \left(p + \frac{\partial P}{\partial x_{1}}\frac{dx_{1}}{2}\right)dx_{2}dx_{3}$$

$$ + \left(\tau_{21} + \frac{dx_{2}}{2}\frac{\partial \tau_{21}}{\partial x_{2}}\right)dx_{1}dx_{3} - \left(\tau_{21} - \frac{dx_{2}}{2}\frac{\partial \tau_{21}}{\partial x_{2}}\right)dx_{1}dx_{3} = 0 $$

$$-\frac{\partial P}{\partial x_{1}}dV + \frac{\partial \tau_{21}}{\partial x_{2}}dV = 0$$

$$-\frac{\partial P}{\partial x_{1}} + \frac{\partial \tau_{21}}{\partial x_{2}} = 0$$

Finally, the governing equation of this kind of flow becomes:

$$ \frac{\partial P}{\partial x_{1}} = \mu\frac{d^{2}U_{1}}{dx_{2}^{2}} $$

with the following boundary conditions:

$$x_{2} = 0 \rightarrow U_{1} = 0$$ and $$x_{2} = D \rightarrow U_{1} = 0$$

Integrating the equation once results in a linear function of $$x_{2}$$:

$$ \mu\frac{dU_{1}}{dx_{2}} = \left(\frac{\partial P}{\partial x_{1}}\right)x_{2} + c_{1} \rightarrow \tau_{21} = \left(\frac{\partial P}{\partial x_{1}}\right)x_{2} + c_{1} $$

The second integration reads:

$$ U_{1} = \frac{1}{2\mu}\left(\frac{\partial P}{\partial x_{1}}\right)x_{2}^{2} + \frac{1}{\mu} c_{1}x_{2} + c_{2} $$

The integration constants is obtained by using the boundary conditions:

$$x_{2} = 0,\ U_{1} = 0 \rightarrow c_{2} = 0$$

$$x_{2} = D,\ U_{1} = 0 \rightarrow c_{1} = -\frac{1}{2}\left(\frac{\partial P}{\partial x_{1}}\right)D$$

Finally, the velocity profile reads:

$$ \begin{array}{lll} U_{1} &=& \frac{1}{2\mu}\left(\frac{\partial P}{\partial x_{1}}x_{2}^{2}\right) - \frac{1}{2\mu}\left(\frac{\partial P}{\partial x_{1}}\right)D{x_{2}}\\ &=&\frac{D^{2}}{2\mu}\left(\frac{\partial P}{\partial x_1}\right) \left[\left(\frac{x_{2}}{D}\right)^{2} - \left(\frac{x_{2}}{D}\right)\right] \end{array} $$

Note that the velocity profile is parabolic!

The shear stress becomes:

$$ \begin{array}{lll} \tau_{21} &=& \left(\frac{\partial P}{\partial x_{1}}\right)x_{2} - \frac{1}{2}\left(\frac{\partial P}{\partial x_{1}}\right)D \\ &=& D\left(\frac{\partial P}{\partial x_{1}}\right)\left[\frac{x_{2}}{D} - \frac{1}{2}\right] \end{array} $$

at the wall i.e. at $$x_{2}$$ = 0 and $$x_{2}$$= D

$$ \tau_{21}(0) = -\frac{1}{2}D\left(\frac{\partial P}{\partial x_{1}}\right) $$

$$ \tau_{21}(D) = \frac{1}{2}D\left(\frac{\partial P}{\partial x_{1}}\right) $$

Note that $$\tau_{21}$$ is maximum near the wall, i.e. momentum loss is maximum near the wall. This is due to the maximum velocity gradient $$\frac{\partial U_{1}}{\partial x_{2}}$$ near the wall!

The volume flow rate is,

$$ Q = \int_{A}U_{i}n_{i}dA = \int_{0}^{D}U_{1}w dx_{2} $$

where $$w$$ is the depth of the channel.

Thus the volume flow rate per depth $$w$$ is given by:

$$ \frac{Q}{w} = \int_{0}^{D}\frac{1}{2\mu}\left(\frac{\partial P}{\partial x_{1}}\right)\left(x^{2} - D{x_{2}}\right)dx_{2} $$

$$ \frac{Q}{w} = -\frac{1}{12\mu}\left(\frac{\partial P}{\partial x_{1}}\right)D^{3} $$

Note that $$\frac{\partial P}{\partial x_{1}}$$ should be constant for the fully developed flow. Hence, for a channel with a finite length $$L$$:

$$\frac{\partial P}{\partial x_{1}} = \frac{p_{2} - p_{1}}{L} = \frac{-\Delta P}{L}$$

Where $$\Delta P$$ is the pressure drop along L.

$$\frac{Q}{w} = -\frac{1}{12\mu}\left(\frac{-\Delta P}{L}\right)D^{3} = \frac{D^{3}}{12\mu L}\Delta P$$

or the pressure drop can be calculated from:

$$ \Delta P = \frac{Q}{l}\frac{12\mu L}{D^{3}} $$

For the same flow rate, increasing the height of the channel would cause a drastic reduction in the pressure drop.

The average velocity $$\bar{U}$$ is:

$$ \bar{U} = U_{0} = \frac{Q}{A} = \frac{Q}{w D} = -\frac{1}{12\mu}\left(\frac{\partial P}{\partial x_{1}}\right)\frac{D^{3}w}{w D} = -\frac{1}{12\mu}\left(\frac{\partial P}{\partial x_{1}}\right)D^{2} $$

The maximum velocity occurs when:

$$ \frac{dU_{1}}{dx_{2}} = 0 = \frac{D^{2}}{2\mu}\left(\frac{\partial P}{\partial x_{1}}\right)\left(\frac{2x_{2}}{D^{2}} - \frac{1}{D}\right) $$

Hence, at $$x_{2} = \frac{D}{2}$$, $$U_{1} = U_{1max}$$

$$U_{1}\left(\frac{D}{2}\right) = U_{1max} = -\frac{1}{8\mu}\left(\frac{\partial P}{\partial x_{1}}\right)D^{2} = \frac{3}{2}\bar{U}$$

The velocity profile can be written as functions of bulk velocity $$\overline{U}$$ or maximum velocity $$U_{1max}$$ by replacing their value the velocity profile equation:

$$ \begin{array}{lll} U_{1} &=& -4 U_{1max} \left[\left(\frac{x_{2}}{D}\right)^{2} - \left(\frac{x_{2}}{D}\right)\right]\\ &=&-6 \overline{U} \left[\left(\frac{x_{2}}{D}\right)^{2} - \left(\frac{x_{2}}{D}\right)\right] \end{array} $$

Same problem can be solved by using moving plates.

Example

Consider the hydrolic control valve coprising a piston, fitted to a cylinder with a mean radial clearance of 0,005mm. Determine the leakage flow rate. The fluid is SAE low oil ($$\rho$$ = 932 $$\frac{kg}{m^{3}}$$, $$\mu$$ =0.018 $$\frac{kg}{m\ sec}$$ at 55ºC). The flow can be assumed to be laminar, steady, incompressible, fully-developed flow. $$\left(\frac{L}{a} = 3000\right)$$

Since $$\frac{D}{a} = \frac{25}{0,005}$$ = 5000 the flow in the clearance can be accepted to be 2-D, with the depth $$w=\pi\cdot D$$, thus:

$$\frac{Q}{w} = \frac{a^{3}\Delta P}{12\mu L}$$

$$Q = \frac{a^{3}\Delta P}{12\mu L}\pi D = \frac{(0.005\cdot 10^{-3})^{3}}{12*(0.018)*15\cdot10^{-3}}*(20-1)\cdot10^{-6}* \pi * 25\cdot10^{-3}$$

$$Q = 57.6\cdot10^{-9} \frac{m^{3}}{s} = 57.6 \frac{mm^{3}}{sec}$$

Check the Reynolds number to ensure that laminar flow assumption is correct.

$$\bar{U} = \frac{Q}{A} = \frac{Q}{\pi D\ a} = 0.147 \frac{m}{sec}$$

$$Re = \frac{\rho \bar{U}a}{\mu} = \frac{932*0.147*0.005\cdot10^{-3}}{0.018} = 0.0375$$

Re $$<< 1800$$, i. e. the flow is laminar.

Layered Channel Flow

This channel flow contains two different and non miscible fluids. Fluids A and B flow at the same time through a channel, which is bounded two flat plates. They both occupy the half height of the channel. The fluid A has a viscosity $$\displaystyle \mu_{A}$$, a density $$\displaystyle \rho_{A}$$ and the mass flow $$\displaystyle \dot{m}_{A}$$. Fluid B, which is located above fluid A, has a viscosity $$\displaystyle \mu_{B}$$, a density $$\displaystyle \rho_{B}$$ and the mass flow $$\displaystyle \dot{m}_{B}$$. The following differential equations correspond to the molecular momentum $$\displaystyle \tau_{21}$$ for each Fluid.

$$ \frac{d\tau_{21}^{A}}{dx_{2}} = \frac{d\Pi}{dx_{1}}$$   and    $$\frac{d\tau_{21}^{B}}{dx_{2}} = \frac{d\Pi}{dx_{1}} $$.

With $$\tau_{21}= \mu \frac{dU_{1}}{dx_{2}}$$ yields the velocity field:

$$\frac{d^{2}U_{1}^{A}}{dx_{2}} = \frac{1}{\mu_{A}}\frac{d\Pi}{dx_{1}} $$   and  $$  \frac{d^{2}U_{1}^{B}}{dx_{2}} = \frac{1}{\mu_{B}}\frac{d\Pi}{dx_{1}}$$

After integration of both equations we obtain:

$$\tau_{21}^{A} = \frac{d\Pi}{dx_{1}}x_{2} + C_{1}^{A}$$ and $$\tau_{21}^{B} = \frac{d\Pi}{dx_{1}}x_{2} + C_{1}^{B}$$

As boundary condition we consider that shear stress on the interface between A and B is the same. Therefore we obtain:

$$\tau_{21}^{A}(x_{2} = 0) = \tau_{21}^{B}(x_{2} = 0)$$

Then,

$$C^{A}_{1} = C_{1}^{B} = C_{1}$$

After the integration for the velocity field:

$$ U_{1}^{A} = \frac{1}{2\mu_{A}}\frac{d\Pi}{dx_{1}}x_{2}^{2} + \frac{C_{1}}{\mu_{A}}x_{2} + C_{2}^{A}$$

and

$$ U_{1}^{B} = \frac{1}{2\mu_{B}}\frac{d\Pi}{dx_{1}}x_{2}^{2} + \frac{C_{1}}{\mu_{B}}x_{2} + C_{2}^{B}$$

The second boundary condition turns out to be on the interface:

$$ U_{1}^{A}(x_{2} = 0) = U_{1}^{B}(x_{2} = 0)$$, i.e. $$C_{2}^{A} \neq C_{2}^{B}$$

therefore,

$$C_{2}^{A} = C_{2}^{B} = C_{2}$$. The integration constants can be calculated with the following boundary conditions:

At $$ x_{2} = -D \ \rightarrow \ \ U_{1}^{A} = 0$$:

$$ 0 = \frac{d\Pi}{dx_{1}}\frac{1}{2\mu_{A}}D^{2} - \frac{C_{1}D}{\mu_{A}} + C_{2}$$

At $$ x_{2} = +D\ \rightarrow \ \ U_{1}^{B} = 0$$:

$$ 0 = \frac{d\Pi}{dx_{1}}\frac{1}{2\mu_{B}}D^{2} + \frac{C_{1}D}{\mu_{B}} + C_{2}$$

Therefore we obtain for the velocity distribution in the fluids A and B:

$$U_{1}^{A} = -\frac{D^{2}}{2\mu_{A}}\frac{d\Pi}{dx_{1}} \left[+\frac{2\mu_{A}}{(\mu_{A} + \mu_{B})} + \left(\frac{\mu_{A} - \mu_{B}}{\mu_{A}+ \mu_{B}}\right)\left(\frac{x_{2}}{D}\right) - \left(\frac{x_{2}}{D}\right)^{2}\right] \ \ \ \ \ $$

and

$$U_{1}^{B} = -\frac{D^{2}}{2\mu_{B}}\frac{d\Pi}{dx_{1}} \left[+\frac{2\mu_{B}}{(\mu_{A} + \mu_{B})} + \left(\frac{\mu_{A} - \mu_{B}}{\mu_{A}+ \mu_{B}}\right)\left(\frac{x_{2}}{D}\right) - \left(\frac{x_{2}}{D}\right)^{2}\right]$$

For the distribution of the shear stress we get:

$$\tau_{21} = D\frac{d\Pi}{dx_{1}}\left[\left(\frac{x_{2}}{D}\right) - \frac{1}{2} \left(\frac{\mu_{A} - \mu_{B}}{\mu_{A}+ \mu_{B}}\right)\right]$$

If we choose $$\displaystyle\mu_{A} = \mu_{B}$$,

$$U_{1} = \frac{-D^{2}}{2\mu_{A}}\frac{d\Pi}{dx_{1}}\left[1 - \left(\frac{x_{2}}{D}\right)^{2}\right]$$

$$\tau_{21} = D\frac{d\Pi}{dx_{1}}\left(\frac{x_{2}}{D}\right)$$

The solution gives that of the channel flow. In other words, velocity has a parabolic profile with the peak in the middle of the channel and a linear shear stress distribution $$\displaystyle\tau_{21}$$, where $$\displaystyle\tau_{21} = 0$$ at the channel's centerline.

If $$\mu_{A} \neq \mu_{B}$$, the position where the maximal velocity occurs can be calculated by introducing $$\displaystyle\tau_{21} = 0$$ on the velocity profile equation:

$$ x_2(U_{1 max}) = \frac{D}{2}\left(\frac{\mu_{A} - \mu_{B}}{\mu_{A}+ \mu_{B}}\right) $$

The shear stress on the upper plate is:

$$ \tau_{W}^{B} = \frac{D}{2}\frac{d\Pi}{dx_{1}}\left[\frac{\mu_{A} + 3\mu_{B}}{\mu_{A} + \mu_{B}}\right]$$

and the shear stress on the lower plate reads:

$$ \tau_{W}^{A} = -\frac{D}{2}\frac{d\Pi}{dx_{1}}\left[\frac{3\mu_{A} + \mu_{B}}{\mu_{A} + \mu_{B}}\right]$$

The average velocities of the fluids A and B are:

$$\tilde{U}_{1}^{A} = -\frac{D^{2}}{12\mu_{A}}\frac{d\Pi}{dx_{1}}\left(\frac{7\mu_{A} + \mu_{B}}{\mu_{A} + \mu_{B}}\right)$$

and

$$\tilde{U}_{1}^{B} = -\frac{D^{2}}{12\mu_{B}}\frac{d\Pi}{dx_{1}}\left(\frac{\mu_{A} + 7\mu_{B}}{\mu_{A} + \mu_{B}}\right)$$

Hence the respectively mass flow rates are:

$$\dot{m}_{A} = BD\tilde{U}_{1}^{A}$$

and

$$\dot{m}_{B} = BD\tilde{U}_{1}^{B}$$

Navier Stokes Equation in Cyclindrical Coordinates
A change of variables on the Cartesian equations will yield the following equations of momentum in r, $$\phi$$, and z directions:

r:\;\;\rho \left(\frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + \frac{u_{\phi}}{r} \frac{\partial u_r}{\partial \phi} + u_x \frac{\partial u_r}{\partial x} - \frac{u_{\phi}^2}{r}\right) = -\frac{\partial p}{\partial r} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_r}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_r}{\partial \phi^2} + \frac{\partial^2 u_r}{\partial x^2}-\frac{u_r}{r^2}-\frac{2}{r^2}\frac{\partial u_\phi}{\partial \phi} \right] + \rho g_r$$

\phi:\;\;\rho \left(\frac{\partial u_{\phi}}{\partial t} + u_r \frac{\partial u_{\phi}}{\partial r} + \frac{u_{\phi}}{r} \frac{\partial u_{\phi}}{\partial \phi} + u_x \frac{\partial u_{\phi}}{\partial x} + \frac{u_r u_{\phi}}{r}\right) = -\frac{1}{r}\frac{\partial p}{\partial \phi} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_{\phi}}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_{\phi}}{\partial \phi^2} + \frac{\partial^2 u_{\phi}}{\partial z^2} + \frac{2}{r^2}\frac{\partial u_r}{\partial \phi} - \frac{u_{\phi}}{r^2}\right] + \rho g_{\phi}$$

x:\;\;\rho \left(\frac{\partial u_x}{\partial t} + u_r \frac{\partial u_x}{\partial r} + \frac{u_{\phi}}{r} \frac{\partial u_x}{\partial \phi} + u_x \frac{\partial u_x}{\partial x}\right) = -\frac{\partial p}{\partial x} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_x}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_x}{\partial \phi^2} + \frac{\partial^2 u_x}{\partial x^2}\right] + \rho g_x.$$

The continuity equation is:

\frac{\partial \rho}{\partial t}+\frac{1}{r}\frac{\partial}{\partial r}\left(\rho r u_r\right) + \frac{1}{r}\frac{\partial }{\partial \phi}\left(\rho u_\phi\right) + \frac{\partial}{\partial x} \left(\rho u_x\right) = 0.$$

Fully Developed Pipe Flow

It is possible to use the same mathematical treatment like before to find and understand the velocity profile for fully developed flow inside a pipe with diameter D and infinite length. To show the flexibility, the same solution for this problem will be approached via 3 different ways.

(i) Infinitesimal Cylinder at the center of the pipe

Applying the RTT to the infinitesimal cylindrical CV along the symmetry axis of horizontal pipe, in which the flow is fully developed, the conservation of mass and the transport side of the conservation of momentum equation drops. Only remaining term governing this kind of flow is the balance of the forces on the CV in $$x$$ direction.

$$\left(p - \frac{\partial P}{\partial x}\frac{dx}{2}\right)\pi r^{2} - \left(p + \frac{\partial P}{\partial x}\frac{dx}{2}\right)\pi r^{2} + \tau_{r_{x}}2\pi rdx = 0$$

$$-\frac{\partial P}{\partial x}dx \pi r^{2} + \tau_{r_{x}}2\pi rdx = 0$$

$$ \tau_{r_{x}} = \frac{r}{2}\frac{\partial P}{\partial x} $$

(ii) Infinitesimal thin hollow cylinder at the center

This time the pressure and viscous force is considered for a concentric hollow cylinder with radius of R and infinitesimally small thickness dr and length dx(as shown in the image besides) along x-direction. Considering pressure term on the cross-section of cylinder

$$p2\pi rdr -\left(p + \frac{\partial P}{\partial x}{dx}\right)2\pi rdr\rightarrow -\left(\frac{\partial P}{\partial x}\right) 2\pi rdrdx$$

considering viscous shear stress on the surface of the cylinder $$\left(\tau_{r_{x}}+\frac{\partial \tau_{r_{x}}}{\partial r}dr\right)2\pi(r+dr)dx - \tau_{r_{x}}2\pi rdx$$ $$=\tau_{r_{x}} 2\pi rdx + \tau_{r_{x}} 2\pi drdx + \left(\frac{\partial \tau_{r_{x}}}{\partial r}\right)2\pi rdrdx + \underbrace{\left(\frac{\partial \tau_{r_{x}}}{\partial r}\right)2\pi r{dr}^2dx}_{{dr}^2\approx 0} - \tau_{r_{x}} 2\pi rdx $$ $$ \rightarrow \tau_{r_{x}} 2\pi drdx + \left(\frac{\partial \tau_{r_{x}}}{\partial r}\right)2\pi rdrdx $$

combining both term,we get the balance equation, $$-\left(\frac{\partial P}{\partial x}\right) 2\pi rdrdx + \tau_{r_{x}} 2\pi drdx + \left(\frac{\partial \tau_{r_{x}}}{\partial r}\right)2\pi rdrdx = 0 $$

it could be rewritten as $$-\left(\frac{\partial P}{\partial x}\right) 2\pi rdrdx +(\tau_{r_{x}} + r \frac{\partial \tau_{r_{x}}}{\partial r})2\pi drdx = 0 $$ dividing both side with $$2\pi drdx $$,we get $$-\left(\frac{\partial P}{\partial x}\right)r + \frac{\partial (r \tau_{r_{x}})}{\partial r} = 0$$ or,$$\left(\frac{\partial P}{\partial x}\right)r = \frac{\partial (r \tau_{r_{x}})}{\partial r} $$ integrating both side ,

$$\int \left(\frac{\partial P}{\partial x}\right)r\, dr = \int \partial (r \tau_{r_{x}}) $$

or,$$\left(\frac{\partial P}{\partial x}\right)\frac{r^2}{2}=r \tau_{r_{x}}$$

or,$$\left(\frac{\partial P}{\partial x}\right)\frac{r}{2}= \tau_{r_{x}}$$

However,in a laminar flow!

$$\mu\frac{dU}{dr} = \frac{1}{2}r\frac{\partial P}{\partial x}$$

integrating,

$$U = \frac{r^{2}}{4\mu}\left(\frac{\partial P}{\partial x}\right) + c_{1}$$

The boundary condition is:

$$U = 0\ at\ r = R$$

Thus $$c_{1}$$ can be calculated from the boundary condition.

$$ c_{1} = -\frac{R^{2}}{4\mu}\left(\frac{\partial P}{\partial x}\right) \Rightarrow U = \frac{1}{4\mu}\left(\frac{\partial P}{\partial x}\right)(r^{2} - R^{2}) $$

or

$$ U = -\frac{R^{2}}{4\mu}\left(\frac{\partial P}{\partial x}\right)\left[1 - \left(\frac{r}{R}\right)^{2}\right] $$

(iii) Using NS in cylindrical co-ordinates:

$$ x:\rho \left(\underbrace{\frac{\partial u_x}{\partial t}}_{steady state} + \underbrace{u_r \frac{\partial u_x}{\partial r}}_{u_r=0} + \underbrace{\frac{u_{\phi}}{r} \frac{\partial u_x}{\partial \phi}}_{u_{\phi}=0}+\underbrace {u_x \frac{\partial u_x}{\partial x}}_{\frac{\partial u_x}{\partial x} =0}\right) = -\frac{\partial p}{\partial x} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_x}{\partial r}\right) + \underbrace{\frac{1}{r^2}\frac{\partial^2 u_x}{\partial \phi^2}}_{=0} + \underbrace{\frac{\partial^2 u_x}{\partial x^2}}_{=0}\right] + \underbrace{\rho g_x}_{g_x=0}$$

thus,  $$\frac{\partial p}{\partial x} = \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_x}{\partial r}\right)\right]$$

or,$$\frac{\partial p}{\partial x} = \frac{\mu}{r}\left[\frac{\partial}{\partial r}\left(r \frac{\partial u_x}{\partial r}\right)\right]$$

or,$$\left(\frac{\partial p}{\partial x}\right)\frac{r}{\mu} = \left[\frac{\partial}{\partial r}\left(r \frac{\partial u_x}{\partial r}\right)\right]$$

Now, integrating both side with respect to $$r$$, we get

$$\int \left(\frac{\partial p}{\partial x}\right)\frac{r}{\mu} dr = \int \partial \left(r \frac{\partial u_x}{\partial r}\right)$$ then, $$ \left(\frac{\partial p}{\partial x}\right)\frac{r^2}{2\mu}+ C = \left(r \frac{\partial u_x}{\partial r}\right)$$ dividing by r in both side, we get then $$ \left(\frac{\partial p}{\partial x}\right)\frac{r}{2\mu}+ \frac{C}{r} = \left(\frac{\partial u_x}{\partial r}\right)$$ integrating again with respect to $$r$$ gives , $$ \left(\frac{\partial p}{\partial x}\right)\frac{r^2}{4\mu}+ C \ln r + D   = U_x(r)$$

Consequently, when $$r = 0$$ then $$U_x =U_{max}$$ and $$\ln r=\infty $$, as a result C=0.

On the other hand, when $$r = R$$ then $$U_x = 0$$, so $$\left(\frac{\partial p}{\partial x}\right)\frac{R^2}{4\mu}+ D = 0 $$ or, $$D = -\left(\frac{\partial p}{\partial x}\right)\frac{R^2}{4\mu}$$

By putting D to the primitive equation, we get, $$ U = -\frac{R^{2}}{4\mu}\left(\frac{\partial P}{\partial x}\right)\left[1 - \left(\frac{r}{R}\right)^{2}\right] $$

Knowing the velocity profile we can evaluate relevant quantities. The shear stress profile will look like:

$$ \tau_{r_{x}} = \mu\frac{dU}{dr} = \frac{r}{2}\left(\frac{\partial P}{\partial x}\right) $$

at r = 0 $$\rightarrow \tau_{r_{x}} = 0$$

at r = R $$\rightarrow \tau_{r_{x}} = \frac{R}{2}\left(\frac{\partial P}{\partial x}\right)$$

The volume flow rate would read

$$\displaystyle Q = \int_{Area} U_{i}n_{i}dA = \int_{Area} U2\pi rdr$$

$$Q = \int{\frac{1}{4\mu}}\left(\frac{\partial P}{\partial x}\right)(r^{2} - R^{2})2\pi rdr$$

$$Q = -\frac{\pi R^{4}}{8\mu}\left(\frac{\partial P}{\partial x}\right)$$

When we approximate $$\frac{\partial P}{\partial x} = -\frac{\Delta P}{L}$$

$$Q = -\frac{\pi R^{4}}{8\mu}\left[\frac{-\Delta P}{L}\right] = \frac{\pi \Delta PR^{4}}{8\mu L} = \frac{\pi \Delta PD^{4}}{128\mu L}$$

$$ \Delta P = \frac{128\mu L}{\pi D^{4}}Q $$

Increase radius to create drastic reduction in the pressure drop.

The mean velocity is:

$$ \bar{U} = \frac{Q}{A} = \frac{Q}{\pi R^{2}} = \frac{-R^{2}}{8\mu}\left(\frac{\partial P}{\partial x}\right) = -\frac{D^{2}}{32\mu}\left(\frac{\partial P}{\partial x}\right) $$

The location where maximum velocity occurs can be found be setting:

$$\frac{dU}{dr} = 0 \rightarrow \frac{dU}{dr} = \frac{1}{2\mu}\left(\frac{\partial P}{\partial x}\right)r$$

at r = 0 $$\rightarrow$$ U = $$U_{max}$$.

$$ U_{max} = U(0) = -\frac{R^{2}}{4\mu}\left(\frac{\partial P}{\partial x}\right) = 2\bar{U} $$

Note that in a channel was $$U_{max}=\frac{3}{2}\bar{U}$$.

$$U$$ can be written as a function of $$U_{max}$$ i.e.

$$ U = \underbrace{-\frac{R^{2}}{4\mu}\left(\frac{\partial P}{\partial x}\right)}_{U_{max}}\left[1 - \left(\frac{r}{R}\right)^{2}\right] $$

$$ \frac{U}{U_{max}} = \left[1 - \left(\frac{r}{R}\right)^{2}\right] $$

Again, the velocity profile becomes parabolic.


text: File:Internal-flow-examples.png

title: File:Internal-flow-examples.png