User:Eml5526.s11.team7/HW5

Problem Description

 * 1) Solve the general 1-D model for an elastic bar with data set 1b (G1DM1.0/D1b) using the weight function with polynomial, Fourier, and exponential basis functions until convergence of $$\displaystyle u^h (0.5) \rightarrow O(10^{-6}) $$.
 * 2) Plot u, uh, and convergence (error versus n).

Given
For G1DM1.0/D1 and D1b, the governing partial differential equation (PDE) and given conditions are,


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

$$\displaystyle \frac{\partial}{\partial x} \left[ a_2(x) \frac{\partial u}{\partial x} \right ] + f(x,t) = \bar{m} \frac{\partial ^s u}{\partial t ^s} (x,t)$$ (5.1.1)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }


 * where
 * $$\displaystyle \Omega = [{0,1}]$$


 * $$\displaystyle a_2 (x) = 2 + 3x$$


 * $$\displaystyle f(x) = 5x$$


 * $$\displaystyle \bar{m} \frac{\partial ^s u}{\partial t ^s} (x,t) = 0 $$ (static case)

The natural and essential boundary conditions are respectively defined below for data set D1b.


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

$$ n(1) a_2(1) \frac{\partial u}{\partial x}(1)=12 \qquad \Rightarrow \qquad \frac{\partial u}{\partial x}(1)= \frac{12}{5}$$ (5.1.2)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle u(0)=4$$ (5.1.3)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Weak Solution
The weak form must first be constructed such that the weight function, $$\displaystyle w(x)$$, vanishes on the essential boundary (x=0). First, the governing equation (5.1.1) and natural boundary condition (5.1.2) are multiplied by the weight function and integrated over the domain. The restriction on the weight function is that $$\displaystyle w(0)=0$$.


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

$$ \int_0^1 w(x) \left[ \frac{d}{dx} \left(a_2(x) \frac{du}{dx} \right) + f(x) \right] dx = 0 \qquad and \qquad w(x)\left( \frac{du}{dx} - \frac{12}{5} \right )_{x=1} = 0 $$     (5.1.4)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

The first equation can be expanded as


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

$$ \int_0^1 w(x) \left[ \frac{d}{dx} \left(a_2(x) \frac{du}{dx} \right) \right] dx + \int_0^1 w(x) f(x) dx= 0 \qquad \forall w. $$ (5.1.5)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Integrating by parts yields


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

$$ \int_0^1 w(x) \frac{d}{dx} \left(a_2(x) \frac{du}{dx} \right) dx = \left( w(x)a_2(x) \frac{du}{dx} \right )_{x=1} - \left( w(x)a_2(x) \frac{du}{dx} \right )_{x=0} - \int_0^1 \frac{dw}{dx}a_2(x) \frac{du}{dx} dx. $$     (5.1.6)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Recall that weight function must vanish at the essential boundary. Then by combining equations 5.1.5 and 5.1.6, the following expression is determined.


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

$$ \int_0^1 \frac{dw}{dx}a_2(x) \frac{du}{dx} dx = \left( w(x) a_2(x) \frac{du}{dx} \right)_{x=1} + \int_0^1 w(x)f(x)dx \qquad \forall w \quad with \quad w(0)=0. $$     (5.1.7)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Then substitute the natural boundary condition (equation 5.1.2) and the given conditions for $$\displaystyle a_2(x)$$ and $$\displaystyle f(x)$$ into the previous equation. The result is the weak form for G1DM1.0/D1b.


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

$$ \int_0^1 \frac{dw}{dx}(2+3x) \frac{du}{dx} dx = 12w(1) + \int_0^1 w(x)5x dx \qquad \forall w \quad with \quad w(0)=0. $$     (5.1.8)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * 
 * }

Exact Solution
The governing partial differential equation (PDE):
 * {| style="width:100%" border="0"

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

\frac{\partial }\left[ {\left( {2 + 3x} \right)\frac} \right] + 5x = 0{\text{ }}\forall {\text{x}} \in \left] {0,1} \right[

$$     (5.1.9) The natural and essential B.C.:
 * 
 * }
 * {| style="width:100%" border="0"

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

\frac\left( 1 \right) = \frac{5}

$$     (5.1.10)
 * 
 * }
 * {| style="width:100%" border="0"

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

u\left( 0 \right) = 4

$$     (5.1.11) Then
 * 
 * }
 * {| style="width:100%" border="0"

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

\int {\frac{\partial }\left[ {\left( {2 + 3x} \right)\frac} \right]} dx = \int { - 5xdx}

$$     (5.1.12)
 * 
 * }
 * {| style="width:100%" border="0"

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

\left( {2 + 3x} \right)\frac + {C_1} = - \frac{5}{2}{x^2}{\text{ where }}{C_1}{\text{ is a constant}}

$$     (5.1.13) By the natural B.C.
 * 
 * }
 * {| style="width:100%" border="0"

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

{C_1} = - \frac{2}

$$     (5.1.14) Then
 * 
 * }
 * {| style="width:100%" border="0"

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

\frac = \frac

$$     (5.1.15)
 * 
 * }
 * {| style="width:100%" border="0"

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

u = - \frac{5}{x^2} + \frac{5}{9}x + \frac\log \left( {3x + 2} \right) + {C_2}{\text{ where }}{C_2}{\text{ is a constant}}

$$     (5.1.16) By the essential B.C.
 * 
 * }
 * {| style="width:100%" border="0"

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

{C_2} = 4 - \frac\log \left( 2 \right)

$$     (5.1.17) Then
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

u = - \frac{5}{x^2} + \frac{5}{9}x + \frac\log \left( {3x + 2} \right) + 4 - \frac\log \left( 2 \right)

$$     (5.1.18)
 * <p style="text-align:right">
 * }

For Polynomial Basis
The corresponding family $${{\bar{F}}_{P}}=\left\{ {{{\bar{b}}}_{j}} \right\}$$ of polynomial basis function satisfying the Constraint Breaking Solution (CBS) for $$\Gamma_g = \left \{ 0 \right \}$$ is:

$${{F}_{P}}=\left\{ {{\left( x \right)}^{j}},\text{ }j=0,1,2,\ldots \right\}$$

Now that the family of basis functions are determined, we recall that the discrete weak form is defined by

$$ \mathbf{\tilde{d}} \cdot \mathbf{b}( \beta ) = g \qquad \text{and} \qquad \mathbf{\tilde{c}} \cdot \left[ \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} + \mathbf{\tilde{K}} \mathbf{\tilde{d}} - \mathbf{\tilde{F}} \right] = 0 $$

where the expression $$\displaystyle \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} = 0 $$ for the static case. Then the discrete weak form is simplified to:

\mathbf{K} \mathbf{d} = \mathbf{F} \qquad \Rightarrow \qquad \mathbf{K} \mathbf{d} = \mathbf{F}_F - \mathbf{K}_{FE} g. $$

The first constant, $$\displaystyle d_0$$, is determined from the formulation of the CBS. It is known that:


 * $$\displaystyle u^h(\beta) = d_0 b_0(\beta)=g$$


 * $$\displaystyle d_0= \frac{g}{b_0(\beta)}=4/1=4$$

Now the matrices used to determine the rest of the coefficients are determined as follows. Note that the number of terms needed in the approximate solution, $$\displaystyle n$$, is determined by the convergence interval given. The first few terms of each matrix computation is shown. A Matlab script was written to iterate these computation for increasing values of $$\displaystyle n$$ until the desired accuracy was reached.



\displaystyle

{F_F} = \left\{ {{F_i},i = 1,2,...,n} \right\} = \left\{ {\tilde f\left( \right),i = 1,2,...,n} \right\}

$$



\displaystyle

