User:Eml5526.s11.team3/Homework 5

 Homework 5

=Problem 5.1 Using appropriate basis functions (Polynomial, Fourier and Exponential basis functions) to obtain the trial solution from Weak Form=

Problem Statement
$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \qquad \forall x \in \left[ 0,1 \right] $$,

Essential Boundary Function  $$ u \left(0\right)=4 $$, and Natural Boundary Function  $$ n(1)a_2(1)\frac{\partial u(x=1)}{\partial x}=12$$ $$\implies (1)(2+3)\frac{\partial u(x=1)}{\partial x}=12$$

$$\implies \frac{\partial u(x=1)}{\partial x}=\frac{12}{5}

$$

$$ \Gamma_g = \left[0\right] \qquad \Gamma_h = \left[1\right] $$

Find
Solve the above problem using Weight Function with appropriate basis function (Polynomial basis function, Fourier basis function and Exponential basis function) until convergence of $$u^h \left(0.5\right)$$ to $$10^{-6}$$.

Plot $$ u, \quad u^h $$, convergence (error vs. n)

Solution
The strong form of $$\text{G1DM1}.0/\text{D1b}$$ is:
 * {| style="width:100%" border="0"

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

\begin{align} & \frac{d}{dx}\left((2+3x)\frac{du}{dx}\right) +5x=0,_ – ^ – \forall x\in (0,1) \\ & u(x=0)=4,_ – ^ – \frac{du}{dx}(x=1)=\frac{12}{5} \\ \end{align}

$$ We can derive the weak form as:
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & \int\limits_{0}^{1}{\frac{\partial w}{\partial x}(2+3x)\frac{\partial u}{\partial x}}dx-\int\limits_{0}^{1}{5xwdx}-(\frac{24}{5}+\frac{36}{5}x)w\left| _{x=1} \right.=0,_ – ^ – \forall w(x=0)=0\\ & u(x=0)=4 \\ \end{align}

$$ The discretized weak form is
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & \underbrace{ [\int\limits_{0}^{1}{\frac{\partial {{b}_{i}}(x)}{\partial x}(2+3x)\frac{\partial {{b}_{j}}(x)}{\partial x}}dx ]}_{{d}_{i}}=\underbrace{[\int\limits_{0}^{1}{5x{{b}_{i}}(x)dx}+(\frac{24}{5}+\frac{36}{5}x){{b}_{i}}(x)\left| _{x=1} \right.]}_ \\ & {{d}_{0}}=4 \\ & i=1,2,\cdots \\ & j=1,2,\cdots \\ \end{align}

$$ The appropriate Polynomial basis function is
 * }
 * {| style="width:100%" border="0"

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

{{f}_{p}}=\{1,{{x}^{j}},j=1,2,\cdots \}

$$ The appropriate Fourier basis function is
 * }
 * {| style="width:100%" border="0"

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

{{f}_{F}}=\{1,\cos (jx)-1,\sin (kx),j=1,2,\cdots ,k=1,2,\cdots \}

$$ The appropriate Exponential basis function is
 * }
 * {| style="width:100%" border="0"

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

{{f}_{e}}=\{1,{{e}^{jx}}-1,j=1,2,\cdots \}

$$ The exact solution of governing equation is
 * }
 * {| style="width:100%" border="0"

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

u(x)=\frac{241}{54}\ln (\frac{2+3x}{2})+\frac{1}{36}(20x-15{{x}^{2}}+144)

$$
 * }

Using Polynomial basis function
The trial solution is
 * {| style="width:100%" border="0"

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

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

$$ Since the trial solution must satisfy the essential boundary condition,$$u \left(x=0\right)=4$$, substituting x=0 into the trial solution, we have   $${{d}_{0}}=4$$.
 * }

In order to solve for the rest coefficient $${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the stiffness matrix and force matrix as follow:
 * {| style="width:100%" border="0"

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

{{K}_{F\times F}}=[{{k}_{ij}}]=[\int\limits_{0}^{1}{i{{x}^{i-1}}(2+3x)j{{x}^{j-1}}}dx]

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

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

{{F}_{F\times 1}}=[{{f}_{i\times 1}}]=[\int\limits_{0}^{1}{5x{{x}^{i}}}+(\frac{24}{5}+\frac{36}{5}*1){{1}^{i}}]

$$ Using Matlab to construct above two 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" |

{{K}_{FF}}=\left( \begin{matrix} 7/2 & 4 & 17/4 & 22/5 & 9/2 & 32/7 & 37/8 & 14/3\\ 4 & 17/3 & 33/5 & 36/5 & 160/21 & 111/14 & 49/6 & 376/45\\ 17/4 & 33/5 & 81/10 & 64/7 & 555/56 & 21/2 & 329/30 & 624/55\\ 22/5 & 36/5 & 64/7 & 74/7 & 35/3 & 188/15 & 728/55 & 152/11\\ 9/2 & 160/21 & 555/56 & 35/3 & 235/18 & 156/11 & 665/44 & 620/39\\ 32/7 & 111/14 & 21/2 & 188/15 & 156/11 & 171/11 & 217/13 & 1608/91\\ 37/8 & 49/6 & 329/30 & 728/55 & 665/44 & 217/13 & 469/26 & 96/5\\ 14/3 & 376/45 & 624/55 & 152/11 & 620/39 & 1608/91 & 96/5 & 308/15\\ \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}_{F\times 1}}={\left( \begin{matrix} 41/3 \\ 53/4 \\ 13 \\ 77/6 \\ 89/7 \\ 101/8 \\ 113/9 \\ 25/2\\ \end{matrix} \right)}

$$ By the equation that $$ Kd=F  $$, we can solve the rest of $$d$$ 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" |

{{d}_{F\times 1}}={\left( \begin{matrix} -0.23479307 \\ 1.2522298 \\ -3.0192656 \\ 4.5191594 \\ -5.0302374 \\ 4.9208409 \\ -5.4294546 \\ 7.2497811\\ \end{matrix} \right)}

$$ Hence, the trial solution of the governing equation is
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \begin{align} & U_{F}^{h}(x)= - 0.23479307\;x^8 + 1.2522298\;x^7 - 3.0192656\;x^6 + 4.5191594\;x^5\\ & - 5.0302374\;x^4 + 4.9208409\;x^3 - 5.4294546\;x^2 + 7.2497811\;x + 4\\ \end{align}
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$
 * }



