text:
User:Eml5526.s11.team7.flater/EGM6611.HW4

=Homework 4=

Description
A material is assumed to have the following constitutive equations:

 * Stress $$\displaystyle \Rightarrow \quad \mathbf{T} = \mathbf{T} (\mathbf{E},\mathbf{\dot{E}},\theta, \nabla \theta)$$
 * Heat flux $$\displaystyle \Rightarrow \quad \mathbf{q} = \mathbf{q} (\mathbf{E},\mathbf{\dot{E}},\theta, \nabla \theta)$$
 * Helmholtz free energy $$\displaystyle \Rightarrow \quad \psi = \psi (\mathbf{E},\mathbf{\dot{E}},\theta, \nabla \theta)$$
 * Entropy $$\displaystyle \Rightarrow \quad s = s (\mathbf{E},\mathbf{\dot{E}},\theta, \nabla \theta)$$

 * where $$\mathbf{E}$$ is Lagrangian finite strain tensor, $$\mathbf{\dot{E}}$$ is time derivative of the Lagrangian finite strain tensor, $$\displaystyle \theta$$ is temperature, and $$\displaystyle \nabla \theta$$ is the gradient of the temperature.

Derive the constitutive equations based on the Clausius-Duhem inequality.

Solution
Clausius-Duhem inequality is defined as

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

$$ \frac{ds}{dt} - \frac{r}{\theta} + \frac{1}{\rho \theta} \nabla \cdot \mathbf{q} - \frac{\mathbf{q}}{\rho \theta^2} \cdot \nabla \theta \ge 0 $$     (1.1)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

 * or in terms of the Helmholtz free energy

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

$$ \frac{\rho}{\theta}(\dot{\psi} +s \dot{\theta}) + \frac{1}{\theta} \mathbf{T}: \dot{\mathbf{E}} - \frac{1}{\theta^2}Q \cdot \rho \theta \ge 0 $$.     (1.2)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

Solve for $$\displaystyle \dot{\psi}$$ and substitute.

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

$$ \frac{\rho}{\theta} \left[ \frac{\partial \psi}{\partial \mathbf{E}} \mathbf{\dot{E}} + \frac{\partial \psi}{\partial \mathbf{\dot{E}}}:\mathbf{\ddot{E}} + \frac{\partial \psi}{\partial \theta} \dot{\theta} + \frac{\partial \psi}{\partial \nabla \theta} \nabla \dot{\theta} +  s \dot{\theta} \right] + \frac{1}{\theta} \mathbf{T}: \dot{\mathbf{E}} - \frac{1}{\theta^2}Q \cdot \rho \theta \ge 0 $$.     (1.3)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

Collect similar terms, and recognize that the above must hold for all values of $$\displaystyle \mathbf{\dot{E}}$$ and $$\displaystyle \mathbf{\ddot{E}}$$. In other words, our constitutive relations are determined.

 * $$\displaystyle \frac{\partial \psi}{\partial \mathbf{\dot{E}}} = 0 \qquad, \qquad \frac{\partial \psi}{\partial \nabla \theta} = 0 \qquad \Rightarrow \qquad \psi = \psi (\mathbf{E}, \theta)$$
 * $$\displaystyle -\frac{\partial \psi}{\partial \mathbf{E}} + \frac{1}{\rho} \mathbf{T} = 0 \qquad \Rightarrow \qquad \mathbf{T}=\rho \frac{\partial \psi}{\partial \mathbf{E}}$$
 * $$\displaystyle \frac{\partial \psi}{\partial \theta} + s = 0 \qquad \Rightarrow \qquad s=- \frac{\partial \psi}{\partial \theta}$$
 * $$\displaystyle -Q \cdot \nabla \theta \ge 0$$

Description
It is often said that for an isotropic material the principle axes of stress coincide with the principle axes of strain. Show that this is true for an elastic solid governed by $$\displaystyle T_{ij}' = 2G \epsilon_{ij}'$$, but not necessarily in a viscoelastic material governed by $$\displaystyle T_{ij}' = 2G \epsilon_{ij}' + 2 \eta \frac{d \epsilon_{ij}'}{dt}$$.

Solution
Principle stresses are defined as

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

$$\displaystyle (T_{ij} - T_{(k)} \delta_{ij})n_j^{(k)} = 0 \qquad \text{for} \quad i,j,k = 1,2,3 $$     (2.1)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

 * where $$\displaystyle T_{(k)}$$ is the principle stress and $$\displaystyle n_j^{(k)}$$ is the principle direction. In matrix form,

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

$$\displaystyle \begin{bmatrix} T_{11}-T_{(k)} & T_{12} & T_{13}\\ T_{21} & T_{22}-T_{(k)} & T_{23}\\ T_{31} & T_{32} & T_{33}-T_{(k)} \end{bmatrix} \cdot \begin{bmatrix} n_1^{(k)}\\ n_2^{(k)}\\ n_3^{(k)} \end{bmatrix} = 0 $$     (2.2)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

The determinant in the characteristic equation now takes the form

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

$$\displaystyle T_{(k)}^3 - \mathrm{I}_T T_{(k)}^2 + \mathrm{II}_T T_{(k)} - \mathrm{III}_T =0 $$     (2.3)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

The stress can be decomposed to

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

$$\displaystyle T_{ij}= \delta_{ij} \frac{T_{kk}}{3} + T_{ij}' $$     (2.4)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

 * where $$\displaystyle T_{ij}'$$ is the deviatoric stress and $$\displaystyle T_{kk}$$ is the mean stress.

Therefore, $$\displaystyle T_{ij}$$ and $$\displaystyle T_{ij}'$$ are related. Since the relationship between deviatoric stress and deviatoric strain is linear (see governing equation for an isotropic material above), the axes for principle strain coincide with principle stress. For this case of a rate-dependent viscoelastic material, the relationship between stress and strain is no longer linear. So the axes of principle stress and principle strain do not necessarily coincide.

Description
If the reference axes are principal axes of stress in an isotropic elastic solid governed by $$\displaystyle \sigma_{ij}' = 2G \epsilon_{ij}'$$ and $$\displaystyle p = -K e$$, combine these equations to obtain three equations relating $$\displaystyle \sigma_1, \ \sigma_2, \ \sigma_3, \ \epsilon_1, \ \epsilon_2, \ \epsilon_3$$.

Solution
From Problem 2 it was shown that the principle stress axes coincide with the principle strain axes for an elastic material that follows Hooke's law. So when $$\displaystyle \sigma_{ij}' = 2G \epsilon_{ij}'$$, we recall equation 2.4 such that

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

$$\displaystyle \sigma_k - p = 2G \left( \epsilon_k - \frac{1}{3}e \right ) $$     (3.1)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

Substituting $$\displaystyle p = -K e$$ into this equation yields

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

$$\displaystyle \sigma_k = 2G \left( \epsilon_k - \frac{1}{3}e \right ) - Ke \qquad \Rightarrow \qquad \sigma_k = 2G \epsilon_k - e \left( \frac{2G}{3} + K\right ) $$     (3.2)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

 * or expanded as

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

$$\displaystyle \sigma_k = 2G \epsilon_k - \left( \frac{2G}{3} + K\right ) ( \epsilon_1 +\epsilon_2 + \epsilon_3) $$     (3.3)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

Therefore the three equations relating $$\displaystyle \sigma_1, \ \sigma_2, \ \sigma_3, \ \epsilon_1, \ \epsilon_2, \ \epsilon_3$$ are

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

$$\displaystyle \begin{matrix} \sigma_1 = \left( \frac{4G}{3} - K\right ) \epsilon_1 - \left( \frac{2G}{3} + K\right ) ( \epsilon_2 + \epsilon_3) \\ \\ \sigma_2 = \left( \frac{4G}{3} - K\right ) \epsilon_2 - \left( \frac{2G}{3} + K\right ) ( \epsilon_1 + \epsilon_3) \\ \\ \sigma_3 = \left( \frac{4G}{3} - K\right ) \epsilon_3 - \left( \frac{2G}{3} + K\right ) ( \epsilon_1 + \epsilon_2) \end{matrix} $$     (3.4)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 *
 * }

Description
Show that the isotropic Hooke's law may be written as

 * $$\displaystyle T_{ij}=2G \left( E_{ij}+\frac{\nu}{1-2\nu} e \delta_{ij}\right )$$

And hence derive the following form for the equilibrium equations of an elastic body in terms of the displacements, neglecting the distinction between spatial and material coordinates.

 * $$\displaystyle G \left(\nabla^2 u_i + \frac{1}{1-2\nu} \frac{\partial^2 u_j}{\partial x_j \partial x_i} \right ) + \rho b_i = 0$$

Solution
Isotropic Hooke's law is defined as

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

$$\displaystyle T_{ij} = \lambda E_{kk} \delta_{ij} + 2 \mu E_{ij} $$.     (4.1)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

 * where Ekk is the volume strain (e).

The two Lamé constants $$ \lambda$$ and $$ \mu$$ are related to the shear and Young's modulii by

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

$$\displaystyle \mu = G = \frac{E}{2(1+\nu)} \qquad \text{and} \qquad \lambda = \frac{\nu E}{(1+\nu)(1-2 \nu)} $$.     (4.2)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

By combining equation 4.1 with the equations of 4.2,

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

$$\displaystyle T_{ij} = \frac{\nu E}{(1+\nu)(1-2 \nu)} e \delta_{ij} + 2G E_{ij} $$     (4.3)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

 * or

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

$$\displaystyle T_{ij} = \quad 2G \left(E_{ij} + \frac{\nu}{1-2\nu} e \delta_{ij} \right) $$     (4.4)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 *
 * }

Now recall the conservation of linear momentum.

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

$$\displaystyle \frac{\partial T_{ji}}{\partial x_j} + \rho b_i = \rho \frac{dv_i}{dt} $$     (4.5)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

Combine equations 4.4 and 4.5.

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

$$\displaystyle \frac{\partial}{\partial x_j} \left[ 2G \left(E_{ij} + \frac{\nu}{1-2\nu} e \delta_{ij} \right) \right] + \rho b_i = \rho \frac{dv_i}{dt} $$     (4.6)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Substitute the gradient expression.

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