\tilde f({b_i}) = {b_i}(\alpha )h + \int\limits_\alpha ^\beta {{b_i}fdx}

$$



\displaystyle

h = {a_2}(\alpha )\frac = (2+3\alpha)\frac = 5*12/5=12

$$



\displaystyle

\tilde f({b_1}) = {b_1}(1)h + \int\limits_0^1 {{b_1}(5x)dx}=1*12+\int\limits_0^1  {{(x)}(5x)dx}=41/3

$$



\displaystyle

\tilde f({b_2}) = {b_2}(1)h + \int\limits_0^1 {{b_2}(5x)dx}=1*12 + \int\limits_0^1  {{(x)^2}(5x)dx}=53/4

$$

$$ \displaystyle

\therefore {F_F} = \left[ {\begin{array}{ccccccccccccccc} {41/3} \\ {53/4} \\  {.} \\  {.} \\  {.} \end{array}} \right]

$$



\displaystyle

{K_{FE}} = \left\{ {{K_{0j}};j = 1,...,n} \right\} = \left\{ {\tilde k\left( {{b_0},{b_j}} \right),j = 1,...,n} \right\}

$$



\displaystyle

\tilde k({b_i},{b_j}) = \int\limits_0^1 {b{'_i}{a_2}b{'_j}dx}

$$



\displaystyle

\tilde k({b_0},{b_1}) = \int\limits_0^1 {b{'_0}{a_2}b{'_1}dx} = \int\limits_0^1 {b{'_0}(2+3x)b{'_1}dx}  = 0 = \tilde k({b_0},{b_2}) = \tilde k({b_0},{b_3}) = 0

$$



\displaystyle

\therefore {K_{EF}} = \left[ {\begin{array}{ccccccccccccccc} 0 \\ 0 \\  . \\ .\\ . \end{array}} \right]

$$

And：

\displaystyle

\because {K_{FE}} = {K_{EF}}^T = \left[ {\begin{array}{ccccccccccccccc} 0 \\ 0 \\  . \\ .\\ . \end{array}} \right]

$$



\displaystyle

{K_{FF}} = \left[ {{K_{ij}};i,j = 1,2,...} \right] = \tilde k({b_i},{b_j})

$$



\displaystyle

\tilde k({b_1},{b_1}) = \int\limits_0^1 {b{'_1}{2+3x}b{'_1}dx} = \int\limits_0^1 {1*(2+3x)*1dx}=7/2

$$



\displaystyle

\tilde k({b_1},{b_2}) = \int\limits_0^1 {b{'_1}{2+3x}b{'_2}dx} = \int\limits_0^1 {1*(2+3x)*2(x)dx}=4

$$



\displaystyle

\tilde k({b_2},{b_1}) = \int\limits_0^1 {b{'_2}{2+3x}b{'_1}dx} = \int\limits_0^1 {2(x)*(2+3x)*1dx}=4

$$



\displaystyle

\tilde k({b_2},{b_2}) = \int\limits_0^1 {b{'_2}{2+3x}b{'_2}dx} = \int\limits_0^1 {2(x)*(2+3x)*2(x)dx}=17/3

$$

\displaystyle \therefore {K_{FF}} = \begin{bmatrix} 7/2&4 &.&. \\ 4&17/3  & .&.\\ . &.  & .&.\\ . &.  & .&. \end{bmatrix} $$

The matrices are then solved for the coefficients:

\displaystyle

Kd = F $$



\displaystyle

F = {F_F} - K_{FE}g = {F_F} $$



\displaystyle

K_{FF}d_F=F_F $$



\displaystyle

d_F=K_{FF}^{-1} F_F $$

The approximate solution, $$\displaystyle u^h(x)$$ is defined as,


 * $$\displaystyle u^h(x)= d_0b_0+\sum_{j=1}^{n}d_jb_j(x) \qquad \text{where} \quad d_0=4 \ \text{and} \ b_0=1$$.

As stated above a Matlab script was written to iterate $$\displaystyle n$$ until the desired accuracy was reached. At $$\displaystyle n=6$$ the error was $$\displaystyle  O(3.5840*10^{-6})$$. The following plots were generated:



For Fourier Basis
The family of Fourier basis functions is of the form
 * $$\displaystyle F_F = \left\{ 1, cos(x-1), sin(x), cos(2x)-1, sin(2x),..., cos\left( \frac{nx}{2}\right)-1, sin\left(\frac{nx}{2}\right) \right\}$$.

Now that the family of basis functions are determined, we recall that the discrete weak form is defined by


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

$$ \mathbf{\tilde{d}} \cdot \mathbf{b}( \beta ) = g \qquad \text{and} \qquad \mathbf{\tilde{c}} \cdot \left[ \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} + \mathbf{\tilde{K}} \mathbf{\tilde{d}} - \mathbf{\tilde{F}} \right] = 0 $$
 * style="width:10%; |
 * style="width:10%; |
 * }


 * where the expression $$\displaystyle \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} = 0 $$ for the static case. Then equation 5.1.11 for the unconstrained case can be rewritten as


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

$$ \mathbf{K} \mathbf{d} = \mathbf{F} \qquad \Rightarrow \qquad \mathbf{K} \mathbf{d} = \mathbf{F}_F - \mathbf{K}_{FE} g. $$
 * style="width:10%; |
 * style="width:10%; |
 * }

This agrees with Vu-Quoc's notes. Now the matrices must be determined. From the derivation of the polynomial basis case above, it holds that $$\displaystyle d_0 = 4 $$, and $$\displaystyle \mathbf{d}_F = \mathbf{K}_{FF}^{-1} \mathbf{F}_F$$.


 * For n=6,


 * $$\displaystyle

\mathbf{F}_F = \tilde f(b_i) = {b_i}(1)h + \int_0^1 {{b_i}(5x)dx} $$


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

$$ \mathbf{F}_F = \begin{bmatrix} \tilde f(b_1) \\ \tilde f(b_2) \\ \tilde f(b_3) \\ \tilde f(b_4) \\ \tilde f(b_5) \\ \tilde f(b_6) \end{bmatrix} \qquad \Rightarrow \qquad \mathbf{F}_F = \begin{bmatrix} -6.108 \\ 11.603 \\ -18.991 \\ -13.089 \\ -27.25 \\ 3.422 \end{bmatrix} $$
 * style="width:10%; |
 * style="width:10%; |
 * }

The stiffness matrix is determined,


 * $$\displaystyle

\mathbf{K}_{FF} = \tilde k(b_i,b_j) = \int_0^1 {b_i}'(2+3x){b_j}'dx $$


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

$$ \mathbf{K}_{FF} = \begin{bmatrix} 1.144 & -1.361 & 3.256 & -0.541 & 3.701 & 2.321 \\ -1.361 & 2.356 & -4.187 & 2.4 & -5.846 & -0.067 \\ 3.256 & -4.187 & 9.512 & -2.35 & 11.619 & 5.465 \\ -0.541 & 2.4 & -2.35 & 4.488 & -5.473 & 5.35 \\ 3.701 & -5.846 & 11.619 & -5.473 & 16.813 & 2.205 \\ 2.321 & -0.067 & 5.465 & 5.35 & 2.205 & 14.687 \end{bmatrix} $$
 * style="width:10%; |
 * style="width:10%; |
 * }

Now, the matrix $$\displaystyle \mathbf{d}_F$$ can be found.


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

$$ \mathbf{d}_{F} = \mathbf{K}_{FF}^{-1} \mathbf{F}_{F} \qquad \Rightarrow \qquad \mathbf{d}_{F} = \begin{bmatrix} 17.629 \\ 23.592 \\ 0.383 \\ -9.94 \\ -0.952 \\ 1.177 \end{bmatrix} $$
 * style="width:10%; |
 * style="width:10%; |
 * }