Using Fourier basis function
The trial solution is
 * {| style="width:100%" border="0"

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

\begin{align} & u_{F}^{h}(x)={{d}_{0}}[1]+{{d}_{1}}[\cos [1*(x)]-1]+{{d}_{2}}\sin [1*(x)]+\cdots \\ & \begin{matrix} {} & {} \\ \end{matrix}+{{d}_{j-1}}[\cos [\frac{j}{2}(x)]-1]+{{d}_{j}}\sin [\frac{j}{2}(x)],j=0,1,2,\cdots  \\ \end{align}

$$ Since the trial solution must satisfy the essential boundary condition, $$u \left(x=0\right)=4$$, substituting $$x=0$$  into the trial solution, we have $${{d}_{0}}=4$$. In order to solve for the rest coefficient$${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the stiffness matrix and force matrix as follow:
 * }
 * {| style="width:100%" border="0"

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

{{K}_{F\times F}}=[{{k}_{ij}}]= \begin{cases} \int\limits_{0}^{1}{i\cos i(x)(2+3x)j\cos j(x)}dx,_ – ^ – \quad when\quad i, j\quad both \quad are \quad even \quad number \\ \int\limits_{0}^{1}{i\sin i(x)(2+3x)j\sin j(x)}dx,_ – ^ – \quad when\quad i, j\quad both \quad are \quad odd \quad number \\ \int\limits_{0}^{1}{i\cos i(x)(2+3x)j\sin j(x)}dx, \quad others \\ \end{cases}

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

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

{{F}_{F\times 1}}=[{{f}_{i\times 1}}]=\begin{cases} \int\limits_{0}^{1}{5x\sin i(x)}dx-(\frac{24}{5}+\frac{36}{5}1)\sin i(1), \quad \text{when i is even number }  \\ \int\limits_{0}^{1}{5x[\cos i(x)-1]}dx-(\frac{24}{5}+\frac{36}{5}1)[\cos i(1)-1], \quad \text{when i is odd number } \\ \end{cases}

$$ Using Matlab to construct above two 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" |

{{K}_{FF}}=\left( \begin{matrix} 1.1444 & -1.3612  &  3.2564 &  -0.5408  &  3.7008  &  2.3208  &  1.2335 &   4.7249\\ -1.3612 &   2.3556 &  -4.1866  &  2.4001 &  -5.8462 &  -0.0674 &  -4.7299 &  -2.9460\\ 3.2564 &  -4.1866  &  9.5121 &  -2.3503 &  11.6194  &  5.4649  &  5.7892  & 12.8497\\ -0.5408 &   2.4001  & -2.3503  &  4.4879 &  -5.4725 &   5.3501 &  -8.9290 &   3.8999\\ 3.7008 &  -5.8462 &  11.6194 &  -5.4725  & 16.8127 &   2.2054 &  14.2443 &  12.2128\\ 2.3208 &  -0.0674  &  5.4649  &  5.3501  &  2.2054 &  14.6873 &  -9.6620 &  19.6948\\ 1.2335 &  -4.7299  &  5.7892  & -8.9290 &  14.2443  & -9.6620 &  23.4828 &  -3.0985\\ 4.7249 &  -2.9460 &  12.8497  &  3.8999 &  12.2128 &  19.6948 &  -3.0985 &  32.5172\\ \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}_{F\times 1}}={\left( \begin{matrix}  -6.1075 \\ 11.6035 \\ -18.9907 \\  13.0886 \\ -27.2503 \\  3.4218 \\ -23.8065 \\ -8.5011\\ \end{matrix} \right)}

$$ By the equation that$$Kd=F$$, we can solve the rest of $$d$$ 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" |

{{d}_{F\times 1}}={\left( \begin{matrix} 0.90619376 \\ -4.1185023 \\ -3.4991531 \\ -30.621018 \\ -7.205183 \\ -0.42479515 \\ -47.408642 \\ 48.575238 \\ \end{matrix} \right)}

$$ Hence, the trial solution of the governing equation is
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & U_{F}^{h}(x)=\text{0}.\text{9}0\text{619376cos}\left( \text{4x} \right)-\text{4}.\text{1185023cos}\left( \text{3x} \right)-\text{3}.\text{4991531cos}\left( \text{2x}\right) \\ & -\text{30}.\text{621018sin}\left( \text{2x} \right)-\text{7}.\text{205183sin}\left( \text{3x}\right)-\text{0}.\text{42479515sin}\left( \text{4x}\right) \\ & -\text{47}.\text{408642cos}\left( \text{x}\right)+\text{48}.\text{575238sin}\left( \text{x}\right)+\text{36}.\text{69718}\\ \end{align}

$$
 * }



Using Exponential basis function
The trial solution is
 * {| style="width:100%" border="0"

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

u_{e}^{h}(x)={{d}_{0}}(1)+{{d}_{1}}({{e}^{1*(x)}}-1)+{{d}_{2}}({{e}^{2*(x)}}-1)+\cdots +{{d}_{j}}({{e}^{j*(x)}}-1),_ – ^ – j=0,1,2,\cdots

$$ Since the trial solution must satisfy the essential boundary condition, $$u\left(x=0\right)=4$$, substituting      $$x=0$$   into the trial solution, we have $${{d}_{0}}=4$$. In order to solve for the rest coefficient$${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the stiffness matrix and force matrix as follow:
 * }
 * {| style="width:100%" border="0"

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

{{K}_{F\times F}}=[{{k}_{ij}}]=\left[\int\limits_{0}^{1}{i{{e}^{i(x)}}(2+3x)j{{e}^{j(x)}}}dx\right]

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

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

{{F}_{F\times 1}}=[{{f}_{i\times 1}}]=\left[\int\limits_{0}^{1}{5x[{{e}^{i(x)}}-1]dx}-(\frac{24}{5}+\frac{36}{5}*1)[{{e}^{i(1)}}-1]\right]

$$
 * }

Using Matlab, if we make the convergence criterion a little larger, that’s to $${{10}^{-5}}$$, we can get following results.


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

$$\displaystyle {{K}_{FF}}=1.0e+009 *\left( \begin{matrix} 0.0000 &  0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0001 \\ 0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0001  &  0.0002  &  0.0005 \\ 0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0001  &  0.0002  &  0.0006  &  0.0017 \\ 0.0000  &  0.0000  &  0.0000  &  0.0000  &  0.0001  &  0.0002  &  0.0007  &  0.0021  &  0.0058 \\ 0.0000  &  0.0000  &  0.0000  &  0.0001  &  0.0003  &  0.0008  &  0.0023  &  0.0065  &  0.0185 \\ 0.0000  &  0.0000  &  0.0001  &  0.0002  &  0.0008  &  0.0023  &  0.0068  &  0.0197  &  0.0565 \\ 0.0000  &  0.0001  &  0.0002  &  0.0007  &  0.0023  &  0.0068  &  0.0201  &  0.0586  &  0.1684 \\ 0.0000  &  0.0002  &  0.0006  &  0.0021  &  0.0065  &  0.0197  &  0.0586  &  0.1711  &  0.4935 \\ 0.0001  &  0.0005  &  0.0017  &  0.0058  &  0.0185  &  0.0565  &  0.1684  &  0.4935  &  1.4281 \\
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

\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}_{F\times 1}}=1.0e+005 *

$$ By the equation that$$Kd=F$$, we can solve the rest of $$d$$ as Hence, the trial solution of the governing equation is
 * }
 * {| style="width:100%" border="0"

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

\begin{align} &u_{e}^{h}(x) = 507.7233\;e^{3x} - 455.35367\;e^{2x} - 371.33188\;e^{4x}+ 182.56105\;e^{5x} \\ & - 59.974767\;e^{6x} + 12.652506\;e^{7x} - 1.5521958\;e^{8x}+ 0.084261502\;e^{9x} \\ & + 250.24399\;e^{x} - 61.0526\\ \end{align}

$$
 * }



Comment
When we use the exponential basis function($$e^{(i(x))}-1$$) to obtain the trial solution from discretized weak form, the error at point x=0.5 never converges to $$10^{-6}$$. But if we change the accuracy to $$10^{-5}$$, it will converges to this far very fast.Among the two criterion for convergence, the continuity is satisfied, while the completeness is not. It's possible because that the exponential basis function doesn't contain a linear term or it can't used to approximate a linear function. So it's incomplete. There is also another problem rises from this assumption. Since there is certain equivalent relationship between exponential basis function and Fourier basis function, the same may also happen for Fourier basis function. But the Fourier basis function does converge to the accuracy of $$10^{-6}$$ here. Since it's much expensive to run the code for Fourier basis function, especially to a much accurate convergence, we don't have the chance to testify this inference.

=Problem 5.2 Using appropriate basis functions (Polynomial, Fourier and Exponential basis functions) to obtain the trial solution from Weak Form=

Problem Statement
$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \qquad \forall x \in \left[ 0,1 \right] $$, Essential Boundary Function  $$ u \left(1\right)=4 $$, and

Natural Boundary Function $$ -\frac{\partial u(x=0)}{\partial x}=6 $$

$$ \Gamma_g = \left[1\right] \qquad \Gamma_h = \left[0\right] $$

Find
Solve the above problem using Weight Function with appropriate basis function (Polynomial basis function, Fourier basis function and Exponential basis function) until convergence of $$u^h \left(0.5\right)$$ to $$10^{-6}$$.

Plot $$ u, \quad u^h $$, convergence (error vs. n)

Solution
The strong form of $$\text{G1DM1}.0/\text{D1}$$ is:
 * {| style="width:100%" border="0"

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

\begin{align} & \frac{d}{dx}[(2+3x)\frac{du}{dx}]+5x=0,_ – ^ – \forall x\in (0,1) \\ & u(x=1)=4,_ – ^ – -\frac{du}{dx}(x=0)=6 \\ \end{align}

$$ We can derive the weak form as:
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & \int\limits_{0}^{1}{\frac{\partial w}{\partial x}(2+3x)\frac{\partial u}{\partial x}}dx-\int\limits_{0}^{1}{5xwdx}-(12+18x)w\left| _{x=0} \right.=0,_ – ^ – \forall w(x=1)=0 \\ & u(x=1)=4 \\ \end{align}

$$ The discretized weak form is
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & \underbrace{[\int\limits_{0}^{1}{\frac{\partial {{b}_{i}}(x)}{\partial x}(2+3x)\frac{\partial {{b}_{j}}(x)}{\partial x}}dx]}_{{d}_{i}}=\underbrace{[\int\limits_{0}^{1}{5x{{b}_{i}}(x)dx}+(12+18x){{b}_{i}}(x)\left| _{x=0} \right.]}_ \\ & {{d}_{0}}=4 \\ & i=1,2,\cdots \\ & j=1,2,\cdots \\ \end{align}

$$ The appropriate Polynomial basis function is
 * }
 * {| style="width:100%" border="0"

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

{{f}_{p}}=\{1,{{(x-1)}^{j}},j=1,2,\cdots \}

$$ The appropriate Fourier basis function is
 * }
 * {| style="width:100%" border="0"

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

{{f}_{F}}=\{1,\cos j(x-1)-1,\sin k(x-1),j=1,2,\cdots ,k=1,2,\cdots \}

$$ The appropriate Exponential basis function is
 * }
 * {| style="width:100%" border="0"

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

{{f}_{e}}=\{1,{{e}^{j(x-1)}}-1,j=1,2,\cdots \}

$$ The exact solution of governing equation is
 * }
 * {| style="width:100%" border="0"

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

u(x)=-\frac{118}{27}\ln (\frac{2+3x}{5})+\frac{1}{36}(20x-15{{x}^{2}}+139)

$$
 * }