$$\displaystyle G \nabla \left( 2E_{ij} + \frac{2 \nu}{1-2 \nu} e \delta_{ij} \right) + \rho b_i = \rho \frac{dv_i}{dt} $$     (4.7)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Now move the right-hand side of the equation to the left and simplify.

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

$$\displaystyle G \left(\nabla^2 u_i + \frac{1}{1-2\nu} \frac{\partial^2 u_j}{\partial x_j \partial x_i} \right ) + \rho b_i = 0 $$     (4.8)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 *

 * }

Problem 5: Objectivity
Determine the objectivity of the following quantities. Refer to the restrictions for vectors and tensors, respectively:

 * $$\displaystyle a_k' = Q_{kl} a_l$$
 * $$\displaystyle S_{kl}' = Q_{km}S_{mn}Q_{nl}$$

Description
Show that for an orthotropic elastic continuum (three orthogonal planes of elastic symmetry) the elastic coefficient matrix is given as

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

$$ [C_{KM}]=\begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0\\ C_{21} & C_{22} & C_{23} & 0 & 0 & 0\\ C_{31} & C_{32} & C_{33} & 0 & 0 & 0\\ 0 & 0 & 0 & C_{44} & 0 & 0\\ 0 & 0 & 0 & 0 & C_{55} & 0\\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{bmatrix} $$.     (6.1)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Solution
The generalized Hooke's law includes nine equations expressing the stress components as linear homoegeneous functions of the nine strain components,
 * $$\displaystyle T_{IJ} = C_{IJKM} E_{KM}$$.

The matrix form is shown below. Recall that the stress tensor and strain tensor are symmetric.
 * $$\displaystyle

\begin{bmatrix} T_{11}\\ T_{22}\\ T_{33}\\ T_{23}\\ T_{31}\\ T_{12} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16}\\ C_{21} & C_{22} & C_{23} & C_{24} & C_{25} & C_{26}\\ C_{31} & C_{32} & C_{33} & C_{34} & C_{35} & C_{36}\\ C_{41} & C_{42} & C_{43} & C_{44} & C_{45} & C_{46}\\ C_{51} & C_{52} & C_{53} & C_{54} & C_{55} & C_{56}\\ C_{61} & C_{62} & C_{63} & C_{64} & C_{65} & C_{66}\\ \end{bmatrix} \begin{bmatrix} E_{11}\\ E_{22}\\ E_{33}\\ 2E_{23}\\ 2E_{31}\\ 2E_{12} \end{bmatrix} $$.


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

title: File:Egm6341.s11.team2.HW7.1.2.fig3.png


text:
User talk:Miss sanders

hey r

title: User talk:Miss sanders


text:
The Art of Computer Programming/Mathematical Induction

This section defines mathematical induction, and proves the correctness of Euclid's algorithm.


text:
Discussion of the Wikisource Translation of the Bible

Wikisource is producing a new translation of the Bible. This project will review that translation.

Canon
The Wikisource Bible wrongly includes 3 Maccabees, 4 Maccabees and Psalm 151 in the Apocrypha, although these books have never been recognised by the Church of England as part of the Apocrypha or included in the King James Version or any Roman Catholic Bibles. This error appears to have arisen because some printings of the Revised Standard Version and the New Revised Standard Version have an "expanded" collection of "The Apocrypha/Deuterocanonical Books" that includes these books. However, it is clear from The Oxford Companion to the Bible, edited by Bruce Metzger, the chairman of the NRSV translators, that these books are not part of the standard Apocrypha. See Apocrypha.

Genesis

 * 1:1 "At first" should be "In the beginning". The Hebrew is "Bereshis"; the prefix "Be" cannot be translated "At".  Further, the whole thrust of the passage is that this was the beginning, the creation; "at first" does not convey that sense.

 * 1:1 "sky" (Hebrew "shamayim") should be "heaven". "Sky" is "shachaq".

Judges

 * 1:3 "destiny" is not very clear. It means "my territory" or "the territory allotted to me".  It would read better if "likewise" were placed earlier in the verse.


text:
User:Egm6341.s11.team4.KTH/HW7

=Problem 7.2: Optimal Solution of Aircraft Bunt Maneuver =

Refer to lecture slides [[media:nm1.s11.mtg40.djvu|40-1]], [[media:egm6341.s10.mtg40.djvu|40-3]] and [[media:egm6341.s10.mtg41.djvu|41-1]] for the problem statement.

Given
The parameters for this problem are given in both the problem statement in the lectures as well as Subchan and Zbikowski's paper Computational optimal control of the terminal bunt manoeuvre—Part 1: Minimum altitude case.

From Subchan and Zbikowski's paper we have the following:

Equations of Motion:

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

(2.1)
 * $$\displaystyle \dot \gamma = \frac{T - D}{mV}sin(\alpha) + \frac{L}{mV}cos(\alpha) - \frac{g cos(\gamma)}{V}$$
 *
 *
 * }
 * }

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

(2.2)
 * $$\displaystyle \dot V = \frac{T - D}{m}cos(\alpha) - \frac{L}{m}sin(\alpha) - g sin(\gamma)$$
 *
 *
 * }
 * }

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

(2.3)
 * $$\displaystyle \dot x = V cos(\gamma)$$
 *
 *
 * }
 * }

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

(2.4)
 * $$\displaystyle \dot y = V sin(\gamma)$$
 *
 *
 * }
 * }

And,

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

\begin{bmatrix} \gamma \\ v \\ x \\ y \end{bmatrix} $$
 * $$\mathbf{z} =
 * $$\mathbf{z} =
 * }
 * }

Where γ, V, x and y are state variables defined such that γ is the flight path angle, V is the speed, x is the horizontal position and y is the altitude.

Also from Subchan and Zbikowski we have:

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

 * $$\displaystyle D(h,V,\alpha) = \frac{1}{2}C_d \rho V^2 S_{ref}, \,\,\, C_d = A_1 \alpha^2 + A_2\alpha + A_3$$
 * }
 * }
 * }

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

 * $$\displaystyle L(h,V,\alpha) = \frac{1}{2}C_l \rho V^2 S_{ref}, \,\,\, C_l = B_1\alpha + B_2$$
 * }
 * }
 * }

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

 * $$\displaystyle \rho = C_1y^2 + C_2y + C_3$$
 * }
 * }
 * }

Where,

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

\begin{align} & m = 1005 \; kg \\ & g = 9.81 \; m/s^2 \\ & S_{ref} = 0.3376 \; m^2 \\ & A_{1} = -1.9431 \\ & A_{2} = -0.1499 \\ & A_{3} = 0.2359 \\ & B_{1} = 21.9 \\ & B_{2} = 0 \\ & C_{1} = 3.312 \times 10^{-9} \; kg/m^5 \\ & C_{2} = -1.142 \times 10^{-4} \; kg/m^4 \\ & C_{3} = 1.224 \; kg/m^3 \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

The initial conditions are given as:

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

\begin{align} & \gamma_0 = 0 \, rad\\ & V_0 = 272 \, m/s \\ & x_0 = 0 \, m \\ & y_0 = 30 \, m \\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

From the problem statement in the lecture notes we get the profiles for the two control inputs thrust, T(t) and angle of attack α(t):

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

T(t) = \begin{cases} 6000 \; N, & \mbox{for }0 \leq t < 27 \; \mbox{and} \; 33 \leq t \leq 40 \\ 1000 \; N, & \mbox{for }27 \leq t < 33 \end{cases} $$     (2.5)
 * $$\displaystyle
 * $$\displaystyle
 *

 * }
 * }

For part 1:

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

\alpha(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.09}{27-21}(t-21) + 0.09 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-27}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ \frac{-0.13+0.2}{40-33}(t-33) - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$     (2.6)
 * $$\displaystyle
 * $$\displaystyle
 *
 * }
 * }

For part 2:

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

\alpha(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.03}{27-21}(t-21) + 0.03 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-27}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$     (2.7)
 * $$\displaystyle
 * $$\displaystyle
 *
 * }
 * }

Objective
a.) Solve for z(t), for t contained on [0,40] such that successive numerical solutions are of the order 10-6, then plot the state variables as a function of time and compare them to the solution using the MatLab function "ode45".  Also compute

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

 * $$\displaystyle J_1 = \int_0^{t_f} y(t)dt \,\,\, s.t. \,\, y(t_f) = 0$$
 * }
 * }
 * }

Use given alpha for part 1.

b.) Solve for z(t), for t contained on [0,40] such that successive numerical solutions are of the order 10-6, then plot the state variables as a function of time and compare them to the solution using the MatLab function "ode45".

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

 * $$\displaystyle J_2 = \int_0^{t_f} y(t)dt \,\,\, s.t. \,\, y(t_f) = 0$$
 * }
 * }
 * }

Use given alpha for part 2.

Solution
Both parts of the problem will be solved using the same technique. The only difference will be the alpha function input into the Matlab code.

From equation (1) on Lecture Slide [[media:nm1.s11.mtg38.djvu|38-3]] we can derive:

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

$$     (2.8)
 * $$\displaystyle \mathbf{z}_{i+1} = \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right]
 * $$\displaystyle \mathbf{z}_{i+1} = \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right]
 *
 * }
 * }

Where:

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

= \frac{1}{2}\left(\mathbf{z}_{i} + \mathbf{z}_{i+1}\right) + \frac{h}{8}\left(\mathbf{f}(\mathbf{z}_{i}) - \mathbf{f}(\mathbf{z}_{i+1})\right) $$     (2.9)
 * $$\displaystyle \mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}) = \mathbf{z}_{i+1/2}
 * $$\displaystyle \mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}) = \mathbf{z}_{i+1/2}
 *
 * }
 * }

F(zi+1) can then be defined by subtracting the right side of eqn. (2.8) from the left side of eqn. (2.8).

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

$$     (2.10)
 * $$\displaystyle \mathbf{F}(\mathbf{z}_{i+1}) = \mathbf{z}_{i+1} - \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right] = \mathbf{0}
 * $$\displaystyle \mathbf{F}(\mathbf{z}_{i+1}) = \mathbf{z}_{i+1} - \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right] = \mathbf{0}
 *
 * }
 * }