The approximate solution, $$\displaystyle u^h(x)$$ is defined as,


 * $$\displaystyle u^h(x)= d_0b_0+\sum_{j=1}^{n}d_jb_j(x) \qquad \text{where} \quad d_0=4 \ \text{and} \ b_0=1$$.


 * OR (simplified)


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

$$\displaystyle u^h(x)= 17.629*cos(x) \ + \ 23.592*sin(x) \ + \ 0.383*cos(2x) \ - \ 9.94*sin(2x) \ - \ 0.952*cos(3x) \ + \ 1.177*sin(3x) \ - \ 13.06 $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }

For Exponential Basis
The exponential basis begins to converge as n increases to around n=3, but then quickly begins to diverge. At n=10, the solution begins to diverge at the right end, and when n=20, it is clear that the solution is not stable.

Problem Description
Rework the previous problem for data set 1, i.e. G1DM1.0/D1.
 * Note: The difference between data set 1 and data set 1b is where the boundary conditions are applied. Specifically, the essential boundary condition for D1 is applied at x=1 while the natural boundary condition is applied at x=0. The opposite is true for data set 1b.

The essential and natural boundary conditions are respectively defined below for data set D1.
 * $$ \Gamma_g = \left \{ 1 \right \} \quad \text{and} \quad \Gamma_h = \left \{ 0 \right \}$$


 * $$ n(0) a_2(0) \frac{\partial u}{\partial x}(0)=-12 \qquad \Rightarrow \qquad \frac{\partial u}{\partial x}(0)= 6$$
 * $$\displaystyle u(1)=4$$

Weak Solution
The weak form must first be constructed such that the weight function, $$\displaystyle w(x)$$, vanishes on the essential boundary (x=0). First, the governing equation (5.1.1) and natural boundary condition (5.1.2) are multiplied by the weight function and integrated over the domain. The restriction on the weight function is that $$\displaystyle w(1)=0$$.


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

$$ \int_0^1 w(x) \left[ \frac{d}{dx} \left(a_2(x) \frac{du}{dx} \right) + f(x) \right] dx = 0 \qquad and \qquad w(x)\left( \frac{du}{dx} - \frac{12}{5} \right )_{x=0} = 0 $$     (5.2.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The first equation can be expanded as


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

$$ \int_0^1 w(x) \left[ \frac{d}{dx} \left(a_2(x) \frac{du}{dx} \right) \right] dx + \int_0^1 w(x) f(x) dx= 0 \qquad \forall w. $$ (5.2.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Integrating by parts yields


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

$$ \int_0^1 w(x) \frac{d}{dx} \left(a_2(x) \frac{du}{dx} \right) dx = \left( w(x)a_2(x) \frac{du}{dx} \right )_{x=1} - \left( w(x)a_2(x) \frac{du}{dx} \right )_{x=0} - \int_0^1 \frac{dw}{dx}a_2(x) \frac{du}{dx} dx. $$     (5.2.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Recall that weight function must vanish at the essential boundary. Then by combining equations 5.1.5 and 5.1.6, the following expression is determined.


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

$$ \int_0^1 \frac{dw}{dx}a_2(x) \frac{du}{dx} dx =- \left( w(x) a_2(x) \frac{du}{dx} \right)_{x=0} + \int_0^1 w(x)f(x)dx \qquad \forall w \quad with \quad w(0)=0. $$     (5.2.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Then substitute the natural boundary condition (equation 5.1.2) and the given conditions for $$\displaystyle a_2(x)$$ and $$\displaystyle f(x)$$ into the previous equation. The result is the weak form for G1DM1.0/D1.


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

$$ \int_0^1 \frac{dw}{dx}(2+3x) \frac{du}{dx} dx = -12w(0) + \int_0^1 w(x)5x dx \qquad \forall w \quad with \quad w(0)=0. $$     (5.2.5)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Exact Solution
The governing partial differential equation (PDE):
 * {| style="width:100%" border="0"

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

\frac{\partial }\left[ {\left( {2 + 3x} \right)\frac} \right] + 5x = 0{\text{ }}\forall {\text{x}} \in \left] {0,1} \right[

$$     (5.2.6) The natural and essential B.C.:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\frac\left( 0 \right) = 6

$$     (5.2.7)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

u\left( 1 \right) = 4

$$     (5.2.8) Then
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\int {\frac{\partial }\left[ {\left( {2 + 3x} \right)\frac} \right]} dx = \int { - 5xdx}

$$     (5.2.9)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\left( {2 + 3x} \right)\frac + {C_1} = - \frac{5}{2}{x^2}{\text{ where }}{C_1}{\text{ is a constant}}

$$     (5.2.10) By the natural B.C.
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{C_1} = - 12

$$     (5.2.11) Then
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\frac = \frac

$$     (5.2.12)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

u = - \frac{5}{x^2} + \frac{5}{9}x - \frac\log \left( {3x + 2} \right) + \frac{5}{9} + {C_2}{\text{ where }}{C_2}{\text{ is a constant}}

$$     (5.2.13) By the essential B.C.
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{C_2} = 4 + \frac\log \left( 5 \right)-\frac{10}{9}+\frac{5}{12]}

$$     (5.2.14) Then
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

u = - \frac{5}{x^2} + \frac{5}{9}x - \frac\log \left( {3x + 2} \right) + \frac{5}{12} + \frac\log \left( 5 \right)+4-\frac{5}{9}

$$     (5.2.15)
 * <p style="text-align:right">
 * }

For Polynomial Basis
The corresponding family $${{\bar{F}}_{P}}=\left\{ {{{\bar{b}}}_{j}} \right\}$$ of polynomial basis function satisfying the Constraint Breaking Solution (CBS) for $$\Gamma_g = \left \{ 1 \right \}$$ is:

$${{F}_{P}}=\left\{ {{\left( x-1 \right)}^{j}},\text{ }j=0,1,2,\ldots \right\}$$

Now that the family of basis functions are determined, we recall that the discrete weak form is defined by

$$ \mathbf{\tilde{d}} \cdot \mathbf{b}( \beta ) = g \qquad \text{and} \qquad \mathbf{\tilde{c}} \cdot \left[ \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} + \mathbf{\tilde{K}} \mathbf{\tilde{d}} - \mathbf{\tilde{F}} \right] = 0 $$

where the expression $$\displaystyle \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} = 0 $$ for the static case. Then the discrete weak form is simplified to:

\mathbf{K} \mathbf{d} = \mathbf{F} \qquad \Rightarrow \qquad \mathbf{K} \mathbf{d} = \mathbf{F}_F - \mathbf{K}_{FE} g. $$

The first constant, $$\displaystyle d_0$$, is determined from the formulation of the CBS. It is known that:


 * $$\displaystyle u^h(\beta) = d_0 b_0(\beta)=g$$


 * $$\displaystyle d_0= \frac{g}{b_0(\beta)}=4/1=4$$

Now the matrices used to determine the rest of the coefficients are determined as follows. Note that the number of terms needed in the approximate solution, $$\displaystyle n$$, is determined by the convergence interval given. The first few terms of each matrix computation is shown. A Matlab script was written to iterate these computation for increasing values of $$\displaystyle n$$ until the desired accuracy was reached.



\displaystyle

{F_F} = \left\{ {{F_i},i = 1,2,...,n} \right\} = \left\{ {\tilde f\left( \right),i = 1,2,...,n} \right\}

$$



\displaystyle

\tilde f({b_i}) = {b_i}(\alpha )h + \int\limits_\alpha ^\beta {{b_i}fdx}

$$



\displaystyle

h = {a_2}(\alpha )\frac = (2+3\alpha)\frac = 2*6=12

$$



\displaystyle

\tilde f({b_1}) = {b_1}(0)h + \int\limits_0^1 {{b_1}(5x)dx}=1*12+\int\limits_0^1  {{(x-1)}(5x)dx}=-77/6

$$



\displaystyle

\tilde f({b_2}) = {b_2}(0)h + \int\limits_0^1 {{b_2}(5x)dx}=-1*12 + \int\limits_0^1  {{(x-1)^2}(5x)dx}=149/12

$$

$$ \displaystyle

\therefore {F_F} = \left[ {\begin{array}{ccccccccccccccc} {-77/6} \\ {149/12} \\  {.} \\  {.} \\  {.} \end{array}} \right]

$$



\displaystyle

{K_{FE}} = \left\{ {{K_{0j}};j = 1,...,n} \right\} = \left\{ {\tilde k\left( {{b_0},{b_j}} \right),j = 1,...,n} \right\}

$$



\displaystyle

\tilde k({b_i},{b_j}) = \int\limits_0^1 {b{'_i}{a_2}b{'_j}dx}

$$



\displaystyle

\tilde k({b_0},{b_1}) = \int\limits_0^1 {b{'_0}{a_2}b{'_1}dx} = \int\limits_0^1 {b{'_0}(2+3x)b{'_1}dx}  = 0 = \tilde k({b_0},{b_2}) = \tilde k({b_0},{b_3}) = 0

$$



\displaystyle

\therefore {K_{EF}} = \left[ {\begin{array}{ccccccccccccccc} 0 \\ 0 \\  . \\ .\\ . \end{array}} \right]

$$

And：

\displaystyle

\because {K_{FE}} = {K_{EF}}^T = \left[ {\begin{array}{ccccccccccccccc} 0 \\ 0 \\  . \\ .\\ . \end{array}} \right]

$$



\displaystyle

{K_{FF}} = \left[ {{K_{ij}};i,j = 1,2,...} \right] = \tilde k({b_i},{b_j})

$$



\displaystyle

\tilde k({b_1},{b_1}) = \int\limits_0^1 {b{'_1}{2+3x}b{'_1}dx} = \int\limits_0^1 {1*(2+3x)*1dx}=7/2

$$



\displaystyle

\tilde k({b_1},{b_2}) = \int\limits_0^1 {b{'_1}{2+3x}b{'_2}dx} = \int\limits_0^1 {1*(2+3x)*2(x-1)dx}=-3

$$



\displaystyle

\tilde k({b_2},{b_1}) = \int\limits_0^1 {b{'_2}{2+3x}b{'_1}dx} = \int\limits_0^1 {2(x-1)*(2+3x)*1dx}=-3

$$



\displaystyle

\tilde k({b_2},{b_2}) = \int\limits_0^1 {b{'_2}{2+3x}b{'_2}dx} = \int\limits_0^1 {2(x-1)*(2+3x)*2(x-1)dx}=11/3

$$

\displaystyle \therefore {K_{FF}} = \begin{bmatrix} 7/2&-3 &.&. \\ -3&11/3  & .&.\\ . &.  & .&.\\ . &.  & .&. \end{bmatrix} $$

The matrices are then solved for the coefficients:

\displaystyle

Kd = F $$



\displaystyle

F = {F_F} - K_{FE}g = {F_F} $$



\displaystyle

K_{FF}d_F=F_F $$



\displaystyle

d_F=K_{FF}^{-1} F_F $$

The approximate solution, $$\displaystyle u^h(x)$$ is defined as,


 * $$\displaystyle u^h(x)= d_0b_0+\sum_{j=1}^{n}d_jb_j(x) \qquad \text{where} \quad d_0=4 \ \text{and} \ b_0=1$$.

As stated above a Matlab script was written to iterate $$\displaystyle n$$ until the desired accuracy was reached. At $$\displaystyle n=6$$ the error was $$\displaystyle  O(3.4091*10^{-6})$$. The following plots were generated:



For Fourier Basis
The family of Fourier basis functions is of the form
 * $$\displaystyle F_F = \left\{ 1, cos(x-1)-1, sin(x-1), cos[2(x-1)]-1, sin[2(x-1)],..., cos[\frac{n}{2}(x-1)]-1, sin[\frac{n}{2}(x-1) \right\}$$.

Now that the family of basis functions are determined, we recall that the discrete weak form is defined by


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

$$ \mathbf{\tilde{d}} \cdot \mathbf{b}( \beta ) = g \qquad \text{and} \qquad \mathbf{\tilde{c}} \cdot \left[ \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} + \mathbf{\tilde{K}} \mathbf{\tilde{d}} - \mathbf{\tilde{F}} \right] = 0 $$
 * style="width:10%; |
 * style="width:10%; |
 * }


 * where the expression $$\displaystyle \mathbf{\tilde{M}} \mathbf{\tilde{d}}^{(s)} = 0 $$ for the static case. Then equation 5.1.11 for the unconstrained case can be rewritten as


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