Using Polynomial basis function
The trial solution is
 * {| style="width:100%" border="0"

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

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

$$ Since the trial solution must satisfy the essential boundary condition,$$u(x=1)=4$$, substituting $$x=1$$into the trial solution, we have $${{d}_{0}}=4$$. In order to solve for the rest coefficient$${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the stiffness matrix and force matrix as follow:
 * }
 * {| style="width:100%" border="0"

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

{{K}_{F\times F}}=[{{k}_{ij}}]=[\int\limits_{0}^{1}{i{{(x-1)}^{i-1}}(2+3x)j{{(x-1)}^{j-1}}}dx]

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

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

{{F}_{F\times 1}}=[{{f}_{i\times 1}}]=[\int\limits_{0}^{1}{5x{{(x-1)}^{i}}}+(12+18*0){{(0-1)}^{i}}]

$$ Using Matlab to construct above two 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" |

{{K}_{FF}}=\left( \begin{matrix}  3.5 & -3.0 & 2.75 & -2.6 & 2.5 & -2.4285714 & 2.375 & -2.3333333  \\   -3.0 & 3.6666667 & -3.9 & 4.0 & -4.047619 & 4.0714286 & -4.0833333 & 4.0888889  \\   2.75 & -3.9 & 4.5 & -4.8571429 & 5.0892857 & -5.25 & 5.3666667 & -5.4545455  \\   -2.6 & 4.0 & -4.8571429 & 5.4285714 & -5.8333333 & 6.1333333 & -6.3636364 & 6.5454545  \\   2.5 & -4.047619 & 5.0892857 & -5.8333333 & 6.3888889 & -6.8181818 & 7.1590909 & -7.4358974  \\   -2.4285714 & 4.0714286 & -5.25 & 6.1333333 & -6.8181818 & 7.3636364 & -7.8076923 & 8.1758242  \\   2.375 & -4.0833333 & 5.3666667 & -6.3636364 & 7.1590909 & -7.8076923 & 8.3461538 & -8.8  \\   -2.3333333 & 4.0888889 & -5.4545455 & 6.5454545 & -7.4358974 & 8.1758242 & -8.8 & 9.3333333  \\ \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}_{F\times 1}}={{\left( \begin{matrix}  -12.833333 & 12.416667 & -12.15 & 12.166667 & -12.119048 & 12.089286 & -12.069444 & 12.055556  \\ \end{matrix} \right)}^{T}}

$$ By the equation that$$Kd=F$$, we can solve the rest of $$d$$ 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" |

{{d}_{F\times 1}}={{\left( \begin{matrix}  -2.8999142 & 0.37306134 & -0.27987546 & 0.32403056 & 0.43872514 & 0.81068777 & 0.61312521 & 0.22992195  \\ \end{matrix} \right)}^{T}}

$$ Hence, the trial solution of the governing equation is
 * }
 * {| style="width:100%" border="0"

$$\displaystyle \begin{align} &u_{p}^{h}(x)=4-2.8999142(x-1)+0.37306134{{(x-1)}^{2}}-0.27987546{{(x-1)}^{3}}\\ & +0.32403056{{(x-1)}^{4}}+0.43872514{{(x-1)}^{5}}+0.81068777{{(x-1)}^{6}}\\ &+0.61312521{{(x-1)}^{7}}+0.22992195{{(x-1)}^{8}}\\ \end{align} $$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * }



Using Fourier basis function
The trial solution is
 * {| style="width:100%" border="0"

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

\begin{align} & u_{F}^{h}(x)={{d}_{0}}*1+{d}_{1}[\cos [1*(x-1)]-1]+{{d}_{2}}\sin [1*(x-1)]+\cdots \\ & \begin{matrix} {} & {} \\ \end{matrix}+{{d}_{j-1}}[\cos [\frac{j}{2}(x-1)]-1]+{{d}_{j}}\sin [\frac{j}{2}(x-1)],j=0,1,2,\cdots  \\ \end{align}

$$ Since the trial solution must satisfy the essential boundary condition,$$u(x=1)=4$$, substituting $$x=1$$into the trial solution, we have $${{d}_{0}}=4$$. In order to solve for the rest coefficient$${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the stiffness matrix and force matrix as follow:
 * }
 * {| style="width:100%" border="0"

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

{{K}_{F\times F}}=[{{k}_{ij}}]=[\left\{ \begin{matrix} \int\limits_{0}^{1}{i\cos i(x-1)(2+3x)j\cos j(x-1)}dx,_ – ^ – \quad when\quad i, j\quad both \quad are \quad even \quad number  \\ \int\limits_{0}^{1}{i\sin i(x-1)(2+3x)j\sin j(x-1)}dx,_ – ^ – \quad when\quad i, j\quad both \quad are \quad odd \quad number \\ \int\limits_{0}^{1}{i\cos i(x-1)(2+3x)j\sin j(x-1)}dx,\quad others \\ \end{matrix} \right.]

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

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

{{F}_{F\times 1}}=[{{f}_{i\times 1}}]=[\left\{ \begin{matrix} \int\limits_{0}^{1}{5x\sin i(x-1)}dx-(12+18*0)\sin i(0-1),\quad \text{when i is even number } \\ \int\limits_{0}^{1}{5x[\cos i(x-1)-1]}dx-(12+18*0)[\cos i(0-1)-1],\quad \text{when i is odd number } \\ \end{matrix} \right.]

$$ Using Matlab to construct above two 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" |

{{K}_{FF}}=\left( \begin{matrix}  0.76429622 & 1.1170869 & 2.3046242 & 0.88463923 & 3.0596259 & -0.77317723 & 2.11004 & -2.5560312 & 0.05369072  \\   1.1170869 & 2.7357038 & 3.6746041 & 3.8194459 & 5.9293441 & 2.8545967 & 6.5624554 & 0.91960333 & 5.6382087  \\   2.3046242 & 3.6746041 & 7.1366861 & 3.4374446 & 10.078933 & -1.1801292 & 8.2449076 & -6.7904503 & 2.7731659  \\   0.88463923 & 3.8194459 & 3.4374446 & 6.8633139 & 7.1897654 & 8.2932776 & 11.08294 & 7.5263528 & 13.691415  \\   3.0596259 & 5.9293441 & 10.078933 & 7.1897654 & 16.154187 & 2.4145154 & 17.155544 & -5.6179262 & 12.527795  \\   -0.77317723 & 2.8545967 & -1.1801292 & 8.2932776 & 2.4145254 & 15.345813 & 11.121863 & 19.588876 & 21.368795  \\   2.11004 & 6.5624554 & 8.2449076 & 11.08294 & 17.155544 & 11.121863 & 25.591721 & 4.9199907 & 28.996706  \\   -2.5560312 & 0.91960333 & -6.7904503 & 7.5263528 & -5.6179262 & 19.588876 & 4.9199907 & 30.408279 & 21.377477  \\ 0.05369072 & 5.6382087 & 2.7731659 & 13.691415 & 12.527795 & 21.368795 & 28.996706 & 21.377477 & 44.420401 \\ \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}_{F\times 1}}={{\left( \begin{matrix}  -5.7178839 & -10.890297 & -17.723578 & -12.274947 & -25.274359 & -3.2817068 & -21.82696 & 7.5951292 & -10.952786  \\ \end{matrix} \right)}^{T}}

$$ By the equation that$$Kd=F$$, we can solve the rest of $$d$$ 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" |

{{d}_{F\times 1}}={{\left( \begin{matrix}  7.4103829 & -40.592257 & 10.904443 & 35.696817 & -12.476341 & -14.112008 & 4.7711925 & 2.1586481 & -0.6326032  \\ \end{matrix} \right)}^{T}}

$$ Hence, the trial solution of the governing equation is
 * }
 * {| style="width:100%" border="0"

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

\begin{align} & U_{F}^{h}(x)=\text{1}0.\text{9}0\text{4443cos}\left( \text{2x}-\text{2} \right)-\text{4}0.\text{592257sin}\left( \text{x}-\text{1} \right)-\text{12}.\text{476341cos}\left( \text{3x}-\text{3} \right) \\ & +\text{4}.\text{7711925cos}\left( \text{4x}-\text{4} \right)-0.\text{6326}0\text{32cos}\left( \text{5x}-\text{5} \right)+\text{35}.\text{696817sin}\left( \text{2x}-\text{2} \right) \\ & -\text{14}.\text{112}00\text{8sin}\left( \text{3x}-\text{3} \right)+\text{2}.\text{1586481sin}\left( \text{4x}-\text{4} \right)+\text{7}.\text{41}0\text{3829cos}\left( \text{x}-\text{1} \right)-\text{5}.\text{977}0\text{749} \\ \end{align}

$$
 * }



Using Exponential basis function
The trial solution is
 * {| style="width:100%" border="0"

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

u_{e}^{h}(x)={{d}_{0}}*1+{{d}_{1}}({{e}^{1*(x-1)}}-1)+{{d}_{2}}({{e}^{2*(x-1)}}-1)+\cdots +{{d}_{j}}({{e}^{j*(x-1)}}-1),_ – ^ – j=0,1,2,\cdots

$$ Since the trial solution must satisfy the essential boundary condition,$$u(x=1)=4$$, substituting $$x=1$$into the trial solution, we have $${{d}_{0}}=4$$. In order to solve for the rest coefficient$${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the stiffness matrix and force matrix as follow:
 * }
 * {| style="width:100%" border="0"

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

{{K}_{F\times F}}=[{{k}_{ij}}]=[\int\limits_{0}^{1}{i{{e}^{i(x-1)}}(2+3x)j{{e}^{j(x-1)}}}dx]

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

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

{{F}_{F\times 1}}=[{{f}_{i\times 1}}]=[\int\limits_{0}^{1}{5x[{{e}^{i(x-1)}}-1]dx}-(12+18*0)[{{e}^{i(0-1)}}-1]]

$$ Use Matlab to construct above two matrices until satisfying the convergence criterion, the trial solution at point x=0.5 can’t converge to$${{10}^{-6}}$$. The plot of convergence using 26 exponential basis functions as follow:
 * }



If we make the convergence criterion a little larger, that’s to $${{10}^{-5}}$$, we can get following results.



By the equation that$$Kd=F$$, we can solve the rest of $$d$$ as Hence, the trial solution of the governing equation is
 * {| style="width:100%" border="0"

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

\begin{align} & u_{e}^{h}(x)=-\text{55}0.\text{21555}{{\text{e}}^{\left( \text{x}-\text{1} \right)}}+\text{2394}.\text{6264}{{\text{e}}^{\left( \text{2x}-\text{2} \right)}}-\text{5543}.\text{4261}{{\text{e}}^{\left( \text{3x}-\text{3} \right)}}+\text{572}0.\text{8768}{{\text{e}}^{\left( \text{4x}-\text{4} \right)}} \\ & +\text{2321}.\text{5781}{{\text{e}}^{\left( \text{5x }-\text{5} \right)}}-\text{11751}.\text{37}{{\text{e}}^{\left( \text{6x}-\text{6} \right)}}+\text{4594}.\text{3739}{{\text{e}}^{\left( \text{7x}-\text{7} \right)}}+\text{19997}.\text{55}{{\text{e}}^{\left( \text{8x}-\text{8} \right)}} \\ & -\text{36338}.\text{561}{{\text{e}}^{\left( \text{9x}-\text{9} \right)}}+\text{286}0\text{2}.\text{493}{{\text{e}}^{\left( \text{1}0\text{x}-\text{1}0 \right)}}-\text{11373}.\text{731}{{\text{e}}^{\left( \text{11x}-\text{11} \right)}}+\text{1866}.\text{4427}{{\text{e}}^{\left( \text{12x}-\text{12} \right)}}+\text{63}.\text{363186} \\ \end{align}

$$
 * }



Comment
When we use the exponential basis function($$e^{(i(x-1))}-1$$) to obtain the trial solution from discretized weak form, the error at point x=0.5 never converges to $$10^{-6}$$. But if we change the accuracy to $$10^{-5}$$, it will converges to this far very fast.Among the two criterion for convergence, the continuity is satisfied, while the completeness is not. It's possible because that the exponential basis function doesn't contain a linear term or it can't used to approximate a linear function. So it's incomplete. There is also another problem rises from this assumption. Since there is certain equivalent relationship between exponential basis function and Fourier basis function, the same may also happen for Fourier basis function. But the Fourier basis function does converge to the accuracy of $$10^{-6}$$ here. Since it's much expensive to run the code for Fourier basis function, especially to a much accurate convergence, we don't have the chance to testify this inference.

=Problem 5.3 Using Lagrange Interpolation Basis Functions to Obtain Trial and Exact Solution of Torsion Bar=

Given: Data set G1DM1.0/D1b defined on Lecture 26-2
Strong form for the Torsion Bar

where $${a_2}\left( x \right){\text{ =  2 + 3x}}$$  $$\qquad\text{and}\quad $$   $$f\left( {x,t} \right) = 5x $$  so we get Essential boundary conditions, Natural boundary conditions,

Find
1. Explain how Langrange Interpolation Basis Functions (LIBF) are used as Constraint Breaking Solutions

2. Plot all LIBF

3. Use matlab quad, WolframAlpha, ... to integrate

4. Plot $$ u_m^h  \ vs\   u$$ and $$ u_m^h\left(0.5 \right)-u\left(0.5 \right)\ vs\ m$$

1. Explaining how LIBF can be used as CBS
For a basis function, and for $$\displaystyle{{\Gamma }_{g}}=\{\beta \}$$, the corresponding basis functions satisfying

{ $$\displaystyle {{b}_{j}}(x), \quad j=1,2...n $$ }   such that   $$\displaystyle {{b}_{1}}(\beta)\ne 0$$   and     $$\displaystyle {{b}_{j}}(\beta)= 0, \quad  j = 2,3...n $$

Lets show that the statement above can be satisfied to explain how LIBF can be used as CBS

We know that general formula of Lagrange Interpolation Basis Function is

Lets denote m=4,

If $$\displaystyle x={{x}_{1}}, \qquad {b}_{1}=L_{1,4}(x_{1})= \delta_{1,1} =\frac{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})} =1 $$

If $$\displaystyle x={{x}_{2}},\qquad {b}_{2}=L_{1,4}(x_{2})= \delta_{1,2}= \frac{(x_{2}-x_{2})(x_{2}-x_{3})(x_{2}-x_{4})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})} = 0$$

If $$\displaystyle x={{x}_{3}},\qquad {b}_{3}= L_{1,4}(x_{3})= \delta_{1,3}= \frac{(x_{3}-x_{2})(x_{3}-x_{3})(x_{3}-x_{4})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})} = 0$$

If $$\displaystyle x={{x}_{4}}, \qquad {b}_{4}= L_{1,4}(x_{4})= \delta_{1,4}=\frac{(x_{4}-x_{2})(x_{4}-x_{3})(x_{4}-x_{4})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})} = 0$$

So for the Lagrange basis function $$\displaystyle {{b}_{j}}(x) $$ satisfy the Constraint Breaking Solution requirements since $$\displaystyle {b}_{1} =1, \quad and \quad {b}_{2}={b}_{3}={b}_{4} = 0 $$

And we can conclude that

We do same calculation for other m values such as m = 2,3,5,6,7,8,9 in Matlab program. As a result, we can choose Lagrange interpolation function as a basis because it will destroy the constraint.

2.Plotting all LIBF used for m = from 2 to 9
The LIBF used must satisfy the essential boundary conditions $$\displaystyle u\left( x=0 \right) = 4 \quad $$  according to Discrete Weak Form.

The essential boundary conditions to be satisfied are

$$\displaystyle {{b}_{1}}(0)=1 $$ and  $$\displaystyle {{b}_{2...n}}(0)=0 $$

Hence to satisfy the condition, the LIBF basis functions used

Using these basis functions, the LIBF are plotted for m = 2,3,4,5,6,7,8, and 9 until the convergence of the basis functions with the actual solution takes place where the error is less than $$\displaystyle  {{10}^{-6}} $$