From Lecture Slide [[media:nm1.s11.mtg39.djvu|39-1]] we get the following:

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

$$     (2.11)
 * $$\displaystyle \mathbf{z}_{i+1}^{(k+1)} = \mathbf{z}_{i+1}^{(k)} - \left[\frac{\partial \mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{\partial \mathbf{z}}\right]^{-1}\mathbf{F}(\mathbf{z}_{i+1}^{(k)})
 * $$\displaystyle \mathbf{z}_{i+1}^{(k+1)} = \mathbf{z}_{i+1}^{(k)} - \left[\frac{\partial \mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{\partial \mathbf{z}}\right]^{-1}\mathbf{F}(\mathbf{z}_{i+1}^{(k)})
 *
 * }
 * }

In order to solve the problem we need to calculate the Jacobian as follows:

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

$$     (2.12) The values in the Jacobian can be seen in the MatLab code below.  "Main Script"
 * $$\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{z}} = \left[\frac{\partial F_i}{\partial z_j}\right]
 * $$\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{z}} = \left[\frac{\partial F_i}{\partial z_j}\right]
 *
 * }
 * }

 Function "zdot"

 Function "dfdz"

Function "nonlinode1" and "nonlinode2" to be used with the ode45 function.  Function "nonlinode1"   Function "nonlinode2"

a.) Part 1: Initial values for alpha

b.) Part 2: Modified values for alpha

c.) Compute $$\ J_1$$, $$\ J_2$$ and compare.
$$ \ J_1 = 16803.19 $$

$$ \ J_2 = 9024.61 $$

$$ \ J_1 > J_2 $$

So, using modified input makes better results.

title: User:Egm6341.s11.team4.KTH/HW7


title: File:Fe1.s11.team1.HW7.fig3.1.jpg


text:
User:Eml5526.s11.team5.srv/hw7

=Finding the static and dynamic solution of a bi-unit square using 2D LIBF=

P.D.E

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

$$T div (grad u) + f = \rho \ \frac{\partial^2u}{\partial t^2} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                     (Eq. 7.2.1)
 * }

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

$$ \quad T = 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                     (Eq. 7.2.2)
 * }

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

$$ f (\mathbf{x}) = 1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                     (Eq. 7.2.3)
 * }

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

$$ \quad \rho = 3 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                     (Eq. 7.2.4)
 * }

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

$$ \quad Basis \ Function \ : \ 2D \ Lagrange \ Interpolation \ Basis \ Function $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

Essential Boundary Condition

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

$$ \quad g = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                     (Eq. 7.2.5)
 * }

Static(steady state)

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

1. Find the static solution $$ \quad U^h $$ to $$ \quad 10^{-6} $$ accuracy at the center. $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

Dynamic(transient)

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

1. Fundamental Frequency and Fundamental Period $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

2. Plot $$ u^h(0,t) \ vs \ t \quad  for \ t \ \epsilon \ [0,2T_1] $$ for symmetric and Non-Symmetric $$ \quad u^0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

3. Produce a movie of the vibrating membrane for symmetric and Non-Symmetric $$ \quad u^0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

Solution
Substituting the given values in the P.D.E, we get,
 * {| style="width:100%" border="0"

$$\quad 4 \ div (grad u) \ + \ 1 \ = \ 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ 4 \left ( \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2}\right ) + 1 = 0$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} = \frac{-1}{4} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.6)
 * }

2D LIBF

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

$$ N_I (x,y) = L_{i,m}(x).L_{j,n}(y) \quad Where \ I = i+(j-1)m $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.7)
 * }

Where,

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

$$ L_{i,m}(x) = \prod_{k=1\neq i}^{m}\frac{x-x_k}{x_i - x_k} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }
 * {| style="width:100%" border="0"

$$ L_{j,n}(x) = \prod_{k=1\neq j}^{m}\frac{y-y_k}{y_j - y_k} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

Consider the Weak Form

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

$$  \tilde{m}(w,u) + \tilde{k}(w,u) = \tilde{f}(w) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.8)
 * }

where,

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

$$m(w,u) = \int\limits_\Omega {w\rho c\frac} d\Omega $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ k(w,u) = \int\limits_\Omega {\nabla wK\nabla u} d\Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ f(w,u) = - \int_{\Gamma h} {whd{\Gamma _h}}  + \int\limits_\Omega  {wf} d\Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

Static case
For static case, $$ \frac{\partial u}{\partial t} = 0 $$ the above equation can be written as,

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

$$ \quad K d = F $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

where,

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

$$ {K_{ij}} = \int_ {(\nabla N_I^e)I(\nabla N_J^e)d{w^e}} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.9)
 * }

Using the above equations in MATLAB, we get,

K =  33.3257  -13.5386  -13.5386    2.8639 -13.5386  33.3257    2.8639  -13.5386       -13.5386    2.8639   33.3257  -13.5386         2.8639  -13.5386  -13.5386   33.3257

F = [0.5625  0.5625  0.5625  0.5625]^T

d = [0.0617  0.0617  0.0617  0.0617]^T

Dynamic case
For this case, It is given that the value of mass density $$ \quad \rho \ = \ 3$$

Free Vibration: Solve Eigen Value problem
The equation given to solve the Eigen value problem is,

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

$$ \mathbf{K} \mathbf{\phi} = \lambda \mathbf{M}\mathbf{\phi}  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.10)
 * }

Where,

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

$$\quad \mathbf{\phi} - Eigen \ vectors $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \quad \lambda - Eigen \ values $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

Using the above equation and the expression for K, F matrix, we can calculate the fundamental frequency and fundamental time period.

Fundamental Frequency =

Fundamental Time Period =

Transient Analysis with symmetric Uo
Initial conditions

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

$$u(\mathbf{x},t=0) = u^h_s(x) \ \forall \  \mathbf{x} \ \epsilon  \ \Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.11)
 * }

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

$$ \dot{u} (\mathbf{x},t=0) = 0 \ \forall \  \mathbf{x} \ \epsilon  \ \Omega $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.12)
 * }

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

$$T = 4, \qquad f(\mathbf{x},t) = 0 \ in \  \Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.13)
 * }

By using the above initial condition, we can produce the movie of the vibrating membrane.

Movie of the vibrating membrane
The Movie of dynamic membrane has been uploaded in the youtube and it can be seen from the link below

FE Analysis of a Dynamic Membrane - symmetric case

For T = 0,

Transient Analysis with Non-symmetric Uo
In this case, initial u value is given by,

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

$$u(\mathbf{x},t=0) = (x+1)(y+\frac{1}{2}) cos (\frac{\pi}{2}x) cos (\frac{\pi}{2}y) \ \forall \  \mathbf{x} \ \epsilon  \ \Omega$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.14)
 * }

By using the above initial condition, we can produce the movie of the vibrating membrane.

Movie of the vibrating membrane
The Movie of dynamic membrane for the non-symmetric case, has been uploaded in the youtube and it can be seen from the link below

FE Analysis of a Dynamic membrane - Non Symmetric case

title: User:Eml5526.s11.team5.srv/hw7


title: File:Fe1.s11.team1.HW7.fig3.2.jpg


text:
User:Eml5526.s11.team01/Hwk7

=Problem 7.1: Solving on Bi-unit Square=

Given
$$\displaystyle \kappa =I$$, $$\displaystyle ~\rho c=3$$, $$\displaystyle g=2$$ on $$\displaystyle \partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$\displaystyle f =1$$, steady state.

Find
Solve the differential equation

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

$$\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |

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

using LIBF

Solution The weak form:
 * {| style="width:100%" border="0"

$$\mathop{\iint }^{}w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{\kappa }_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)+f-\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }=0$$ Let’s calculate the derivative of the following factor
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.2) $$
 * }
 * }
 * {| style="width:100%" border="0"

$$\frac{\partial }\left( w{{K}_{ij}}\frac{\partial u} \right)=\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+w\frac{\partial }\left( {{K}_{ij}}\frac{\partial u} \right)$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.3) $$
 * }
 * }

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

$$w\frac{\partial }\left( {{K}_{ij}}\frac{\partial u} \right)=\frac{\partial }\left( w{{K}_{ij}}\frac{\partial u} \right)-\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.4) $$
 * }
 * }

Replacing in the integral we obtain
 * {| style="width:100%" border="0"

$$\mathop{\iint }^{}\left( \frac{\partial }\left( w{{K}_{ij}}\frac{\partial u} \right)-\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+wf\left( x,t \right)-w\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.5) $$
 * }
 * }

Appling the Gauss theorem on the first term we obtain
 * {| style="width:100%" border="0"

$$\mathop{\int }^{}w{{K}_{ij}}\frac{\partial u}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }+\mathop{\iint }^{}\left( -\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+wf\left( x,t \right)-w\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }=0$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.6) $$
 * }
 * }

Recognizing that the heat flux is given by

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

$${{q}_{i}}=-{{\text{ }\!\!\Kappa\!\!\text{ }}_{ij}}\frac{\partial u}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.7) $$
 * }
 * }

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

$$\mathop{\int }^{}-w{{q}_{i}}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }+\mathop{\iint }^{}\left( -\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+wf\left( x,t \right)-w\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }=0$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.8) $$
 * }
 * }

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

$$\mathop{\iint }^{}\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}-\mathop{\int }^{}wH(u-{{u}_{\infty }})d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.9) $$
 * }
 * }

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

$$\mathop{\iint }^{}\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}wHud{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}+\mathop{\int }^{}wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.10) $$
 * }
 * }

For the case when the conductivity is the identity matrix we have
 * {| style="width:100%" border="0"

$$\mathop{\iint }^{}\frac{\partial w}\left( {{\delta }_{ij}}\frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}wHud{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}+\mathop{\int }^{}wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.11) $$
 * }
 * }

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

$$\mathop{\iint }^{}\frac{\partial w}\left( \frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}wHud{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}+\mathop{\int }^{}wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.12) $$
 * }
 * }

But in this case we do not have natural boundary condition, therefore we have

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