$$ \mathbf{K} \mathbf{d} = \mathbf{F} \qquad \Rightarrow \qquad \mathbf{K} \mathbf{d} = \mathbf{F}_F - \mathbf{K}_{FE} g. $$
 * style="width:10%; |
 * style="width:10%; |
 * }

This agrees with Vu-Quoc's notes. Now the matrices must be determined. From the derivation of the polynomial basis case above, it holds that $$\displaystyle d_0 = 4 $$, and $$\displaystyle \mathbf{d}_F = \mathbf{K}_{FF}^{-1} \mathbf{F}_F$$.


 * For n=6,


 * $$\displaystyle

\mathbf{F}_F = \tilde f(b_i) = {b_i}(0)h + \int_0^1 {{b_i}(5x)dx} $$


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

$$ \mathbf{F}_F = \begin{bmatrix} \tilde f(b_1) \\ \tilde f(b_2) \\ \tilde f(b_3) \\ \tilde f(b_4) \\ \tilde f(b_5) \\ \tilde f(b_6) \end{bmatrix} \qquad \Rightarrow \qquad \mathbf{F}_F = \begin{bmatrix} -5.718 \\ -10.89 \\ -17.724 \\ -12.275 \\ -25.274 \\ -3.282 \end{bmatrix} $$
 * style="width:10%; |
 * style="width:10%; |
 * }

The stiffness matrix is determined,


 * $$\displaystyle

\mathbf{K}_{FF} = \tilde k(b_i,b_j) = \int_0^1 {b_i}'(2+3x){b_j}'dx $$


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

$$ \mathbf{K}_{FF} = \begin{bmatrix} 0.764 & 1.117 & 2.305 & 0.885	& 3.06	& -0.773 \\ 1.117 & 2.736 & 3.675 & 3.819 & 5.929 & 2.855 \\ 2.305 & 3.675 & 7.137 & 3.437 & 10.079 & -1.18 \\ 0.885 & 3.819 & 3.437 & 6.863 & 7.19 & 8.293 \\ 3.06 & 5.929 & 10.079 & 7.19 & 16.154 & 2.415 \\ -0.773 & 2.855 & -1.18 & 8.293 & 2.415 & 15.346 \end{bmatrix} $$
 * style="width:10%; |
 * style="width:10%; |
 * }

Now, the matrix $$\displaystyle \mathbf{d}_F$$ can be found.


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

$$ \mathbf{d}_{F} = \mathbf{K}_{FF}^{-1} \mathbf{F}_{F} \qquad \Rightarrow \qquad \mathbf{d}_{F} = \begin{bmatrix} -26.379 \\ 1.214 \\ 8.796 \\ -3.542 \\ -1.074 \\ 0.991 \end{bmatrix} $$
 * style="width:10%; |
 * style="width:10%; |
 * }

