User:Eml5526.s11.team3.akj/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 $$, $$ \forall x \in \left[ 0,1 \right] $$, $$ u(0)=4 $$, $$ 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 = [0] $$, $$ \Gamma_h = [1] $$

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(0.5)$$ to $$10^{-6}$$.

Plot $$ u, 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=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=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}}=\{{{x}^{j}},j=0,1,2,\cdots \}

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

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

{{f}_{F}}=\{\cos (jx)-1,\sin (kx),j=0,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}}=\{{{e}^{jx}}-1,j=0,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}}{{x}^{0}}+{{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(x=1)=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}}[\cos [0*(x)]-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(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)(2+3x)j\cos j(x)}dx,_ – ^ – when_ – ^ – i,j_ – ^ – both_ – ^ – are_ – ^ – even_ – ^ – number \\ \int\limits_{0}^{1}{i\sin i(x)(2+3x)j\sin j(x)}dx,_ – ^ – when_ – ^ – i,j_ – ^ – both_ – ^ – are_ – ^ – odd_ – ^ – number \\ \int\limits_{0}^{1}{i\cos i(x)(2+3x)j\sin j(x)}dx,_ – ^ – 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)}dx-(12+18*0)\sin i(1),_ – ^ – when_ – ^ – i_ – ^ – is_ – ^ – even_ – ^ – number \\ \int\limits_{0}^{1}{5x[\cos i(x)-1]}dx-(12+18*0)[\cos i(1)-1],_ – ^ – 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} 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}}({{e}^{0*(x)}}-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(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)}}(2+3x)j{{e}^{j(x)}}}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]dx}-(12+18*0)[{{e}^{i(1)}}-1]]

$$
 * }

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.