$$\mathop{\iint }^{}\frac{\partial w}\left( \frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.13) $$
 * }
 * }

Performing the discretization on the dominium we have $$u=\left[ N \right]{{u}_{e}}$$ and $$w=\left[ N \right]{{w}_{e}}$$, where $$\left[ N \right]$$is the vector of interpolation functions (Lagrange polynomial) and $$~{{u}_{e}}$$ and $${{w}_{e}}$$ are the nodal values o u and w respectively. The derivatives with respect to $${{x}_{1}}$$ and $${{x}_{2}}$$ are
 * {| style="width:100%" border="0"

$$~\left[ \begin{matrix} \frac{\partial u}{\partial {{x}_{1}}} \\ \frac{\partial u}{\partial {{x}_{2}}} \\ \end{matrix} \right]=\left[ \begin{matrix} \frac{\partial N}{\partial {{x}_{1}}} \\ \frac{\partial N}{\partial {{x}_{2}}} \\ \end{matrix} \right]{{u}_{e}}=\left[ B \right]{{u}_{e}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.14) $$
 * }
 * }

Replacing on the semi-discrete weak form we obtain

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

$${{w}^{eT}}\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }{{u}^{e}}+{{w}^{eT}}\mathop{\int }^{}\rho c{{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }{{\dot{u}}_{e}}={{w}^{eT}}\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.15) $$
 * }
 * }

From the last equation we can obtain the conductivity matrix, the heat source vector and the capacitance matrix. The conductivity matrix is given by

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

$$\overset{}{\mathop{k}}\,=\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.16) $$
 * }
 * }

The capacitance matrix is given by:

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

$$\overset{}{\mathop{M}}\,=~\mathop{\int }^{}\rho c{{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.17) $$
 * }
 * }

And the heat source is

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

$$\overset{}{\mathop{f}}\,=\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.18) $$
 * }
 * }

To solve the equation:

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

$$Ku+M\dot{u}=f$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.19) $$
 * }
 * }

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

$$\dot{u}={{M}^{-1}}\left( f-Ku \right)$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.20) $$
 * }
 * }

To solve to u from u ̇ we used the Crank-Nicolson scheme where:
 * {| style="width:100%" border="0"

$${{u}_{i+1}}={{u}_{i}}+\text{ }\!\!\Delta\!\!\text{ }t{{\dot{u}}_{i+\alpha }}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.21) $$
 * }
 * }

And:

For this scheme $$\alpha =\frac{1}{2}$$, therefore solving we have

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

$${{u}_{i+1}}={{u}_{i}}+\frac{\text{ }\!\!\Delta\!\!\text{ }t}{2}\left( {{{\dot{u}}}_{i}}+{{{\dot{u}}}_{i+1}} \right)$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.23) $$
 * }
 * }

Part 1: Steady state
$$\displaystyle \frac{\partial u}{\partial t}=0$$ and $$f\left( x \right)=1$$ using LIBF. To solve the problem using LIBF we divide the x axis and y axis in 2 parts initially. In that way the approximated solution using LIBF is given by:

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

$$u=\frac{\left( x-{{x}_{2}} \right)\left( x-{{x}_{3}} \right)}{\left( {{x}_{1}}-{{x}_{2}} \right)\left( {{x}_{1}}-{{x}_{3}} \right)}\frac{\left( y-{{y}_{2}} \right)\left( y-{{y}_{3}} \right)}{\left( {{y}_{1}}-{{y}_{2}} \right)\left( {{y}_{1}}-{{y}_{3}} \right)}{{u}_{1}}+...+\frac{\left( x-{{x}_{1}} \right)\left( x-{{x}_{2}} \right)}{\left( {{x}_{3}}-{{x}_{1}} \right)\left( {{x}_{3}}-{{x}_{2}} \right)}\frac{\left( y-{{y}_{1}} \right)\left( y-{{y}_{2}} \right)}{\left( {{y}_{3}}-{{y}_{1}} \right)\left( {{y}_{3}}-{{y}_{2}} \right)}{{u}_{9}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.24) $$
 * }
 * }

Where $${{u}_{i}}$$ is the nodal value at the node i. the numbering is given in the figure. The last equation can be written as $$u=\left[ N \right]$$.

And [N] can be replaced in the weak form to obtain the solution.

The obtained solution is shown in the next figure (obtained with the routine 1.m) Because lack of memory, the tolerance was only 10^-2

Part 2: Dynamic (Transient)
Solve the differential equation,

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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.25) $$
 * }
 * }

$$\displaystyle \kappa =I$$, $$\displaystyle ~\rho c=3$$, $$\displaystyle  g=2$$ on $$\displaystyle  \partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\displaystyle  \text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where  $$\displaystyle  f =1$$, steady state.

Similar results to the first part were obtained as shown in the figure. The routine 2.m was used to solve the problems. This code can be found in problem 1 appendix.

Part 3: Solving
Solve the differential equation,

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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.26) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, steady state, using LLEBF with NON-uniform mesh.

In this case the analysis was made using 100 elements and no convergence analysis was made. The problem was solved by using the matlab routine 3.m (can be found in the appendix) Inside that routine is included the file Nonunoformmesh.m which contains the connectivity, the coordinates of the nodes and the boundary conditions. The results are shown in the figure 1.4; it can be observed that the results are similar to those obtained with uniform mesh.

Part 4
Solve the differential equation

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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.27) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, transient; $$\text{u}\left( x,0 \right)\text{ }=\text{ xy}$$  using LIBF

Because of memory problems only 16 elements and a tolerance of 10^(-3) was used. Figure 1.5 shows the temperature of the center as function of time obtained with the Crank-Nicholson scheme. It is observed that the temperature tend to the value at the boundaries (2). It is observed also that the temperature decreased at the beginning; this error is probably because of the time step size which probably is too long.

5) Solve the differential equation

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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.28) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, transient; $$\text{u}\left( x,0 \right)\text{ }=\text{ xy}$$   using LLEBF; uniform mesh The initial condition is shown in figure 1.6

We used the Crank-Nicolson scheme before mentioned to solve the transient part of the problem. Initially we used a time step of 0.01 and the solution present some errors as can be observed in the following two figures.

Then we decreased the time step to 0.001 and better results were obtained as shown in the next 3 figures.

The next figure shows the temperature as function of time for the center of the domain point (0,0); it is observed that the temperature converge to the value of 2.07 which is the value at steady state.

The next figure shows the error before the convergence to the steady state. It was found that the error increases at the beginning. This instability could be probably eliminated decreasing more the time step; however that increase the computation time.

Part 6: Solving
Solve the differential equation,

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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 1.29) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, transient; $$\text{u}\left( x,0 \right)\text{ }=\text{ xy}$$  using LLEBF; NON-uniform mesh Solution The next figure shows the mesh used

The initial condition x*y is shown in the next figure

Because the value at steady state with f = 1 is close to the value at the borders (2.07) we increase the source term to f=20; in that way the steady state solution was as shown in the next figure

The next figure shows the evolution of the temperature at the center

The code 6.m, as found in the appendix, was used to solve this problem

Appendix
=Problem 2: Membrane=

Given
Given the membrane square; $$\displaystyle T = 4$$; $$\displaystyle f = 1$$; $$\displaystyle \rho =3$$; $$\displaystyle g = 0$$ on all the boundaries; bi-unit square; using LIBF

Find
Find the steady state and dynamic solution

Part 1

 * The general PDE is:

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

$$Tdiv\left( grad\left( u \right) \right)+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.1) $$
 * }
 * }

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

$$T\nabla .\left( \nabla u \right)+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.2) $$
 * }
 * }

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

$$T{{\nabla }^{2}}(u)+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.3) $$
 * }
 * }

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

$$T\frac{{{\partial }^{2}}u}{\partial x_{i}^{2}}+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.4) $$
 * }
 * }

The weak form is similar to the latter problem. The obtained result is

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

$${{w}^{eT}}\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }{{u}^{e}}+{{w}^{eT}}\mathop{\int }^{}\left( \frac{\rho }{T} \right){{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }{{\dot{u}}_{e}}={{w}^{eT}}\frac{1}{T}\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.5) $$
 * }
 * }

Therefore, the stiffness matrix is given by:

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

$$\overset{}{\mathop{k}}\,=\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.6) $$
 * }
 * }

The mass matrix is:

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

$$\overset{}{\mathop{M}}\,=~\mathop{\int }^{}\left( \frac{\rho }{T} \right){{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.7) $$
 * }
 * }

And the heat source is

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

$$\overset{}{\mathop{f}}\,=\frac{1}{T}\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.8) $$
 * }
 * }

In the static case, the equation is the same as that for heat transfer in steady state; therefore we used the same routine 1.m for solve this problem. In the next figures is shown the solution and the error Convergence of 10^-6 was achieved with 1600 elements after 20 iterations

Part 2

 * Free vibration eigenvalue problem. Find the three first fundamental eigenpairs. For free vibration the factor f is zero, therefore the PDE is reduced to:
 * {| style="width:100%" border="0"

$$\rho \frac{{{\partial }^{2}}u}{\partial x_{i}^{2}}=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.9) $$
 * }
 * }

In 2D we can write:

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

$$\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}u}{\partial {{y}^{2}}}=\frac{\rho }{T}\frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.10) $$
 * }
 * }

This is case where analytical solution can be found; thererfor, we did not use any numerical method; instead, using separation of variables we can write the solution as:
 * {| style="width:100%" border="0"

$$u\left( x,y,t \right)=F\left( x \right)G\left( y \right)H(t)$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.11) $$
 * }
 * }

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

$$F(x)G\left( y \right)H\left( t \right)+F\left( x \right){{G}^{\left( y \right)}}H\left( t \right)=\frac{\rho }{T}F\left( x \right)G\left( y \right)H''(t)$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.12) $$
 * }
 * }

Where ‘’ means second derivative with respect to x, y or t depending if the derivative is of F, of G or of H respectively. If the latter equation is divided by u we obtain
 * {| style="width:100%" border="0"

$${{F}^{}}(x)G\left( y \right)H\left( t \right)+F\left( x \right){{G}^{}}(y)H\left( t \right)=\frac{\rho }{T}F\left( x \right)G\left( y \right)H''(t)$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.13) $$
 * }
 * }

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