The convergence becomes less than $$\displaystyle  {{10}^{-6}} $$ for  $$\displaystyle  m=9 $$. It is 1.0266 e-07

The LIBF are plotted for various values of $$\displaystyle m $$ until convergence.

















The plots of the basis functions were developed using Matlab. These LIBF plots were produced by the following Matlab code below.

i. Constracting of Exact Solution
We have general 1-Dimensional Model Data set 1.b from Eq 3.1 to Eq 3.5

Integrating the strong form Eq 3.2 with respect to x, we get,  -

From the natural boundary condition,

If we denote natural boundary condition, we obtain

If we substitute $$ \displaystyle C_{1} $$ into Eq 3.11

Integrating Eq 3.15 again with respect to x using integration by parts with constant $$ \displaystyle C_{2}$$ involved, we have

We already know this from Homework 2.4

If we take it out of integral, we have

Now using essential boundary condition with $$ \displaystyle u \left( {x = 0} \right) =4$$ ,we can update above equation as -

So we find

Substituting this constant values gives It is also equal to

ii. Constructing of Weak Form
As the weight function must vanish on the essential boundaries, we consider all smooth weight functions $$w(x)$$ such that $$w(0)= 0$$.

Multiplying the given differential equation and the natural boundary condition over the domain specified, by an arbitrary weight function

Now we integrate the Equation 3.24 by parts,

Integration by parts:

So

w(0) =0 therefore, the first term on the RHS of the above vanishes at x = 0. Substituting Equation 3.27 in Equation 3.24,

Substituting Equation 3.25 in Equation 3.28,

iii. Discretization of Weak Form
Approximated solution $$ u^{h} \left ( x \right ) $$ and $$ w^{h} \left ( x \right ) $$ :

where $$ \left ( c_{i} \right ) $$ and $$ \left ( d_{j} \right ) $$ are constants and $$ \left ( b_{i} \right ) $$ and $$ \left (b_{j}\right ) $$ are the LIBF. We divide equally to discretizate on our domain for m=2,3,4,5,6,7,8,9. We will use these LIBF as a shape function satisfying the CBS.

When we expand (Eq 3.30) for = 4, we get

where,

We get matrices K and F calculated in Matlab code below.

The convergence occurs at $$\displaystyle m=9 $$

Using Matlab to construct above two 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" |

{{K}_{FF}}=\left( \begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\   -48.861 & 167.423 & -276.316 &  345.621 & -349.434 & 260.149 & -142.253 & 55.498 & -11.827 \\    58.104 & -276.316 & 593.033 & -828.416 & 880.332 & -718.901 & 440.528 & -193.210 & 44.845 \\   -74.821 & 345.622 & -828.416 & 1349.981 & -1583.123 & 1382.496 & -919.426 & 432.555 & -104.867 \\    72.070 & -349.434 & 880.332 & -1583.123 & 2080.253 & -1954.709 & 1351.427 & -660.485 & 163.668\\    -49.752 & 260.150 & -718.901 & 1382.497 & -1954.709 & 2000.226 & -1453.569 & 712.071 & -178.013 \\    24.315 & -142.253 & 440.528 & -919.426 & 1351.427 & -1453.569 & 1152.701 & -595.975 &  142.252 \\    -8.185 &  55.498 & -193.211 &  432.555 & -660.485 &  712.071 & -595.975 & 376.137 & -118.406\\    1.520 & -11.827 & 44.845 &  -104.867 & 163.668 & -178.013 &  142.252 & -118.407 &  60.828\\ \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}_{F*1}}=\left( \begin{matrix}  4\\   0.130\\  -0.041\\    0.694\\  -0.400\\   1.157\\ -0.123\\   0.909\\   12.174\\ \end{matrix} \right)

$$
 * }

By the equation $$\displaystyle Kd =F $$, we 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" |

{{d}_{F*1}}=\left( \begin{matrix}  4.0000\\    4.8299\\    5.5341\\    6.1415\\    6.6712\\    7.1363\\    7.5463\\    7.9085\\    8.2283\\ \end{matrix} \right)

$$
 * }

4. Plotting Exact vs Trial Solution and the errors at x=0.5 for different m




=Problem 5.4 Using Lagrangian Interpolation basis functions to obtain the trial solution from Weak Form=

Problem statement
Please refer to lecture note 29-6 for more details.

Similar to HW 5.2, but using LIBF with uniform discretization (equidistant nodes) for n = 4, 6, 8

1.) Explain how LIBF are used as CBS

2.) Plot all LIBF used

3.) Use matlab quad, WA, ... to integrate the basis function

4.) Plot $$\displaystyle u_m^h $$ vs $$\displaystyle u $$  and    $$\displaystyle \quad $$ $$ u_m^h(0.5) - u(0.5) $$ vs.$$\displaystyle n $$

Using LIBF as CBS
i> For a basis function, and for $$\displaystyle{{\Gamma }_{g}}=\{\beta \}$$, the corresponding basis functions satisfying

{ $$\displaystyle {{b}_{j}}(x), \quad j=1,2...n $$ }   such that   $$\displaystyle {{b}_{1}}(\beta)\ne 0$$   and     $$\displaystyle {{b}_{j}}(\beta)= 0, \quad  j = 2,3...n $$

Let the LIBF be defined by a function


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

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

L_{i,n}(x) = \prod_{k = 1}^{m}\frac{x-x_{k}}{x_{i}-x_{k}}\;\; k\neq i $$ ; where $$\displaystyle n $$ is number of nodes.


 * }

Now for the given function, at $$\displaystyle x={{x}_{0}}, L_{0,n}(x)= \frac{({{x}_{0}}-{{x}_{1}})({{x}_{0}}-{{x}_{2}})({{x}_{0}}-{{x}_{3}})...({{x}_{0}}-{{x}_{n}})}{({{x}_{0}}-{{x}_{1}})({{x}_{0}}-{{x}_{2}})({{x}_{0}}-{{x}_{3}})...({{x}_{0}}-{{x}_{n}})} =1 $$

Now for the given function, at $$\displaystyle x={{x}_{1}}, L_{1,n}(x)= \frac{({{x}_{1}}-{{x}_{1}})({{x}_{1}}-{{x}_{2}})({{x}_{1}}-{{x}_{3}})...({{x}_{1}}-{{x}_{n}})}{({{x}_{0}}-{{x}_{1}})({{x}_{0}}-{{x}_{2}})({{x}_{0}}-{{x}_{3}})...({{x}_{0}}-{{x}_{n}})}= 0$$

Now for the given function, at $$\displaystyle x={{x}_{2}}, L_{2,n}(x)= \frac{({{x}_{2}}-{{x}_{1}})({{x}_{2}}-{{x}_{2}})({{x}_{2}}-{{x}_{3}})...({{x}_{2}}-{{x}_{n}})}{({{x}_{0}}-{{x}_{1}})({{x}_{0}}-{{x}_{2}})({{x}_{0}}-{{x}_{3}})...({{x}_{0}}-{{x}_{n}})}= 0$$

Now for the given function, at $$\displaystyle x={{x}_{n}},L_{n,n}(x)= \frac{({{x}_{n}}-{{x}_{1}})({{x}_{n}}-{{x}_{2}})({{x}_{n}}-{{x}_{3}})...({{x}_{n}}-{{x}_{n}})}{({{x}_{0}}-{{x}_{1}})({{x}_{0}}-{{x}_{2}})({{x}_{0}}-{{x}_{3}})...({{x}_{0}}-{{x}_{n}})}= 0$$

At $$\displaystyle x = x_{0}$$, the value of the basis function is 1 at node 1 and zero for all other nodes, which is required condition for any function to be the CBS

So for the basis function, $$\displaystyle {b}_{1} =1 , {b}_{2},{b}_{3},{b}_{n} = 0 $$

Hence the LIBF can be used as a CBS as it satisfies the conditions for the CBS.

Plotting the LIBF used
ii> For the LIBF to satisfy the Constraint Breaking Solution(CBS), it should satisfy the conditions of the CBS.

The CBS conditions to be satisfied are:-


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

$$\displaystyle {{b}_{1}}(1)=1 $$ and  $$\displaystyle {{b}_{2...n}}(1)=0 $$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * }

Hence to satisfy the condition, the LIBF basis functions used are:-


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

$$\displaystyle L_{i,n}(x) = \prod_{k = 1}^{m}\frac{1-x-x_{k}}{x_{i}-x_{k}}\;\; k\neq i $$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * }

The exact solution of governing equation is
 * {| style="width:100%" border="0"

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

u(x)=-\frac{118}{27}\ln (\frac{2+3x}{5})+\frac{1}{36}(20x-15{{x}^{2}}+139)

$$
 * }

Using these basis functions, the LIBF are plotted until the convergence of the basis functions with the actual solution takes place where the error is less than $$\displaystyle {{10}^{-6}} $$

The difference between the LIBF and the exact solution becomes less than $$\displaystyle  {{10}^{-6}} $$ for  $$\displaystyle  n=9 $$

The LIBF are plotted for various values of $$\displaystyle n $$ until convergence.

















Using Matlab to integrate the basis functions
iii> Using the above expressions, the plots of the basis functions were developed using Matlab

Plotting the LIBF and Exact Solution and the errors at x=0.5 for different values of n
iv> The convergence occurs at $$\displaystyle n=9 $$

Using Matlab to construct above two 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" |

{{K}_{FF}}=\left( \begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\   -118.406 & 376.137 & -595.975 &  712.071 & -660.485 & 432.555 & -193.210 & 55.498 & -8.185 \\   142.2520 & -595.975 & 1152.701 & -1453.567 & 1351.427 & -919.4265 & 440.528 & -142.253 & 24.315 \\   -178.012 & 712.071 & -1453.569 & 2000.226 & -1954.708 & 1382.496 & -718.900 & 260.149 & -49.751 \\   163.668 & -660.485 & 1351.427 & -1954.709 & 2080.253 & -1583.127 & 880.332 & -349.433 & 72.069\\   -104.8667 & 432.555 & -919.427 & 1382.497 & -1583.123 & 1349.980 & -828.416 & 345.622 & -74.822 \\   44.845 & -193.210 & 440.529 & -718.900 & 880.332 & -828.415 & 593.033 & -276.316 &  58.104 \\    -11.827 &  55.498 & -142.253 & 260.1492 & -349.433 &  345.622 & -276.3164 & 167.423 & -48.861\\   1.519 & -8.185 & 24.315 & -49.752 & 72.069 & -74.822 &  58.104 & -48.861 &  25.609\\ \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}_{F*1}}=\left( \begin{matrix}  4\\   0.908\\  -0.122\\    1.157\\  -0.400\\   0.694\\ -0.0409\\   0.129\\   12\\ \end{matrix} \right)

$$
 * }

By the equation $$\displaystyle Kd =F $$, we 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" |

{{d}_{F*1}}=\left( \begin{matrix}  4.0000\\    4.3689\\    4.7537\\    5.1595\\    5.5935\\    6.0649\\    6.5867\\    7.1775\\    7.8656\\ \end{matrix} \right)

$$
 * }





=Problem 5.5 Continue work with Calculix=

Problem statement
Please refer to lecture note [[media:fe1.s11.mtg29.djvu|29-7]] for more details.

Objective
1.) For the disk problem, extract:

node info: node numbers are coordinates

element info: element numbers and element nodes

2.) Generate 3 meshes of same disk with triangular elements (increase no. of elements)

3.) Install ccx, run examples, write report for "dummies" (explain commands, screenshots...)

1.)Extract node info and element info
Using the methods for creating a disc discussed in HW 4.7, we can begin by extracting node and element numbers from initial mesh. In order to get the node numbers and coordinates and element numbers and nodes, we need to create a mesh file. The command to generate a mesh file is


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

$$\displaystyle Send abq $$ Which is attached behind the mesh generating command. The node number and coordinates for this disc is as follow:
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * }

To do this graphically, we can simply type a few simple commands at the end of file we use to create out disc.

For showing element numbers on the plot use:
 * {| style="width:100%" border="0"

$$\displaystyle PLOT ea all $$ For showing node number on the plot use:
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * }
 * {| style="width:100%" border="0"

$$\displaystyle PLOT na all $$ Selecting Tools>Preprocess with display the window with the meshed geometry showing these numbers.
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * }

Left clicking on the space outside the display area will open a menu where you can select Viewing>Toggle Element Edges to show the elements.



Element numbers are show in the above picture, while nodal numbers are show below.



2.)Generating 3 meshes of same disk with triangular element type
Changing the type of elements and mesh density is not difficult once we have our initial geometry figured out. To change the type of elements to triangular elements, we simply must edit the line prefaced by “ELTY”. Previously we used quad meshes, which are created with by writing “ELTY QU8 All”. In order to change this to triangles, we can simple use “ELTY TR3 ALL”. Keep in mind that using a different element type might affect how the geometry is meshed and could cause bad elements to be created. Elements with volumes less than or equal to zero are considered bad and should be fixed before proceeding any further. Elements that are disproportionate to their neighboring elements in size and length may also cause inaccuracies in output data. To increase the mesh density, simply increase the number of divisions in the lines of the geometry. Be sure to take note of where the lines are, and how they interact with their neighboring lines. If you choose divisions that are wildly different, you could create an uneven mesh that could produce bad elements.

See the below input and output display for variations in mesh density.



3).Run CCX and Post-processing
To run ccx, you must have already installed Calculix. Save your .fdb file as an Input Deck File (.inp) to access the solver functions.In addition to the geometry and mesh, you must prescribe material properties, boundary conditions, and loads.

ELEMENT SETS:

This option is used to assign elements to an element set.

The parameter ELSET containing the name of the set is required (maximum 80 characters), whereas the parameter GENERATE (without value) is optional. If present, element ranges can be expressed by their initial value, their final value, and an increment. If a set with the same name already exists, it is reopened and complemented. The name of a set is case insensitive. Internally, it is modified into upper case and a 'E' is appended to denote it as element set.

First line:

•	*ELSET

•	Enter any needed parameters and their values.

Following line if the GENERATE parameter is omitted:

•	List of elements and/or sets of elements previously defined to be assigned to this element set (maximum 16 entries per line).

Repeat this line if needed.

Following line if the GENERATE parameter is included:

•	First element in set.

•	Last element in set.

•	Increment in element numbers between elements in the set. Default is 1.

Repeat this line if needed.

Example:


 * ELSET,ELSET=E1,GENERATE

20,25


 * ELSET,ELSET=E2

E1,50,51

assigns the elements with numbers 20, 21, 22, 23, 24 and 25 to element set E1 and the elements with numbers 20, 21, 22, 23, 24, 25 (= set E1), 50 and 51 to element set E2.

MATERIAL PROPERTIES: First line: •	*MATERIAL •	Enter the NAME parameter and its value. Example:


 * MATERIAL,NAME=EL

starts a material block with name EL This option is used to assign material properties to 3D, plane stress, plane strain and axisymmetric element sets. The parameters ELSET and MATERIAL are required, the parameter ORIENTATION is optional. The parameter ELSET defines the element set to which the material specified by the parameter MATERIAL applies. The parameter ORIENTATION allows to assign local axes to the element set. If activated, the material properties are applied to the local axis. This is only relevant for non isotropic material behavior. For plane stress and plane strain elements the thickness can be specified on the second line. Default is 1.

ELASTIC PROPERTIES: This option is used to define the elastic properties of a material. There is one optional parameter TYPE. Default is TYPE=ISO First line: •	*ELASTIC •	Enter the TYPE parameter and its value, if needed Following line for TYPE=ISO: •	Young's modulus. •	Poisson's ratio. •	Temperature. Example:

28000000, 0.3
 * ELASTIC

APPLYING MATERIAL PROPERTIES TO ELEMENTS: First line: •	*SOLID SECTION •	Enter any needed parameters. Second line (only relevant for plane stress, plane strain and axisymmetric elements; can be omitted for 3D elements): •	thickness (plane stress and plane strain elements) Example:


 * SOLID SECTION,MATERIAL=EL,ELSET=Eall,ORIENTATION=OR1

assigns material EL with orientation OR1 to all elements in (element) set Eall.

STEP: This card describes the start of a new STEP, and precedes the following.

STATIC: This procedure is used to perform a static analysis. The load consists of the sum of the load of the last *STATIC step and the load specified in the present step with replacement of redefined loads.

HOMOGENOUS BOUNDARY CONDITION: Homogeneous conditions should be placed before the first *STEP keyword card. First line: •	*BOUNDARY •	Enter any needed parameters and their value. Following line: •	Node number or node set label •	First degree of freedom constrained •	Last degree of freedom constrained. This field may be left blank if only one degree of freedom is constrained. Repeat this line if needed. Example:

73,1,3
 * BOUNDARY

fixes the degrees of freedom one through three (global if no transformation was defined for node 73, else local) of node 73.