The approximate solution, $$\displaystyle u^h(x)$$ is defined as,


 * $$\displaystyle u^h(x)= d_0b_0+\sum_{j=1}^{n}d_jb_j(x) \qquad \text{where} \quad d_0=4 \ \text{and} \ b_0=1$$.


 * OR (simplified)


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

$$\displaystyle u^h(x)= - 26.379*cos(x-1) \ + \ 1.214*sin(x-1) \ + \ 8.796*cos(2x-2) \ - \ 3.542*sin(2x-2)] \ - \ 1.074*cos(3x-3) \ + \ 0.991*sin(3x-3) \ + \ 22.657 $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }

For Exponential Basis
The exponential basis function is given by:
 * {| style="width:100%" border="0"

$$ \displaystyle {u^{h}}(x)={d_{j}}\left (e^{j(x-1)}-1 \right ) $$       (5.2.exp.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \frac{{u^{h}}(x)}{dx} ={d_{j}}\left (e^{j(x-1)} \right ) $$       (5.2.exp.2) The general form of the stiffness matrix becomes:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle K=\int_{0}^{1}i*e^{i(x-1)}*(2+3x)*j*e^{j(x-1)}dx $$       (5.2.exp.3) The general form of the force matrix becomes:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle F=12(e^{i(0-1)}-1)+\int_{0}^{1} 5x*i(e^{i(x-1)}-1)dx $$       (5.2.exp.4) Using MatLab, the K and F matrix with N=10 is:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle K=\begin{bmatrix} 1.716166179 &	2.633475288	& 3.170329089	& 3.512453499	&3.74690156	&3.917139098	&4.046398014	&4.147965319	&4.229930538	&4.297494435\\ 2.633475288	&4.227105451	&5.268680249	&5.995042496	&6.528565163	&6.93668231	&7.258939308	&7.519876512	&7.735489984	&7.916648746\\ 3.170329089	&5.268680249	&6.744422808	&7.834278196	&8.670852887	&9.332921967	&9.869837922	&10.31398664	&10.68747581	&11.00590793\\ 3.512453499	&5.995042496	&7.834278196	&9.248909746	&10.3699133	&11.27981477	&12.03298442	&12.66663799	&13.20708952	&13.67346515\\ 3.74690156	&6.528565163	&8.670852887	&10.3699133	&11.74980705	&12.89248331	&13.85413531	&14.67454391	&15.38264829	&15.99999816\\ 3.917139098	&6.93668231	&9.332921967	&11.27981477	&12.89248331	&14.24996774	&15.4082711	&16.40815817	&17.27999802	&18.04687424\\ 4.046398014	&7.258939308	&9.869837922	&12.03298442	&13.85413531	&15.4082711	&16.7499948	&17.91999794	&18.94921795	&19.86159138\\ 4.147965319	&7.519876512	&10.31398664	&12.66663799	&14.67454391	&16.40815817	&17.91999794	&19.24999918	&20.42906542	&21.48148136\\ 4.229930538	&7.735489984	&10.68747581	&13.20708952	&15.38264829	&17.27999802	&18.94921795	&20.42906542	&21.74999987	&22.93628804 \end{bmatrix} $$       (5.2.exp.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle F=\begin{bmatrix} -194.4238\\ 556.4552\\ -794.5674\\ 388.5998\\ 355.4029\\ -550.4317\\ 202.6984\\ 13.1421\\ 77.9578\\ -68.7552\\ -267.2934\\ 294.8667\\ 224.2897\\ -420.1257 150.8844 \end{bmatrix} $$       (5.2.exp.5) Solving for the general solution with n=6 gives an error of 2.83E-06 at u(.5) and a general solution of:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle u^{h}=-167.2599*e^{1*(x-1)}	487.6591e^{2*(x-1)}	-808.7098e^{3*(x-1)}	770.2175e^{4*(x-1)}	-393.7902e^{5*(x-1)}	83.8821e^{6*(x-1)}+32.0012 $$       (5.2.exp.5) The boundary condition u(1)=4 is kept and the exact solution is very close to the trial solution:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

1. Explain how the basis functions satisfy the constraint breaking solution (CBS).
The following are the basis functions:


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

$$\displaystyle \sum_{i=1}^m \prod_{k=1,i\ne k}^m \frac{(x-x_k)}{(x_i-x_k)} {m=2,4,6,8....}$$ (5.3.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

3. Solve the following equation using these basis functions and matlab.

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

$$\displaystyle \frac{\partial}{\partial x} \left[ a_2(x) \frac{\partial u}{\partial x} \right ] + f(x,t) = \bar{m} \frac{\partial ^s u}{\partial t ^s} (x,t)$$ (5.3.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

with boundary conditions:


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

$$ n(1) a_2(1) \frac{\partial u}{\partial x}(1)=12 \qquad \Rightarrow \qquad \frac{\partial u}{\partial x}(1)= \frac{12}{5}$$ (5.3.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$\displaystyle u(0)=4$$ (5.3.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

4. Make Plots of u^h vs u and u^h(.5)-u(.5) vs m

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

$$\displaystyle u_m^h$$ vs $$\displaystyle u$$ (5.3.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

and,


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

$$\displaystyle u_m^h(0.5)-u(0.5)$$ vs $$\displaystyle m$$ (5.3.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Part 1
These basis functions have the property where at each node there is a single function equal to one and the rest are zero, as seen in figure on part 2. So, all that needs to be done is set the coefficient of the single nonzero basis function to zero to meet the constraint breaking solution requirement.

Part 2
Using the following Matlab code for m=2:

With the Following Plot:



For M=4:

With Plot:



For m=6:

With Plot:



For m=8:

With Plot for m=8:



To get a better picture of the properties of these 8 basis functions, the following is a plot where the basis functions are only evaluated at the nodes:



Part 3
Starting with the weak form of equation 5.3.3 (see derivation of weak form in problem 5.1):
 * {| style="width:100%" border="0"

$$ \int_0^1 \frac{dw}{dx}(2+3x) \frac{du}{dx} dx = 12w(1) + \int_0^1 w(x)5x dx \qquad \forall w \quad with \quad w(0)=0. $$     (5.1.8)
 * <p style="text-align:right">
 * }

Now approximate the solution by:
 * {| style="width:100%" border="0"

$$\displaystyle u\simeq u^h = \sum_{i=1}^md_i\prod_{k=1,k\neq i}^m \frac{(x-x_k)}{(x_i-x_k)}=\sum_{i=1}^md_i*b_i(x) $$     (5.3.7)
 * <p style="text-align:right">
 * }

and approximate the weighting function by:


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

$$\displaystyle w\simeq w^h = \sum_{i=1}^mc_i\prod_{k=1,k\neq i}^m \frac{(x-x_k)}{(x_i-x_k)}\sum_{i=1}^mc_i*b_i(x) $$     (5.3.8)
 * <p style="text-align:right">
 * }

Where bi(x) has been substituted in for the basis functions for ease of finding the solution.

Substituting Equation 5.3.7 and 5.3.8 into 5.1.8 and using linearity of the integral operator gives after some rearrangement:


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

$$\displaystyle \sum_{i=1}^m[\int_0^1(\frac{d}{dx}c_ib_i(x))(2+3x)(\frac{d}{dx}d_ib_i(x)-12c_ib_i(1)-\int_0^1c_ib_i(x)5xdx] $$     (5.3.9)
 * <p style="text-align:right">
 * }

Factoring out the constant ci:


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

$$\displaystyle \sum_{i=1}^mc_i[\int_0^1\frac{d}{dx}b_i(x)(2+3x)(\frac{d}{dx}d_ib_i(x)-12b_i(1)-\int_0^1b_i(x)5xdx] $$     (5.3.10)
 * <p style="text-align:right">
 * }

To solve for the unknown coefficients simultaneously, we will move to write 5.3.10 in tensor form.


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

$$\displaystyle \vec C=\begin{bmatrix} \vec C_E \\ \vec C_F\end{bmatrix}=\begin{bmatrix} 0 \\ \vec C_F\end{bmatrix} $$     (5.3.11)
 * <p style="text-align:right">
 * }

The zero in equation comes from satisfying the constraint breaking solution where the coefficient at the natural boundary condition for the only nonzero basis function coefficient on the weighting function is zero. Hence the subscript E in equation 5.3.11 for essential. Looking at the coefficients for uh:


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

$$\displaystyle \vec d=\begin{bmatrix} \vec d_E \\ \vec d_F\end{bmatrix}=\begin{bmatrix} g \\ \vec d_F\end{bmatrix} $$     (5.3.12)
 * <p style="text-align:right">
 * }

Likewise here, g is the value at the essential boundary condition. Since there is only one nonzero basis function here, the coefficient (since the basis function is one) must be equal to the essential boundary condition value. The subscript F can be thought of as free and will be determined.

Now defining the following


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

$$\displaystyle \tilde k(b_i,b_j)=\int_0^1(\frac{d}{dx}b_i)(2+3x)(\frac{d}{dx}b_j)dx $$     (5.3.13)
 * <p style="text-align:right">
 * }

This shows how the far left integral term of equation 5.3.10 can be written in components of a matrix. Additionally defining:


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

$$\displaystyle K_{EE}=[\tilde k(b_1,b_1)] $$     (5.3.14)
 * <p style="text-align:right">
 * }

Where subscript one signifies the first basis functions.
 * {| style="width:100%" border="0"

$$\displaystyle K_{EF}=[\tilde k(b_1,b_j)] $$     (5.3.15)
 * <p style="text-align:right">
 * }


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

$$\displaystyle K_{FE}=[\tilde k(b_i,b_1)] $$     (5.3.16)
 * <p style="text-align:right">
 * }


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

$$\displaystyle K_{FF}=[\tilde k(b_{remaining i},b_{remaining j})] $$     (5.3.17)
 * <p style="text-align:right">
 * }

Giving:


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

$$\displaystyle \mathbf K=\begin{bmatrix} \mathbf K_{EE} \mathbf K_{EF} \\ \mathbf K_{FE} \mathbf K_{FF}\end{bmatrix} $$     (5.3.18)
 * <p style="text-align:right">
 * }

Lastly, combining the constants into:
 * {| style="width:100%" border="0"

$$\displaystyle \mathbf F = [\tilde f_i] $$     (5.3.19)
 * <p style="text-align:right">
 * }

Where,


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

$$\displaystyle \tilde f_i=12b_i(1)+\int_0^1b_i(x)5xdx $$     (5.3.20)
 * <p style="text-align:right">
 * }