$$\frac{F}{F}+\frac{G}{G}={{c}^{2}}\frac{H''}{H}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.14) $$
 * }
 * }

As every fraction in the latter equation depends only on x, y or t respectively they should be equal to a constant. Let’s call that constant w^2
 * {| style="width:100%" border="0"

$$\frac{F}{F}+\frac{G}{G}={{c}^{2}}\frac{H''}{H}=-{{w}^{2}}$$ Now we have the equations:
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.15) $$
 * }
 * }
 * {| style="width:100%" border="0"

$${{c}^{2}}\frac{H''}{H}=-{{w}^{2}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.16) $$
 * }
 * {| style="width:100%" border="0"
 * {| style="width:100%" border="0"

$$\frac{F}{F}=-\frac{G}{G}-{{w}^{2}}$$ In the last equation we again observe that the left hand side depends only on x and the left hand side depends on y therefore should be constant. Let’s call $$k_{1}^{2}$$and $$k_{2}^{2}$$ then
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.17) $$
 * }
 * }
 * {| style="width:100%" border="0"

$$\frac{F''}{F}=-k_{1}^{2}$$ And
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.18) $$
 * }
 * }
 * {| style="width:100%" border="0"

$$\frac{G''}{G}=~-k_{2}^{2}$$ Replacing equations
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.19) $$
 * }
 * }
 * {| style="width:100%" border="0"

$${{w}^{2}}=k_{1}^{2}+k_{2}^{2}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.20) $$
 * }
 * }

From above we have
 * {| style="width:100%" border="0"

$${{F}^{''}}+k_{1}^{2}F=0$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.21) $$
 * }
 * }

This equation has the solution $$F\left( x \right)=Acos\left( {{k}_{1}}x \right)+Bsin\left( {{k}_{1}}x \right)$$ To apply easily the boundary conditions we place the origin of coordinates at the lower left corner (Figure)

From the figure we observe that F(0) = 0, and F(L) = 0 With the first condition we have:
 * {| style="width:100%" border="0"

$$F\left( 0 \right)=Acos\left( 0 \right)+Bsin\left( 0 \right)=0$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.22) $$
 * }
 * }

Therefore A = 0. Also:
 * {| style="width:100%" border="0"

$$F\left( L \right)=Bsin\left( {{k}_{1}}L \right)=0$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.23) $$
 * }
 * }

In the last equation B cannot be zero because if B is zero, as A is zero F(x) would be zero and therefore u would be zero; therefore the $$\sin \left( {{k}_{1}}L \right)$$ has to be zero, which means that $${{k}_{1}}L=n\pi $$ where $$n=1,~2,~3,~\ldots $$ As in our case the length of the domain is 1 (L = 1) we have $${{k}_{1}}=n\pi $$ In similar way we can obtain that $${{k}_{2}}=m\pi $$ From above we have
 * {| style="width:100%" border="0"

$${{H}^{''}}\left( t \right)+\fracH\left( t \right)=0$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.24) $$
 * }
 * }

Therefore the solution is:
 * {| style="width:100%" border="0"

$$H\left( t \right)=Dcos\left( \frac{w}{c}t \right)+Esin\left( \frac{w}{c}t \right)$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.25) $$
 * }
 * }

From the latter equation we observe that the natural frequency is
 * {| style="width:100%" border="0"

$${{w}_{n}}=\frac{w}{c}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.26) $$
 * }
 * }

Replacing we have

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

$${{w}_{mn}}=\frac{1}{c}\sqrt{k_{1}^{2}+k_{2}^{2}}=\sqrt{\frac{T}{\rho }\left( {{\left( \frac{n\pi }{L} \right)}^{2}}+{{\left( \frac{m\pi }{L} \right)}^{2}} \right)}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.27) $$
 * }
 * }

The first 3 natural frequencies are when n = m=1; n = 1 and m = 2; n = 2 and m = 2

With T = 4, $$\rho =3$$, and L = 1 we obtain: w1 = 5.13; w2 = 8.11;  w3 = 10.26

Part 3

 * Transient analysis with f = 0; u(x,0) = u-static

In the dynamic case, the Newmark scheme was used to solve the equation:

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

$$\left[ M \right]\ddot{u}+\left[ K \right]u=F$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.28) $$
 * }
 * }

According to the Newmark scheme the function u and its first derivative are approximated by
 * {| style="width:100%" border="0"

$${{u}_{i+1}}={{u}_{i}}+\Delta t{{\dot{u}}_{i}}+\frac{1}{2}\Delta {{t}^{2}}{{\ddot{u}}_{i+\gamma }}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.29) $$
 * }
 * }
 * {| style="width:100%" border="0"

$${{\dot{u}}_{i+1}}={{\dot{u}}_{i}}+\Delta t{{\ddot{u}}_{i+\alpha }}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.30) $$
 * }
 * }

Where

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

$${{\ddot{u}}_{i+\gamma }}=\left( 1-\alpha \right){{\ddot{u}}_{i}}+\alpha {{\ddot{u}}_{i+1}}$$
 * style="width:95%" |
 * style="width:95%" |
 *
$$\displaystyle (Eq. 2.31) $$
 * }
 * }

The values of $$\alpha =\gamma =\frac{1}{2}$$where used. The problem was solved by using the routine 3.m in the appendix of problem 1. The next figure shows the displacement of the center of the membrane as function of time.

Part 4

 * Transient analysis with f = 0;

The next figure is the displacement vs time of the center.

'''Above are the Images from the Movie of the VIBRATION MEMBRANE. To see the full 1:51 minutes Movie go to the Following link.'''

=Problem 7.3 Static Solution for Unit Circle=

Given
Let $$ \Omega = $$ circle with unit radius.

$$ T = 4 $$

$$ f (x) = 1 in \Omega $$

$$ g = 0 on \Gamma_g = \partial \Omega to 10^{-6} $$ accuracy at center.

Find
For both quad and triangle elements.

sym $$ \Rightarrow $$ use $$~ \frac{1}{4} ~$$ of meshes.

Plot deformed shape in 3-D.

Quad)
Using SolidWorks software. A unit circle was created, with a mesh of 9 elements, all quadrilateral. It was from the data points that the mesh nodes were created by hand using the measurement tool. The picture below shows the mesh construction on SolidWorks.

After obtaining the coordinates from the unit circle and numbering the nodes from 1 to 16 since there are 9 elements, the boundary conditions were defined.

The boundary on the bottom is defined as

and for the left hand side of the circle, the boundary condition is as follows:

This means that these boundaries are "free" and the nodes are open to deflection. However the arc of the circle is not "free" these boundary conditions are locked and the nodes stay in one place.

After obtaining this information, a Matlab code was created to re-create the nodes and boundary conditions for the structure. The code is listed below.

The following code then calls the mesh file and makes a 3D plot of the unit circle.

Here is the result of the 3D plot. The arc of the circle does not deflect, but the free ends of the circle are clearly displaced.

Triangular)
The procedure is done the same way except with a different Matlab code, shown below for the mesh and plot.

=Problem 7.4: Newtons Law of Colling=

Given
The partial Differential Equation is

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

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

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

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

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

Further the Capacitance Matrix is Given by,

And the Heat Source is

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

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

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

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

Comment:

We can see from the above figures that there also exists a negative value, which is the same with our results from HW6.10.

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

title: User:Eml5526.s11.team01/Hwk7


text:
User:Eml5526.s11.team3.hylon/Homework 7

=Problem 7.2 Membrane Analysis=

Problem Statement
Refer to Lecture notes [[media:fe1.s11.mtg40.djvu|Lecture 40-5]] [[media:fe1.s11.mtg41.djvu|Lecture 41-1,41-2]] for more detailed description.

title: User:Eml5526.s11.team3.hylon/Homework 7


title: File:Fe1.s11.team3.HW7.2.fig1.jpg


title: File:Fe1.s11.team3.HW7.2.fig2.jpg


title: File:Fe1.s11.team3.HW7.2.fig3.jpg


title: File:Fe1.s11.team3.HW7.2.fig5.jpg


title: File:Fe1.s11.team3.HW7.2.fig4.jpg


title: File:Fe1.s11.team3.HW7.2.fig6.jpg


text:
User:Egm6341.s11.team4.kurth/hw7

Exact Solution
Starting with the logistic equation and separating the variables, it becomes

 * $$\displaystyle \frac{1}{x(1-\frac{x}{x_{max}})}\;dx=r\;dt$$

Now, multiply the left-hand side of the equation by $$\frac{x_{max}}{x_{max}}$$  to simplify the denominator

 * $$\displaystyle \frac{x_{max}}{x(x_{max}-x)}\;dx=r\;dt$$

The left-hand side can now be deconstructed into the partial fractions

 * $$\displaystyle \left[\frac{1}{x}+\frac{1}{x_{max}-x}\right]dx=r\;dt$$

Integrating both sides yields

 * $$\displaystyle ln|x|-ln|x_{max}-x|-ln|x_0|+ln|x_{max}-x_0|=r(t-t_0)$$

Taking the inverse natural log of both sides

 * $$\displaystyle \frac{x(x_{max}-x_0)}{x_0(x_{max}-x)}=e^{r(t-t_0)}$$

Now algebraically solving for x

 * $$\displaystyle x(x_{max}-x_0)=x_0(x_{max}-x)e^{r(t-t_0)}$$

 * $$\displaystyle x(x_{max}-x_0)=x_0x_{max}e^{r(t-t_0)}-x_0xe^{r(t-t_0)}$$

 * $$\displaystyle x(x_{max}-x_0+x_0e^{r(t-t_0)})=x_0x_{max}e^{r(t-t_0)}$$

 * $$\displaystyle x = \frac{x_0x_{max}e^{r(t-t_0)}}{x_{max}-x_0+x_0e^{(t-t_0)}}$$

title: User:Egm6341.s11.team4.kurth/hw7


text:
User:Eml5526.s11.team3.ushnish12/Homework 7

Problem statement
Please refer to lecture note 37-2 and 38 for more details.