CONCENTRATED LOAD: This option allows concentrated forces to be applied to any node in the model which is not fixed by a single or multiple point constraint. First line: •	*CLOAD •	Enter any needed parameters and their value. Following line: •	Node number or node set label. •	Degree of freedom. •	Magnitude of the load Repeat this line if needed. Example:

1000,3,10.3
 * CLOAD,OP=NEW,AMPLITUDE=A1,TIME DELAY=20.

removes all previous point load forces and applies a force with magnitude 10.3 and amplitude A1 (shifted in positive time direction by 20 time units) for degree of freedom three (global if no transformation was defined for node 1000, else local) of node 1000.

DISTRIBUTED LOAD: This option allows the specification of distributed loads. These include constant pressure loading on element faces and mass loading (load per unit mass) either by gravity forces or by centrifugal forces. For surface loading the faces of the elements are numbered as follows (for the node numbering of the elements see Section 3.1): for hexahedral elements: •	face 1: 1-2-3-4 •	face 2: 5-8-7-6 •	face 3: 1-5-6-2 •	face 4: 2-6-7-3 •	face 5: 3-7-8-4 •	face 6: 4-8-5-1 for tetrahedral elements: •	Face 1: 1-2-3 •	Face 2: 1-4-2 •	Face 3: 2-4-3 •	Face 4: 3-4-1 for wedge elements: •	Face 1: 1-2-3 •	Face 2: 4-5-6 •	Face 3: 1-2-5-4 •	Face 4: 2-3-6-5 •	Face 5: 3-1-4-6 for quadrilateral plane stress, plane strain and axisymmetric elements: •	Face 1: 1-2 •	Face 2: 2-3 •	Face 3: 3-4 •	Face 4: 4-1 for triangular plane stress, plane strain and axisymmetric elements: •	Face 1: 1-2 •	Face 2: 2-3 •	Face 3: 3-1 for beam elements: •	Face 1: pressure in 1-direction •	Face 2: pressure in 2-direction First line: •	*DLOAD •	Enter any needed parameters and their value Following line for surface loading: •	Element number or element set label. •	Distributed load type label. •	Actual magnitude of the load (for Px type labels) or fluid node number (for PxNU type labels) Repeat this line if needed. Example:

Se1,P3,10. assigns a pressure loading with magnitude 10. times the amplitude curve of amplitude A1 to face number three of all elements belonging to set Se1. PRINTING VARIABLES WITH NODE FILE: This option is used to print selected nodal variables in file jobname.frd for subsequent viewing by CalculiX GraphiX. The following variables can be selected: •	Displacements (key=U) •	Displacements: magnitude and phase (key=PU, only for *STEADY STATE DYNAMICS calculations) and *FREQUENCY calculations with cyclic symmetry). •	Maximum displacements orthogonal to a given vector at all times for *FREQUENCY calculations with cyclic symmetry. The components of the vector are the coordinates of a node stored in a node set with the name RAY. This node and node set must have been defined by the user (key=MAXU). •	Temperatures (key=NT). This includes both structural temperatures and total fluid temperatures in a network. •	Temperatures: magnitude and phase (key=PNT, only for *STEADY STATE DYNAMICS calculations). •	External forces (key=RF) •	External concentrated heat sources (key=RFL) •	Static temperatures in networks and 3D fluids (key=TS) •	Total temperatures in networks and 3D fluids (key=TT) •	Mass flows in networks (key=MF). The mass flow through a fluid element is stored in the middle node of the element. In the end nodes the mass flow is not unique, since more than two element can be connected to the node. For end nodes the sum of the mass flow leaving the node is stored. Notice that at nodes where mass flow is leaving the network the value will be wrong if no proper exit element (with node number 0) is attached to that node. •	Velocities in 3D fluids (key=V) •	Mach numbers in 3D compressible fluids (key=MACH) •	Static pressures in liquid networks and 3D fluids (key=PS) •	Total pressures in gas networks and 3D fluids (key=PT) •	Fluid depth in channel networks (key=DEPT) •	Critical depth in channel networks (key=HCRI) •	Pressure coefficient in 3D compressible fluids (key=CP) •	Turbulence variables in 3D compressible fluids (key=TURB)
 * DLOAD,AMPLITUDE=A1

First line: •	*NODE FILE •	Enter any needed parameters and their values. Second line: •	Identifying keys for the variables to be printed, separated by kommas. Example:

RF,NT
 * NODE FILE,FREQUENCY=2,TIME POINTS=T1

requests the storage of reaction forces and temperatures in the .frd file every second increment. In addition, output will be stored for all time points defined by the T1 time points sequence.

PRINTING VARIABLES WITH EL FILE: This option is used to save selected element variables averaged at the nodal points in a frd file (extension .frd) for subsequent viewing by CalculiX GraphiX. The following element variables can be selected: •	true (Cauchy) stress (key=S). For beam elements this tensor is replaced by the section forces if SECTION FORCES is selected. •	stress: magnitude and phase (key=PHS, only for *STEADY STATE DYNAMICS calculations and *FREQUENCY calculations with cyclic symmetry). •	maximum worst principal stress at all times for *FREQUENCY calculations with cyclic symmetry. It is stored for nodes belonging to the node set with name STRESSDOMAIN. This node set must have been defined by the user with the *NSET command. The worst principal stress is the maximum of the absolute value of the principal stresses (key=MAXS). •	strain (key=E). This is the Lagrangian strain for (hyper)elastic materials and incremental plasticity and the Eulerian strain for deformation plasticity. •	equivalent plastic strain (key=PEEQ) •	equivalent creep strain (key=CEEQ; is converted internally into PEEQ since the viscoplastic theory does not distinguish between the two; consequently, the user will find PEEQ in the frd file, not CEEQ) •	the energy density (key=ENER) •	the internal state variables (key=SDV) •	heat flux (key=HFL) •	Zienkiewicz-Zhu improved stress (key=ZZS) [66], [67] First line: •	*EL FILE •	Enter any needed parameters and their values. Second line: •	Identifying keys for the variables to be printed, separated by kommas. Example:

S,PEEQ requests that the (Cauchy) stresses and the equivalent plastic strain is stored in .frd format for subsequent viewing with CalculiX GraphiX
 * EL FILE

ENDING THE STEP: Type *END STEP at the end of the file A list of input deck formats can be found here http://www.bconverged.com/calculix/doc/ccx/html/index.html

SOLVING: In the same manner that you preprocess the .fbd file. The .inp file can be solved from the Tools Menu.

Samples are from CalculiX:

Steps to run the CCX:

1.Start the SciTe software.

2.Create a new file with intention of .inp.

3.Coding the CCX file in the pane(here we just paste the code copied from CalculiX), then Save

4.On the menu bar, select Tools--Solve

5.Again on the menu bar, select Tools--Post process

After doing step 5, a window will pop up.

'''6.Left click on the new interface, select Datasets. If you want to show the displacement plot, then select DISP. If you want to show the stress plot, select STRESS.'''

'''7.Again, Left click on the new interface, select Datasets. If you select DISP in step 6, then expand the Entity, you have four different choice of displacement to be displayed. If you select STRESS in the step 6, expand the Entity, you have more choice of stress to be displayed.'''

'''8.You also can alternate between DISP and STRESS. '''

Following are some sample result:



=Problem 5.6 Plots of Quadratic Lagrangian Element Basis Function (QLEBF)=

Problem Statement
Plot the Lagrangian Element basis function in case of a Quadratic element for the individual node and a combined plot.

Solution
A quadratic element has 3 nodes. Hence,

$$\;\;n = 3\;\;$$



For convenience we take

$$x_{1} = -1\; ; x_{2} = 0\; ;x_{3} = 1$$

The Lagrangian Interpolation Basis Function 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" |

L_{i,n}(x) = \prod_{k = 1}^{m}\frac{x-x_{k}}{x_{i}-x_{k}}\;\; k\neq i

$$
 * }

$$L_{1,3}(x) = \frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})} = \frac{(x-0)(x-1)}{(-1-0)(-1-1)} = \frac{x(x-1)}{2} = \frac{x^2-x}{2}$$

$$L_{2,3}(x) = \frac{(x-x_{1})(x-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})} = \frac{(x-(-1))(x-1)}{(0+1)(0-1)} = \frac{(x+1)(x-1)}{-1} = 1-x^2$$

$$L_{3,3}(x) = \frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})} = \frac{(x-(-1))(x-0)}{(1+1)(1-0)}  = \frac{(x+1)(x)}{2} = \frac{x^2+x}{2}$$