Putting it all together:


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

$$\displaystyle \vec C \cdot \mathbf K\mathbf d=\begin{bmatrix}0\\\vec C_F\end{bmatrix}\cdot \Big (\begin{bmatrix} \mathbf K_{EE} \mathbf K_{EF} \\ \mathbf K_{FE} \mathbf K_{FF}\end{bmatrix}\begin{bmatrix}\mathbf g\\\vec d_F\end{bmatrix}-\begin{bmatrix}\mathbf F_E\\ \mathbf F_F\end{bmatrix} \Big )=0 $$     (5.3.21) Now, our weighting function is arbitrary so what is in parenthesis corresponding to CF must equal zero. Writing this out explicitly for four basis functions gives the following:
 * <p style="text-align:right">
 * }


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

$$\displaystyle \begin{bmatrix}K_{11}K_{12}K_{13}K_{14}\\K_{21}K_{22}K_{23}K_{24}\\K_{31}K_{32}K_{33}K_{34}\\K_{41}K_{42}K_{43}K_{44}\end{bmatrix}\begin{bmatrix}d_1\\d_2\\d_3\\d_4\end{bmatrix}=\begin{bmatrix}f_1+r_1\\f_2\\f_3\\f_4\end{bmatrix} $$     (5.3.22)
 * <p style="text-align:right">
 * }

In order to solve four our unknowns d2-d4, in terms of our knowns, we rearrange to get the following:


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

$$\displaystyle \begin{bmatrix}K_{22}K_{23}K_{24}\\K_{32}K_{33}K_{34}\\K_{42}K_{43}K_{44}\end{bmatrix}\begin{bmatrix}d_2\\d_3\\d_4\end{bmatrix}=\begin{bmatrix}f_2-K_{21}d_1\\f_3-K_{31}d_1\\f_4-K_{41}d_1\end{bmatrix} $$     (5.3.23)
 * <p style="text-align:right">
 * }

This solution method was applied in Matlab to give the following:

m=2

m=4

m=6

m=8

Part 4
The following are the resulting plots from Part 3: m=2



m=4



m=6



m=8



This Matlab code was used to plot uh-u vs m:

Giving:



Problem Description
Solve G1D1M1.D/D1d, for$$\frac{d}((2 + 3x)\frac) + 5x = 0$$, using weak form with LIBF (m=4, 6, 8, ..), until convergence of $${u^h}(0.5)to{\rm O}({10^{ - 6}})$$.
 * 1) Explain how LIBF are used as CBS
 * 2) Plot all LIBF used
 * 3) Use Matlab quad or WolfAlpha to integrate
 * 4) Plot $${u^h}$$ vs $$u$$, $${u^h}(0.5) - u(0.5)$$ vs m

Solution
(1)
 * {| style="width:100%" border="0"

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

{L_{i,m}}({x_j}) = {\delta _{i,m}} = \left\{ {\begin{array}{ccccccccccccccc} {1,forj = i} \\ {0,forj \ne i} \end{array}} \right.

$$ So for series of LIBF, at the essential boundary condition, there will be only one term does not equal to zero, the other terms equal to zero. So the CBS has been satisfied.
 * <p style="text-align:right">(Eq 5.4.1)
 * }

(2) Basis function:
 * {| style="width:100%" border="0"

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

{L_{i,m}}(x) = \prod\limits_{\begin{array}{ccccccccccccccc} {k = 1} \\ {k \ne i} \end{array}}^m {\frac}

$$     (Eq 5.4.2)
 * <p style="text-align:right">
 * }

Figs of the basis functions

m=4:



m=6:

m=8：

(3) For weak form
 * {| style="width:100%" border="0"

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

w(x) = \sum\limits_{i = 1}^m {{c_i} \cdot {L_{i,m}}}

$$     (Eq 5.4.3)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

u(x) = \sum\limits_{i = 1}^m {{b_i} \cdot {L_{i,m}}}

$$     (Eq 5.4.4) To boundary condition
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

w(1) = 0

$$     (Eq 5.4.5)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\alpha _m} = 0

$$     (Eq 5.4.6)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

u(1) = 4

$$     (Eq 5.4.7)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{b_m} = 4

$$     (Eq 5.4.8) The weak form is:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\left. {w(2 + 3x)\frac} \right|_0^1 - \int\limits_0^1 {\frac(2 + 3x)\fracdx} + \int\limits_0^1 {w \cdot 5x}  = 0

$$     (Eq 5.4.9)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

- w(0)(2 + 3*0)\frac(0) - \int\limits_0^1 {\frac(2 + 3x)\fracdx} + \int\limits_0^1 {w \cdot 5x}  = 0

$$     (Eq 5.4.10)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

- w(0)*12 + \int\limits_0^1 {w \cdot 5x} = \int\limits_0^1 {\frac(2 + 3x)\fracdx}

$$     (Eq 5.4.11)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\Rightarrow \sum\limits_i {{C_i}\left[ {\sum\limits_j  - {F_i}} \right]}  = 0

$$     (Eq 5.4.12) Where
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{K_{ij}} = \int\limits_0^1 {{L_{i,m}}'(2 + 3x){L_{j,m}}'}dx

$$     (Eq 5.4.13)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{F_i} = \int\limits_0^1 {{L_{i,m}} \cdot 5x}dx + w(0)*12

$$     (Eq 5.4.14) In this static case
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

Kd = Fwhere{F_F} - {K_{FE}}g

$$     (Eq 5.4.15) Solve for d
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{u^h} = \sum\limits_{i = 1}^m  \cdot {L_{i.m}}

$$     (Eq 5.4.16)
 * <p style="text-align:right">
 * }