Solve the heat conduction problem in 2D with boundary convection.

1> Construct the weak form for heat conduction in 2D with boundary convection.

2> Develop semi-discrete equations (ODEs) similar to (3) 23-3 Give detailed expression of all quantities.

3>Solve G2DM2.0/D1 using 2D LIBF till $$\displaystyle 10^{-6} $$ accuracy at the center. Plot solution in 3D with contour.

Developing the Weak Form of Newton's law of Cooling
i> From Lecture 33-2, the PDE of the problem is given by:-

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

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

\frac{\partial }\left( {{K}_{ij}}\frac{\partial u} \right)+f\left( x,t \right)=\rho c\frac{\partial u}{\partial t}

$$ (Eq )<1>
 *
 * }

From the PDE the weak form can be written as:-

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

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

\int{w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)+f\left( x,t \right)-\rho c\frac{\partial u}{\partial t} \right)}d\text{ }\!\!\Omega\!\!\text{ =0}

$$ (Eq )<2>
 *
 * }

Now from the theory of Integration by parts, we can write

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

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

\int{w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right) \right)}d\text{ }\!\!\Omega\!\!\text{ =}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial }{\partial {{x}_{i}}}\left( w{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}}d}\text{ }\!\!\Omega\!\!\text{ }

$$
 * }

Hence equation 2 can be now written as:-

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

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

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial }{\partial {{x}_{i}}}\left( w{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0

$$ (Eq )<3>
 *
 * }

Now we can write $$\displaystyle d\Omega $$ as $$\displaystyle d\Omega ={{\Gamma }_{g}}\cup {{\Gamma }_{h}}\cup {{\Gamma }_{H}}$$ , where

$$\displaystyle {{\Gamma }_{g}}$$ = essential boundary condition. $$\displaystyle {{\Gamma }_{h}}$$ = natural boundary condition. $$\displaystyle {{\Gamma }_{H}}$$ = Mixed boundary condition.

So by applying the Gauss theorem on the first term we obtain,

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

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

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}w{{K}_{ij}}\frac{\partial u}{\partial{{x}_{j}}}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0

$$ (Eq )<4>
 *
 * }

Now the heat flux is given by:-

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

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

{{q}_{i}}=-{{\text{ }\!\!\Kappa\!\!\text{ }}_{ij}}\frac{\partial u}{\partial {{x}_{j}}}

$$ (Eq )<5>
 *
 * }

So by rearranging Equation 5, we can write

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

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

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{-w{{q}_{i}}{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0

$$ (Eq )<6>
 *
 * }

We can write the first term as-

$$\displaystyle \int\limits_{\Gamma }{-w{{q}_{i}}{{n}_{i}}d\Gamma =}\int\limits_{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{g}}}-\int\limits_{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{h}}}-\int\limits_{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{H}}}$$

We select $$\displaystyle w $$ such that $$\displaystyle w=0 $$ on $$\displaystyle {{\Gamma }_{g}}$$

On $$\displaystyle {{\Gamma }_{h}},\quad  q(x,t)n(x)=h(x,t) ,\forall x\in {{\Gamma }_{h}}$$

$$\displaystyle {{\Gamma }_{H}},\quad q(x,t)n(x)=H[u(x,t)-{u}_{\infty}] ,\forall x\in {{\Gamma }_{H}}$$

So we can write Equation 6 as :-

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

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

-\int\limits_{whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}-\int\limits_{wH(u-{{u}_{\infty }})d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)d\text{ }\!\!\Omega\!\!\text{ }}+\int\limits_{\Omega }{wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-}\int\limits_{\Omega }{w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=0}

$$ (Eq )<7>
 *
 * }

This can be re-arranged in the form:-

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

$$\displaystyle \int\limits_{\Omega }{w\rho c\frac{\partial u}{\partial t}d\Omega +\int\limits_{\Omega }{\nabla w.\kappa .\nabla ud\Omega +}}\int\limits_{wHu}d{{\Gamma }_{H}}=-\int\limits_{wh}d{{\Gamma }_{h}}+\int\limits_{wH{{u}_{\infty }}}d{{\Gamma }_{H}}+\int\limits_{\Omega }{wf(x,t)}d\Omega $$ (Eq )<8>
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 *
 * }

This Equation is of the form:-

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

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

\tilde{m}(w,u)+\tilde{k}(w,u)=\tilde{f}(w) \quad \forall w $$ such that $$\displaystyle w=0 \quad on \quad {{\Gamma }_{g}}

$$
 * }

where,

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

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

\tilde{m}(w,u)=\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }}$$ $$\displaystyle\tilde{k}(w,u)=\int\limits_{\Omega }{\nabla w.\kappa .\nabla ud\Omega }+\int\limits_{wHu}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$ $$\displaystyle\tilde{f}(w)=-\int\limits_{whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}+\int\limits_{wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}+\int\limits_{\Omega }{wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }}

$$ (Eq )<9>
 *
 * }

Non-discrete

To form the semi-discrete equations, we take $$\displaystyle w \quad and \quad u$$ as:-

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

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

w^h=\sum\limits_{i} \quad and \quad u^h=\sum\limits_{j}

$$ (Eq )<10>
 *
 * }

Accordingly the mass, capacitance and heat source equations can be written as:- Capacitance matrix

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

$$\displaystyle \overset{\tilde{\ }}{\mathop{m}}\,(w^h,u^h)=\int\limits_{\Omega }{w^h\rho c\frac{\partial u^h}{\partial t}d\Omega } $$ $$\displaystyle=\int\limits_{\Omega }{\sum{{{c}_{i}}{{b}_{i}}\rho c\sum{{{d}_{j}}{{b}_{j}}d\Omega }}}$$ $$\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}\left( \int\limits_{\Omega }{{{b}_{i}}\rho c{{b}_{j}}}d\Omega \right)}}d_{j}^{'} $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \overset{\tilde{\ }}{\mathop{M}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{{{b}_{i}}\rho c{{b}_{j}}}d\Omega $$ (Eq )<11> Conductivity matrix
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 *
 * }

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

$$\displaystyle k(w^h,u^h)=\int\limits_{\Omega }{\nabla w^h.\kappa .\nabla u^hd\Omega +\int\limits_{w^hHu^hd{{\Gamma }_{H}}}} $$ $$\displaystyle=\int\limits_{\Omega }{\frac{\partial \sum}{\partial {{x}_{i}}}}.\kappa .\frac{\partial \sum}{\partial {{x}_{j}}}d\Omega +\int\limits_{\sumH\sumd{{\Gamma }_{\Eta }}} $$ $$\displaystyle=\sum\limits_{i}{\sum\limits_{j}\int\limits_{\Omega }{\left( \frac{\partial {{b}_{i}}}{\partial {{x}_{i}}}.\kappa \frac{\partial {{b}_{j}}}{\partial {{x}_{j}}}d\Omega \right)}}{{d}_{j}}+\sum\limits_{i}{\sum\limits_{j}\int\limits_{\left( {{b}_{i}}H{{b}_{j}}d{{\Gamma }_{\Eta }} \right)}}{{d}_{j}}
 * style="width:95%" |
 * style="width:95%" |

$$
 * }

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

$$\displaystyle \overset{\tilde{\ }}{\mathop{K}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{b_{i}^{'}\kappa b_{j}^{'}d\Omega +\int\limits_{{{b}_{i}}H{{b}_{j}}d}}{{\Gamma }_{\Eta }} $$ (Eq )<12> Force matrix
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 *
 * }

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

$$\displaystyle \overset{\tilde{\ }}{\mathop{f}}\,(w^h)=-\int\limits_{w^hhd{{\Gamma }_{h}}}+\int\limits_{w^hH{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{w^hf(x,t)d\Omega }$$ $$\displaystyle =-\int\limits_{\sum\limits_{i}hd{{\Gamma }_{h}}}+\int\limits_{\sum\limits_{i}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{\sum\limits_{i}f(x,t)d\Omega }$$ $$\displaystyle -\sum\limits_{i}\int\limits_{{{b}_{i}}hd{{\Gamma }_{h}}}+\sum\limits_{i}\int\limits_{{{b}_{i}}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\sum\limits_{i}\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega } $$
 * style="width:95%" |
 * style="width:95%" |
 * }

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

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

\overset{\tilde{\ }}{\mathop{F}}\,(b_i)=-\int\limits_{{{b}_{i}}hd{{\Gamma }_{h}}}+\int\limits_{{{b}_{i}}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }

$$ (Eq )<13>
 *
 * }

3D Contour plots
iii>

consider the case of n=m=3 the number of nodes is 9 The trial solution is
 * {| style="width:100%" border="0"

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

u_{p}^{h}(x)={{d}_{1}}*N_{1}+{{d}_{2}}N_{2}+{{d}_{3}}N_{3}+\cdots +{{d}_{9}}N_{9},_ – ^ – j=0,1,2,\cdots ,9

$$
 * }

Since the trial solution must satisfy the essential boundary condition,, we should have  $${{d}_{1}}=2,{{d}_{2}}=2,{{d}_{3}}=2,{{d}_{6}}=2,{{d}_{9}}=2$$. , because it is specified that through out the boundary the temperature is 2 i.e for all the nodes on the boundary temperature is 2

In order to solve for the rest coefficient $${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the Conductance matrix, Heat Source matrix and Capacitance matrix as given by Equations 12,13 and 14.

Using Matlab to construct above matrices until satisfying the convergence criterion, we finally have
 * {| style="width:100%" border="0"

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

=\left( \begin{matrix} 151/200 &  -1/5 &  -1/30 & -133/1000 & -16/45 &    1/9 & -33/500 &    1/9 &  -1/45\\ -1/5 &  88/45 &   -1/5 & -16/45 & -16/15 & -16/45 & 1/9 & 0 & 1/9\\ -1/30 & -1/5 & 28/45 & 1/9 & -16/45 & -1/5 &  -1/45 & 1/9 & -1/30\\ -133/100 & -16/45 & 1/9 & 112/45 & -16/15 & 0 & -133/1000 & -16/45 & 1/9\\ -16/45 & -16/15 & -16/45 & -16/15  & 256/45 & -16/15 & -16/45 & -16/15 & -16/45\\ 1/9 & -16/45 & -1/5 & 0 & -16/15 & 88/45 & 1/9 & -16/45 & -1/5\\ -33/500 & 1/9 & -1/45 & -133/1000 & -16/45 & 1/9 & 1511/2000 & -1/5 & -1/30\\ 1/9 & 0 & 1/9 & -16/45 & -16/15 & -16/45 & -1/5 & 88/45 & -1/5\\ -1/45 & 1/9 & -1/30 & 1/9 & -16/45 & -1/5 & -1/30 & -1/5 & 28/45\\ \end{matrix} \right)

$$
 * }

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

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

{F}=\left( \begin{matrix} 166/1000\\ 0\\ 0\\ 33/500\\ 0\\ 0\\ -59/60\\ -4\\ -1\\ \end{matrix} \right)

$$
 * }

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

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

{d}=\left( \begin{matrix} 2.000\\ 2.000\\ 2.000\\ 0.335\\ 0.847\\ 2.000\\ -1.460\\ -1.331\\ 2.000\\ \end{matrix} \right)

$$
 * }