Using the above expressions, the plots of the basis functions were developed using Matlab









=Problem 5.7 Using Linear Lagrangian element basis functions to obtain the trial solution from Weak Form=

Problem statement
From 30-7:

Similar to HW 5.1 and HW 5.3, but using Linear Lagrangian element basis function with uniform discretization (equidistant element nodes, nel = 4, 6, 8, ...)

1.) Explain how Linear Lagrangian element basis function are used as CBS

2.) Plot all Lagrangian element basis function used

3.) Use matlab quad, WA, ... to integrate

4.) Plot $$ u_{nel}^h $$ versus $$ \displaystyle u $$ and $$\displaystyle u_{nel}^h(0.5)-u(0.5) $$ versus nel

Using LLEBF as CBS
i> Let the LLEBF be defined by a function


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

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

L_{i,2}(x) = \frac{x-x_{k}}{x_{i}-x_{k}}\;\; k\neq i $$ where n is the number of elements
 * }



Consider the number of elements =2, where each element is linear, i.e, each element will have only two nodes.

Now for the given function for element 1, LLEBF can be written as:-


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

$$\displaystyle L_{1,2}(x)= \frac{{x}-{{x}_{2}}}{{{x}_{1}}-{{x}_{2}}} $$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$\displaystyle L_{2,2}(x)= \frac{{x}-{{x}_{1}}}{{{x}_{2}}-{{x}_{1}}} $$


 * }

For element 2, LLEBF can be written as:-


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

$$\displaystyle L_{2,2}(x)= \frac{{x}-{{x}_{3}}}{{{x}_{2}}-{{x}_{3}}} $$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$\displaystyle L_{3,2}(x)= \frac{{x}-{{x}_{2}}}{{{x}_{3}}-{{x}_{2}}} $$


 * }

For element 1,


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

at $$\displaystyle x={{x}_{1}}, L_{1,2}(x)= \frac{{x}_{1}-{{x}_{2}}}{{{x}_{1}}-{{x}_{2}}} = 1$$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

and $$\displaystyle L_{2,2}(x)= \frac{{x}_{1}-{{x}_{1}}}{{{x}_{1}}-{{x}_{2}}}= 0$$

at $$\displaystyle x={{x}_{2}}, L_{1,2}(x)= \frac{{x}_{2}-{{x}_{2}}}{{{x}_{1}}-{{x}_{2}}}= 0$$

and $$\displaystyle L_{2,2}(x)= \frac{{x}_{2}-{{x}_{1}}}{{{x}_{2}}-{{x}_{1}}} = 1$$
 * }

For element 2,


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

at $$\displaystyle x={{x}_{2}}, L_{2,2}(x)= \frac{{x}_{2}-{{x}_{3}}}{{{x}_{2}}-{{x}_{3}}} = 1$$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

and $$\displaystyle L_{3,2}(x)= \frac{{x}_{2}-{{x}_{2}}}{{{x}_{3}}-{{x}_{2}}}= 0$$

at $$\displaystyle x={{x}_{3}}, L_{1,2}(x)= \frac{{x}_{3}-{{x}_{3}}}{{{x}_{1}}-{{x}_{2}}}= 0$$

and $$\displaystyle L_{3,2}(x)= \frac{{x}_{3}-{{x}_{2}}}{{{x}_{3}}-{{x}_{2}}} = 1$$ In element 1, when $$\displaystyle x = x_{1}$$, the value of the basis function is 1 at node 1 and zero at node 2. Similarly, when $$\displaystyle x = x_{2}$$ the value of the basis function is 1 at node 2 and zero at node 1., which is required condition for any function to be the CBS
 * }

Hence the LLEBF can be used as a CBS as it satisfies the conditions for the CBS.

Plotting the LLEBF used
ii> Liner Lagrangian Interpolation Basis Function plotted for 4 elements is shown as below.

It is clear from the figure that for element 1 at node 1, the value of the basis function is 1 and 0 at node 2, which is in accordance with the CBS.

Same is observed for rest of the elements.



Using Matlab to integrate the basis functions
iii> Using the above expressions, the plots of the basis functions were developed using Matlab

Plotting the LIBF and Exact Solution and the errors at x=0.5 for different values of n
iv> The convergence occurs at $$\displaystyle n=59 $$, the numbers of elements is found to be a high value since, the convergence is checked till the tolerance reaches $$10^{-4}$$,





=Problem 5.8 Using Linear Lagrangian element basis functions to obtain the trial solution from Weak Form=

Problem statement
From 30-7:

Similar to HW 5.2 and HW 5.4, but using Linear Lagrangian element basis function with uniform discretization (equidistant element nodes, nel = 4, 6, 8, ...)

1.) Explain how Linear Lagrangian element basis function are used as CBS

2.) Plot all Lagrangian element basis function used

3.) Use matlab quad, WA, ... to integrate

4.) Plot $$ u_{nel}^h $$ versus $$ \displaystyle u $$ and $$\displaystyle u_{nel}^h(0.5)-u(0.5) $$ versus nel

Using LLEBF as CBS
i> Let the LLEBF be defined by a function


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

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

L_{i,2}(x) = \frac{x-x_{k}}{x_{i}-x_{k}}\;\; k\neq i $$ where n is the number of elements
 * }



Consider the number of elements =2, where each element is linear, i.e, each element will have only two nodes.

Now for the given function for element 1, LLEBF can be written as:-


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

$$\displaystyle L_{1,2}(x)= \frac{{x}-{{x}_{2}}}{{{x}_{1}}-{{x}_{2}}} $$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$\displaystyle L_{2,2}(x)= \frac{{x}-{{x}_{1}}}{{{x}_{2}}-{{x}_{1}}} $$


 * }

For element 2, LLEBF can be written as:-


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

$$\displaystyle L_{2,2}(x)= \frac{{x}-{{x}_{3}}}{{{x}_{2}}-{{x}_{3}}} $$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$\displaystyle L_{3,2}(x)= \frac{{x}-{{x}_{2}}}{{{x}_{3}}-{{x}_{2}}} $$


 * }

For element 1,


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

at $$\displaystyle x={{x}_{1}}, L_{1,2}(x)= \frac{{x}_{1}-{{x}_{2}}}{{{x}_{1}}-{{x}_{2}}} = 1$$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

and $$\displaystyle L_{2,2}(x)= \frac{{x}_{1}-{{x}_{1}}}{{{x}_{1}}-{{x}_{2}}}= 0$$

at $$\displaystyle x={{x}_{2}}, L_{1,2}(x)= \frac{{x}_{2}-{{x}_{2}}}{{{x}_{1}}-{{x}_{2}}}= 0$$

and $$\displaystyle L_{2,2}(x)= \frac{{x}_{2}-{{x}_{1}}}{{{x}_{2}}-{{x}_{1}}} = 1$$
 * }

For element 2,


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

at $$\displaystyle x={{x}_{2}}, L_{2,2}(x)= \frac{{x}_{2}-{{x}_{3}}}{{{x}_{2}}-{{x}_{3}}} = 1$$
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

and $$\displaystyle L_{3,2}(x)= \frac{{x}_{2}-{{x}_{2}}}{{{x}_{3}}-{{x}_{2}}}= 0$$

at $$\displaystyle x={{x}_{3}}, L_{1,2}(x)= \frac{{x}_{3}-{{x}_{3}}}{{{x}_{1}}-{{x}_{2}}}= 0$$

and $$\displaystyle L_{3,2}(x)= \frac{{x}_{3}-{{x}_{2}}}{{{x}_{3}}-{{x}_{2}}} = 1$$
 * }

In element 1, when $$\displaystyle x = x_{1}$$, the value of the basis function is 1 at node 1 and zero at node 2. Similarly, when $$\displaystyle x = x_{2}$$ the value of the basis function is 1 at node 2 and zero at node 1., which is required condition for any function to be the CBS

Hence the LLEBF can be used as a CBS as it satisfies the conditions for the CBS.

Plotting the LLEBF used
ii> Liner Lagrangian Interpolation Basis Function plotted for 4 elements is shown as below.

It is clear from the figure that for element 1 at node 1, the value of the basis function is 1 and 0 at node 2, which is in accordance with the CBS.

Same is observed for rest of the elements.



Using Matlab to integrate the basis functions
iii> Using the above expressions, the plots of the basis functions were developed using Matlab

Plotting the LIBF and Exact Solution and the errors at x=0.5 for different values of n
iv> The convergence occurs at $$\displaystyle n=23 $$





=Contributing Members & Referenced Lecture=