Plot $${u^h}$$ vs $$u$$

m=4:



m=6:



m=8：



Calculate for error

Error vs m :



Problem 5.5.1 and 5.5.2
'''5.5.1 For the disc problem find node information : nodes numbers and their coordinates, Element numbers and their respective node numbers. 5.5.2 For the same disc problem generate 3 meshes with triangle elements by changing number of elements. ''' Eml5526.s11.team7.Srilalithkumar 16:03, 22 March 2011 (UTC)

Problem 5.5.3: CalculiX CrunchiX (CCX) for Dummies
Adopted from CCX User's Manual and Mechanical Hack's website. This walk through provides the basic overview of CCX. Experimentation and actual development of input files will be required on one's on.

Note: Previous installation of Calculix in Windows additionally installed CCX, so one can immediately begin to follow this example.

CCX is computational tool based on the information the user puts into the input file. CCX is constructed from FORTRAN77 and C programing languages, with a definite emphasis on FORTRAN. In relation to Calculix, CCX is considered the main subroutine of the total program Calculix.

CCX uses what it calls an "input deck" which specifies what will be the input. This input deck is made up of two parts, the model and step portions. Inside of these portions keywords are used. These keywords can be thought of as commands. These keywords belong either to the model or the step portion, and the model portion MUST be used before the step portion. To get a visual, look at the following figure built off of Mechanical Hack's website:



The *WORD are the keywords, or commands. The *HEADING command allows for data to be entered on the following line for the benefit of the user. The next set of commands describes the material that will be analyzed with CCX. This data can either be entered through *INCLUDE, which copies the data from the specified file (who must be in the designated path for CCX) to the input, or it can be entered explicitly as seen in the following code from the test folder on the Calculix website:



A search of the users manual provides additional methods of inputting data if needed. The remaining parts of this model section adds the necessary attributes relating to the physics of the model (boundary condition, physical properties, etc). This can be seen with the *BOUNDARY command, which allows for the proper formatting of the data being read from the input FIXED_123.BOU to be read as boundary conditions. Next the *MATERIAL is used to allow for assignment through NAME=EL of a name to the model being created. The keyword *ELASTIC allows for assignment of the elastic properties of the material. There are a multitude of additional keywords to completely specify the material of the problem. These are in the user's manual as well, and an example of them can be seen here:



Now that the material has been specified the keyword step can be used. Step allows for a quasi-static analysis. For example, in the picture a dynamic load is applied, so a step in time can be used to find the final results. End step signifies the end of the step, and there is a limit to the number of steps the program is capable of running.

With the basics of the input file created, the results are found by clicking on Tools, Calculate:



This splits the screen in two and shows the command portion of Calculix. If the input file was created correctly, the following "Job finished, Exit code:0" will be seen as well as the output.



Additionally output can be specified to a file for post processing with CGX (as seen in the figure above with the *NODE FILE and *EL FILE commands. If this happens, the data is sent to a .frb file which is post processed simply by Tools->Post Process, or Shift-F10.

Problem 5.6: Quadratic Lagrangian Element Basis Functions (QLEBF)
Plot three figures, similar to those in lecture 30, for the quadratic Lagrangian element basis functions. Consider the case where the distance between nodes is not equal.

Solution


Comments:
 * The first three quadratic basis functions are defined as:
 * $$\displaystyle b_1(x)=\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}$$,
 * $$\displaystyle b_2(x)=\frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)}$$,
 * $$\displaystyle b_3(x)=\frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}$$.
 * where $$\displaystyle x_N$$ represents the x value of the global node 'N'. Subsequent basis functions (i.e. $$\displaystyle b_4, \ b_5, \ b_6, \ \text{etc.}$$) follow the same form.


 * Each element is not discretized uniformly. In other words each node is not equidistant from adjacent nodes. The x-position for each node was arbitrarily chosen and are given as:
 * $$\displaystyle \qquad x_1=1, \qquad x_2=1.75, \qquad x_3=3, \qquad x_4=3.5, \qquad \text{and} \qquad x_5=5$$.


 * The sum of all basis functions at any point along 'x' is one.
 * The first three plots of Figure 5.6 illustrate the global representation. The last plot of Figure 5.6 shows the case for the individual element $$\displaystyle \omega^{(1)}$$ and is defined between global nodes 1 and 3. A second element, say $$\displaystyle \omega^{(2)}$$, would then be defined between global nodes 3 and 5.
 * While each element contains three nodes, each node is assigned a "local" node number from 1 to 3. Specifically for the case of element 1 in Figure 5.6, the global node numbers are the same as the local node number. However, for element 2, global node 3 corresponds to local node 1, global node 4 corresponds to local node 2, and global node 5 corresponds to local node 3.
 * Since it can easily be confusing when correlating global node numbers and local node numbers, the relationship between the two can be generalized using the equation below.


 * $$\displaystyle I = i+(3-1)(e-1) \qquad \text{(for the case of 3 nodes per element)}$$
 * where $$\displaystyle I$$ is the global node number, $$\displaystyle i$$ is the local node number, and $$\displaystyle e$$ is the element number.

Problem Description
Solve the general 1-D model for an elastic bar with data set 1b (G1DM1.0/D1b) using the Lagrangian Interpolation Basis Function until convergence of $$\displaystyle u^h (0.5) \rightarrow O(10^{-6}) $$.
 * 1. Explain how the basis functions satisfy the constraint breaking solution (CBS).
 * 2.Plot these basis functions.
 * 3.Use matlab to integrate.
 * 4.Make Plots of u^h vs u and u^h(.5)-u(.5) vs m

Solution
The Lagrangian Element Basis Function (LEBF) is given by the following equation and generates a polynomial of degree $$\displaystyle n-1$$:


 * $$L_{i,n}\left (x\right )=\prod_{k=1}^{n}\frac{x-x_{k}}{x_{i}-x_{k}},k\neq i$$

For a linear function the number of nodes is 2 ($$\displaystyle n=2$$). For uniform discretization with $$\displaystyle nel=4$$ the domain $$\displaystyle \Omega$$ will be subdivided into 4 equal length elements each having two nodes. The LEBFs for the first element $$\displaystyle n=1$$ over the elemental domain $$\displaystyle \omega^{(1)}=[0,.25]$$ can be written as follows:


 * $$L_{1,2}^{(1)}\left (x\right)=b_{1}^{(1)}=\frac{(x-x_2)}{x_1-x_2}=\frac{(x-.25)}{0-.25}=-4(x-.25)$$


 * $$L_{2,2}^{(1)}\left (x\right)=b_{2}^{(1)}=\frac{(x-x_1)}{x_2-x_1}=\frac{(x-0)}{.25-0}=4(x)$$

Recall the constraint breaking solution (CBS) states the basis functions must satisfy $$\displaystyle \left\{ b_i(x), i=0,1,2,...\right\}$$ such that $$\displaystyle b_0(\beta) \neq 0$$ and $$\displaystyle b_i(\beta) = 0$$ for $$\displaystyle i=1,2,...,n$$.