The trial solution using 2D LIBF obtained is as follows:-

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

$$\displaystyle \begin{align} &U_{F}^{h}=(x)0.50*x*y*(x - 1)*(y - 1) - 1*x*(1.0*x + 1.0)*(1.0*y + 1.0)*(y - 1)\\& \quad + 0.665*y*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1) - 0.1678*x*(1.0*y + 1.0)*(x - 1)*(y - 1)\\& \quad - 1.000*y*(1.0*x + 1.0)*(x - 1)*(y - 1) + 0.847*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1)*(y - 1) \\& \quad + 0.500*x*y*(1.0*x + 1.0)*(1.0*y + 1.0) + 0.500*x*y*(1.0*x + 1.0)*(y - 1) - 0.365*x*y*(1.0*y + 1.0)*(x - 1)\\ \end{align} $$ (Eq )<14>
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 *
 * }

The Matlab code to generate the plots and LIBF is:

The 3D contour plot generated from MATLAB is dislpayed below.


text: File:Eml5526.s11.team7.hw7.2a n3.jpg

title: File:Eml5526.s11.team7.hw7.2a n3.jpg


text:
User:Egm6341.s11.team4.hylon/hw7

= Problem 7.3 Solve logistic equation using inconsistent Trapezoidal-Simpson Algorithm =

Given
Refer to Lecture note S10 41-2 for more detailed problem statement

The logistic equation for population dynamics is given in S10 38-3,

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

$$\displaystyle \dot{x}(t) = f(x) = rx \left(1-\frac{x}{x_{max}}\right) $$ $$
 * $$\displaystyle
 * }.
 * }.

Where, $$\displaystyle x_{max}=10;\ r = 1.2,\ $$ and $$\displaystyle t \in [0,10] $$

And the analytical solution is given by,

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

$$\displaystyle x(t) = \frac{x_0x_{max}e^{rt}}{x_{max} + x_0(e^{rt}-1)} $$ $$
 * $$\displaystyle
 * }.
 * }.

Find
Integrate the population function using the inconsistent Trapezoidal-Simpson Algorithm given by,

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

$$\displaystyle z_{i+1/2} = \frac{1}{2} [z_{i} + z_{i+1}] $$
 * }.
 * }.

Two initial conditions, $$\displaystyle x_0=2\ and\ x_0 = 7,\ $$.

Solution
The Hermite-Simpson algorithm behaves better than inconsistent Trapezoidal-Simpson algorithm, that the error is much smaller comparing to inconsistent Trapezoidal-Simpson algorithm. Below are the plots for comparing Hermite-Simpson algorithm with Trapezoidal-Simpson algorithm.

Reference
S10.Team4

title: User:Egm6341.s11.team4.hylon/hw7


text: File:Eml5526.s11.team7.hw7.2a n3a.jpg

title: File:Eml5526.s11.team7.hw7.2a n3a.jpg


text: File:Eml5526.s11.team7.hw7.2a n4.jpg

title: File:Eml5526.s11.team7.hw7.2a n4.jpg


text: File:Eml5526.s11.team7.hw7.2a n4a.jpg

title: File:Eml5526.s11.team7.hw7.2a n4a.jpg


text: File:Eml5526.s11.team7.hw7.2a n5.jpg

title: File:Eml5526.s11.team7.hw7.2a n5.jpg


text: File:Eml5526.s11.team7.hw7.2a n5a.jpg

title: File:Eml5526.s11.team7.hw7.2a n5a.jpg


title: File:Nonuniform6.jpg


title: File:Nonuniform8.jpg


title: File:Nonuniform10.jpg


title: File:Fe1.s11.team1.hw7.fig7.4.1.jpg


title: File:Fe1.s11.team1.hw7.fig7.4.2.jpg


title: File:Nonuniform4.jpg


title: File:Nm1.s11.team5.HW7.fig3a.png


title: File:Nm1.s11.team5.HW7.fig3b.png


title: File:Fe1.s11.team1.hw7.1.Figure1.1.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.2.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.3.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.4.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.5.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.6.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.7.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.8.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.9.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.10.jpg


title: File:Uniform4.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.11.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.12.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.13.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.14.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.15.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.16.jpg


title: File:Fe1.s11.team1.hw7.1.Figure1.17.jpg


title: File:Uniform6.jpg


title: File:Uniform8.jpg


text:
Phenomenology

Phenomenology
Phenomenology studies human experience to discover underlying aspects of this experience. The aim is not find cause and effect, but to understand how things are experienced, lived, interpreted, what people see as significant. Phenomenology focuses on life situations and sees meaning in routine human behaviour.

Phenomenology collects views from participants and describes what participants have in common as they describe a phenomenon (Cresswell et al, 2007). The emphasis is on people's subjective experiences and interpretations, the ‘lived experience’ of the phenomenon.

As well as a reflective research approach, phenomenology is also an influential and complex philosophic tradition (Van Manen, 2007). Different versions of phenomenology developed over time, but they all centre on investigating the essence of human experience. Denscombe (2010) distinguishes a European and a North American approach to phenomenology

The European approach includes (Van Manen, 2007)
 * Transcendental phenomenology - knowledge as coming from insight, interpretative rather than objective (Husserl)
 * Existential phenomenology - observer cannot separate her/himself from observed world (Sartre, De Beauvoir)
 * Hermeneutical (study of interpretation) phenomenology - knowledge reached through language and understanding (Heidegger)
 * Linguistic phenomenology - role of language and discourse in revealing relations between culture, historicality, identity  (Derrida, Foucault)

Moustakas (1990 cited in Richards and Morse, 2007) also distinguish Heuristic phenomenology – experience-based techniques when researchers seek to understand themselves and their lived worlds. This is largely autobiographical research.

The North American approach emanates from the social phenomenology of Alfred Schutz (Denscombe, 2010) and focuses on mental processes by which people make sense of experiences. 'This type of phenomenology is less concerned with revealing the essence of experience, and more concerned with describing the ways in which humans give meaning to their experiences' (Denscombe, 2010, p 101). It has been criticised for accepting to limit itself to describe what is being experienced, without seeking the essence of what is described, and for seeing value in describing experiences just for the sake of describing experience (Denscombe, 2010)

Assumptions
Phenomenology assumes that perceptions provide evidence of the world and that phenomena can be understood, but this requires putting aside already established ideas, beliefs. By going putting them aside, phenomenology aims bring to light the limitations of our culture. While culture help us understand the world, this understanding is limiting and it is important to go beyond.Phenomenology seeks to derive the essence of experience, and to expose taken for granted assumptions (Sokolowski, 2000 cited in Starks and Brown, 2007).

The researcher
Generally phenomenology emphasises subjectivity rather than objectivity. Phenomenology encourages questioning what is taken for granted. The researcher attempts to look at phenomena with fresh eyes, setting aside what we take for granted about these (Crotty, 1998). Existing interpretations should not interfere nor filter through (or as little as possible). Preconceptions have to be consciously reduced (if not eliminated).

Research problems addressed
Understanding the lived experience of a particular phenomenon (e.g in health – exploration of experience of illness – child of someone with illness)

Methods and data collection
Data collection aims to find out how participants experience a phenomenon. A variety of data can be used. Often more than one source is used. Phenomenology typically uses unstructured, conversational interviews, where participants are asked to describe their experience, to pay attention to a phenomenon, talk about it, describe it openly. Methods can be empirical as well as reflective. Reflective methods involves ‘reduction’ or the suspension of everyday attitude and the ‘vocatio’ (letting language, a text speak to us)

Results
Reflection on personal experience, experience of others, and observation that gives insight into meanings.

Resources
Center for Advanced Research in Phenomenology http://www.phenomenologycenter.org/ Colaizzi, Pr. (1973)Marton, F., & Booth, S. (1997). Learning and Awareness. New Jersey: Lawerence Erlbaum Associates. Denscombe,Martyn (2010)(4th ed). The good research guide for small scale social research projects. Maidenhead: Open University Pres McGraw Hill Langdridge, D. (2007). Phenomenological psychology: theory, reasearch and method. Prentice Hall: Pearson. Moustakas, C. (1990). Heuristic research: design, methodology and applications. Newbury Park: Sage. Moustakas, C. (1994). Phenomenological research methods. Thousand Oaks: Sage Phenomenology online http://www.phenomenologyonline.com/inquiry/9.html Starks, Helen and Grown, Susan (2007), ‘Choose your method: a comparison of phenomenology. Discourse analysis and grounded theory’, Qualitative Health Research  17/10, pp 1372-1380 Van Manen, M. (1990). Researching lived experience: human science for an action sensitive pedagogy. State University of New York Press: New York. Van Manen (2007) http://www.phenomenologyonline.com/

Phenomenology and Practice journal is hosted by University of Alberta.

title: Phenomenology


text: File:Egm6341 s11 team5 hw5 graph41.png

title: File:Egm6341 s11 team5 hw5 graph41.png


title: File:Unifrom10.jpg


title: File:Fe1.s11.team1.hw7.2.Figure1.1.jpg


title: File:Fe1.s11.team1.hw7.2.Figure1.2.jpg


title: File:Fe1.s11.team1.hw7.2.Figure1.3.jpg


title: File:Fe1.s11.team1.hw7.2.Figure1.4.jpg


title: File:Fe1.s11.team1.hw7.2.Figure1.5.jpg


title: File:Fe1.s11.team1.hw7.2.Figure1.6.jpg


text:
Discourse analysis

Discourse analysis involves an ‘analysis of the ways in which discourses – which can be read in texts and talk – constitute the social world (Mason, 2006). Developed from linguistics, literary criticism, semiotics, discourse analysis looks at meaning behind ‘text’ or implied meanings’. Discoure analysis situates a text within a context and unpacks what people are implicitly trying to do in a text. It is largely concerned with language, but text can also refer to images, film...The methodology assumes that words and images do not depict reality, but create reality, that words are chosen to have an effect to readers

Unlike grounded theory, discourse analysis works with prior assumptions, since existing knowledge about society informs the analysis.

Discourse analysis can mean different things, since many strands developed over the years. It is therefore important to define how the term is used. Some see discourse analysis as method rather than a research framework or strategy. Indeed there is some overlap between linguistic phenomenology and discourse analysis.

Critical discourse analysis explores how texts serves the interest of powerful groups and how discourse achieves power. Discourse analysis can also examine the blending together of different texts. This assumes that knowledge and meaning are produced through interaction with multiple discourses (Starks and Brown 2007, p 1373).

Assumption
Language is itself meaningless – a system of signs but agreed meaning generates meaning. Words are not determined by what they represent, they are chosen to have an effect (same for images and photo). Language thus reveals background assumptions and has to be examined within the context in which it is produced.

Data collection and methods
Written and spoken texts, images. Method involves examining how language is used to accomplish certain objectives and positions in relation to others. Involves deconstruction of data to show how texts sustain particular ideas about social life, to find out what a text is trying to do and how this is achieved (Denscombe, 2010). Discourse analysis is generally qualitative, but can be empirical and quantitative, e.g. linguistic study of texts scrutinised as separate from their author.

Resources

 * Denscombe, Martyn (2010). The Good Research Guide. (4th ed).Maidenhead: McGraw-Hill
 * Lazaraton, Anne (2002). ‘Quantitative and Qualitative Approaches to Discourse’, Annual Review of Applied Linguistics Vol 22, pp 32 – 51
 * Lazaraton, A. (2009). 'Discourse analysis'. In J. Heigham & R.A. Croker (Eds.). Qualitative research in applied linguistics: A practical introduction (pp. 242-259). Hampshire, UK: Palgrave Macmillan.
 * Potter, J. and Wetherell, M. (1994). ‘Analysing discourse’, in A Bryman and R.G. Burgess (eds). Analysing qualitative data. London: Routledge
 * Woofit, R. (2005).Conversation analysis and discourse analysis: a comparative and critical introduction .London: Routledge

title: Discourse analysis


text:
User:Eml5526.s11.team2.sandhu/hw7

=Problem 7.4=

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

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

$$  \displaystyle Q_{1}=-\int_{\partial w}\mathbf{q}\mathbf{n}d(\partial w)=\int_{w}div\mathbf{q}dw $$      (7.4.1)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle Q_{2}=\int_{w}f(\mathbf{x},t)dw $$      (7.4.2)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle Q_{3}=\int_{w}\rho c\frac{\partial u}{\partial t}dw $$      (7.4.3)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle Q_{1}+Q_{2}=Q_{3} $$      (7.4.4)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle \mathbf{q}=-\mathbf{K}grad(u) $$      (7.4.5)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle div(\mathbf{K}grad(u))+f=\rho c\frac{\partial u}{\partial t} $$ (7.4.6)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

Boundry Conditions

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

$$  \displaystyle u(\mathbf{x},t)]_{\Gamma _{g}}=g(\mathbf{x},t) $$      (7.4.7)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle (q.n)]_{\Gamma _{h}}=h(\mathbf{x},t) $$      (7.4.8)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle (q.n)]_{\Gamma _{H}}=H(u-u_{\infty }) $$      (7.4.9)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

WRF

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

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

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

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

$$  \displaystyle K=-\int_{\Gamma _{g}}wnqd\Gamma _{g}-\int_{\Gamma _{h}}wnqd\Gamma _{h}-\int_{\Gamma _{H}}wnqd\Gamma _{H}+\int_{\Omega }grad(w).K.grad(u)d\Omega $$      (7.4.12)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

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

$$  \displaystyle \tilde{m}(w,u)=\int_{\Omega}w\rho c\frac{\partial u}{\partial t}d\Omega $$      (7.4.13)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \tilde{k} (w,u)=\int_{\Omega}grad(w)Kdiv(u)d\Omega+\int_{\Gamma_{H}}wHud\Gamma_{H} $$      (7.4.14)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

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

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

Discrete WF

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

$$  \displaystyle w^{app}=\sum_{I=1}^{n}N(x,y)_{I}c_{I} $$      (7.4.16)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle u^{app}=\sum_{J=1}^{n}N(x,y)_{J}d_{J} $$      (7.4.17)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

$$  \displaystyle grad(u)=grad(N(x,y)_{I})d_{I}=B(x,y)_{I}d_{I} $$      (7.4.18)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

$$  \displaystyle \tilde{m}(w,u)+\tilde{k}(w,u)=\tilde{f}(w,u) $$      (7.4.19) Now lets substitute everything in the equation above.
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

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

 * }

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

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

\end{bmatrix} $$     (7.4.21)
 *

 * }

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

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

\end{bmatrix} $$     (7.4.22)
 *

 * }

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

$$  \displaystyle \tilde{\mathbf{K}}=\begin{bmatrix} K_{E} &K_{EF} \\ K_{EF}^{T} & K_{F} \end{bmatrix} $$     (7.4.23)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

Example of calculation of elements of K & M matrices when only one node lie on essential boundary.
 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} & K_{E}=\int_{\Omega }B_{1}.K.B_{1}d\Omega+\int\limits_ {{N_1}H{N_1}d{\Gamma _H}} &\\ & K_{EF}=\int_{\Omega }B_{1}.K.B_{J}d\Omega+\int\limits_ {{N_1}H{N_j}d{\Gamma _H}} & J=2,n\\ & K_{FF}=\int_{\Omega }B_{I}.K.B_{J}d\Omega+\int\limits_ {{N_i}H{N_j}d{\Gamma _H}} & I,J=2,n
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$     (7.4.24)
 *

 * }

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

$$  \displaystyle \tilde{\mathbf{M}}=\begin{bmatrix} M_{E} &M_{EF} \\ M_{EF}^{T} & M_{F} \end{bmatrix} $$     (7.4.25)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

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

\end{align} $$     (7.4.26)
 *

 * }

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

$$ \underline {\tilde c} = \left[ {\begin{array}{ccccccccccccccc} {\underline } \\ \hline \end{array}} \right] $$     (7.4.27)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

$$  \displaystyle \underline {\tilde c} = \left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right] $$     (7.4.28)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

Structure of d matrix
 * {| style="width:100%" border="0"

$$  \displaystyle \underline {\tilde d} = \left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right] $$     (7.4.29)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

Structure of K matrix for this type of problem
 * {| style="width:100%" border="0"

$$  \displaystyle {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T{\left[ {\begin{array}{ccccccccccccccc} &\vline & && \\ \hline &\vline & && \\ &\vline & && \\ &\vline & && \end{array}} \right]_{\tilde n \times \tilde n}}{\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{\tilde n \times 1}}}+{\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T{\left[ {\begin{array}{ccccccccccccccc} {}&\vline & {}&{}&{} \\ \hline {}&\vline & {}&{}&{} \\ {}&\vline & {}&{\underline K _{HH}^H}&{} \\ {}&\vline & {}&{}&{} \end{array}} \right]_{\tilde n \times \tilde n}}{\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{\tilde n \times 1}}} $$     (7.4.30)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

Components of force matrix
 * {| style="width:100%" border="0"

$$  \displaystyle {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T\left[ {\begin{array}{ccccccccccccccc} {\underline F _E^f} \\ \hline {\underline F _h^f} \\ {\underline F _H^f} \\ {\underline F _R^f} \end{array}} \right] + {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T\left[ {\begin{array}{ccccccccccccccc} {} \\ \hline {\underline F _h^h} \\ {} \\  {} \end{array}} \right] + {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T\left[ {\begin{array}{ccccccccccccccc} {} \\ \hline {} \\  {\underline F _H^H} \\ {} \end{array}} \right] $$     (7.4.31)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

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

 * }

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

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

 * }

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

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

 * }

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

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

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

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

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

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

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

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

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

1. Static case (steady state) $$f\left(\boldsymbol{x} \right)=1$$ in $$\Omega =\square $$

1.(a) using 2D LIBF.

1(b). 2D LLEBF

1b1). uniform mesh

1b2). non-uniform mesh

2. Dynamic Case(Transient)

$$\rho c=3$$

initial condition: $$u\left( x,t=0 \right)=xy$$, $$\forall x\in \Omega $$

2(a). using 2D LIBF.

$$f\left(\boldsymbol{x},t \right)=0$$, $$\forall \boldsymbol{x}\in \Omega $$, $$\forall t>0$$

$$g\left( \boldsymbol{x},t \right)=2$$ on $${{\Gamma }_{g}}=\partial \Omega $$

2(b).$$f\left( \boldsymbol{x},t \right)=1$$, $$\forall \boldsymbol{x}\in \Omega $$, $$\forall t>0$$

$$g\left( \boldsymbol{x},t \right)=2$$ on $${{\Gamma }_{g}}=\partial \Omega $$

2b1) 2D LIBF

2b2) 2D LLEBF

2b2a) uniform mesh

2b2a) non uniform mesh

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

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

For LLEBF we will have only four basis function for each element, which can be constructed from the coordinates of that particular element in the same way as LIBF are constructed.

title: User:Eml5526.s11.team2.sandhu/hw7


text: File:Sui blackbox.png

title: File:Sui blackbox.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.fig1.png