It can be seen that these basis functions for the first element satisfy the CBS at each node. For the first node, the basis function has a value of one at that node and a value of zero at every other node (there is only one other node for the linear LEBF). In other words the basis function is equal to the kronicker delta at each node.


 * $$b_{1}^{(1)}(\beta^{(1)}_1)=L_{1,2}^{(1)}\left (0\right)=b_{1}^{(1)}(0)=-4(0-.25)=1$$


 * $$L_{1,2}^{(1)}\left (.25\right)=b_{1}^{(1)}(.25)=-4(.25-.25)=0$$

Similarly for the second node, the basis function has a value of one at that second node and a value of zero at the first node.

A Matlab script was written to plot the LEBFs at each node for each element. The Matlab script is given followed by the plots. The first plot shows each of the LEBFs plotted across the entire domain. The second plot is more informative as the plot is divided into the elemental domains and lines corresponding to the values zero and one are present. The legend states the node number, $$n$$, and element number, $$e$$, of each LEBF.



As with the CBS, the coefficient at the essential boundary condition is equal to the essential boundary condition because the basis function is equal to one at the essential boundary condition:


 * $$g = \mathbf{d_{1}^{(1)}}\mathbf{b_{1}^{(1)}}(0)=\mathbf{d_{1}^{(1)}}*1=\mathbf{d_{1}^{(1)}}=4$$

The element matrices are formed as in problems 5.1 and 5.2 as follows:


 * $$\displaystyle \mathbf{F}^{e}_{F}=b_{i}^{(1)}(\alpha_{i}^{e} )h + \int_{\omega^{e}}^{.} {{b_{i}^{(1)}}fdx}$$



\displaystyle

h = {a_2}(\alpha )\frac = (2+3\alpha)\frac = 5*12/5=12

$$


 * $$\displaystyle \mathbf{F_{F}^{(1)}}=\mathbf{F}^{(1)}_{2}=b_{2}^{(1)}(.25)h + \int_{0}^{.25} {{b_{2}^{(1)}}fdx}=1*12+\int_{0}^{.25} 4x(5x)=12.104167$$

$$ \displaystyle

\therefore {F_F} = \left[ {\begin{array}{ccccccccccccccc} {12.104167} \end{array}} \right]

$$


 * $$\displaystyle \mathbf{K}^{e}_{ij}=\int_{\omega^{e}}^{.} b_{i}^{e'}(x^e)a_2(x^e)b_{j}^{e'}(x^e)dx$$



\displaystyle

\mathbf{K^{(1)}_{FE}} = \left\{ {\mathbf{K^{(1)}_{1j}};j = 2,...,n} \right\}

$$



\displaystyle

\mathbf{K^{(1)}_{12}}= \int\limits_0^{.25} b_{1}^{(1)'}a_2b_{2}^{(1)'}=\int\limits_0^{.25}-4*(2+3x)*4=-32(.25)-24(.25)^2=-9.5

$$



\displaystyle \therefore \mathbf{K^{(1)}_{FF}} = \begin{bmatrix} -9.5 \end{bmatrix} $$



\displaystyle

\mathbf{K^{(1)}_{FF}} = \left\{ {\mathbf{K^{(1)}_{ij}};i,j = 2,...,n} \right\}

$$



\displaystyle

\mathbf{K^{(1)}_{22}}= \int\limits_0^{.25} b_{2}^{(1)'}a_2b_{2}^{(1)'}=\int\limits_0^{.25}4*(2+3x)*4=32(.25)+24(.25)^2=9.5

$$



\displaystyle \therefore \mathbf{K^{(1)}_{FF}} = \begin{bmatrix} 9.5 \end{bmatrix} $$


 * $$\mathbf{K^{(1)}} \mathbf{d^{(1)}} = \mathbf{F_{F}^{(1)}} - \mathbf{K_{FE}^{(1)}} g$$


 * $$\mathbf{d_{2}^{(1)}} = \mathbf{K}_{FF}^{(1)-1}( \mathbf{F_{F}^{(1)}} - \mathbf{K_{FE}^{(1)}} g)$$


 * $$\therefore \mathbf{d_{2}^{(1)}} = 5.5208$$

Now the process is repeated for each element. The matrices are assembled into global matrices and all of the coefficients can be determined. The approximate solution is formed for each element by taking the coefficients times the basis functions. The approximate solutions for each element make up the global approximate solution.


 * $$\displaystyle \mathbf{K}_{ij} = \sum_{e=1}^{nel}\mathbf{K}^{e}_{ij}$$


 * $$\displaystyle \mathbf{F}_{i} = \sum_{e=1}^{nel}\mathbf{F}^{e}_{i}$$


 * $$\displaystyle u^{he}(x)= \sum_{j=1}^{n}d_jb_j(x),\quad \quad u^{h}(x)=\sum_{j=1}^{nel}u^{he}(x) $$

A Matlab script was written to generate the plots. The results and script are given below:



Problem Description
Solve the general 1-D model for an elastic bar with data set 1 (G1DM1.0/D1) using the Linear Lagrangian element basis function (LLEBF) with uniform discretization(equidistant element nodes, nel=4, 6, 8).
 * 1.Explain how the basis functions satisfy the constraint breaking solution (CBS).
 * 2.Plot these basis functions.
 * 3.Use matlab to integrate.
 * 4.Make Plots of u^h vs u and u^h(.5)-u(.5) vs m

Explain how LLEBF are used as CBS
This problem is similar to problem 5.7 except switch the location of natural and essential boundary conditions.

The Lagrangian Element Basis Function (LEBF) is given by the following equation and generates a polynomial of degree $$\displaystyle n-1$$:


 * $$L_{i,n}\left (x\right )=\prod_{k=1}^{n}\frac{x-x_{k}}{x_{i}-x_{k}},k\neq i$$

For a linear function the number of nodes is 2 ($$\displaystyle n=2$$). For uniform discretization with $$\displaystyle nel=4$$ the domain $$\displaystyle \Omega$$ will be subdivided into 4 equal length elements each having two nodes. The LEBFs for the first element $$\displaystyle n=1$$ over the elemental domain $$\displaystyle \omega^{(1)}=[0,.25]$$ can be written as follows:


 * $$L_{1,2}^{(1)}\left (x\right)=b_{1}^{(1)}=\frac{(x-x_2)}{x_1-x_2}=\frac{(x-.25)}{0-.25}=-4(x-.25)$$


 * $$L_{2,2}^{(1)}\left (x\right)=b_{2}^{(1)}=\frac{(x-x_1)}{x_2-x_1}=\frac{(x-0)}{.25-0}=4(x)$$

Recall the constraint breaking solution (CBS) states the basis functions must satisfy $$\displaystyle \left\{ b_i(x), i=0,1,2,...\right\}$$ such that $$\displaystyle b_0(\beta) \neq 0$$ and $$\displaystyle b_i(\beta) = 0$$ for $$\displaystyle i=1,2,...,n$$.

It can be seen that these basis functions for the first element satisfy the CBS at each node. For the first node, the basis function has a value of one at that node and a value of zero at every other node (there is only one other node for the linear LEBF). In other words the basis function is equal to the kronicker delta at each node.


 * $$b_{1}^{(1)}(\beta^{(1)}_1)=L_{1,2}^{(1)}\left (0\right)=b_{1}^{(1)}(0)=-4(0-.25)=1$$


 * $$L_{1,2}^{(1)}\left (.25\right)=b_{1}^{(1)}(.25)=-4(.25-.25)=0$$

Similarly for the second node, the basis function has a value of one at that second node and a value of zero at the first node.

Therefore, LLEBF satifies CBS