University of Florida/Eml5526/s11.team2.reiss.HW3 HW3

Obtaining Solution to the Weak Form from a Trial Solution

Given

${\frac {d}{dx}}\left({AE{\frac {du}{dx}}}\right)+2x=0$ on $1<x<3$

(10.1)

$\sigma (1)={\left({E{\frac {du}{dx}}}\right)_{x=1}}=0.1$

(10.2)

$u(3)=0.001$

(10.3)

where the weak form is

$\int \limits _{1}^{3}{{\frac {dw}{dx}}AE{\frac {du}{dx}}dx=-0.1{{(wA)}_{x=1}}+\int \limits _{1}^{3}{2xwdx}}$ $\forall w$ with $w(3)=0$

(10.4)

A trial solution for u(x) is suggested as

$u(x)={\alpha _{0}}+{\alpha _{1}}(x-3)+{\alpha _{2}}{(x-3)^{2}}$

(10.5)

with w(x) of the same form.

1) determine the solution to the weak form. 2) verify that the solution satisfies the equilibrium equation of the strong form 3) verify that the natural boundary condition is satisfied

Solution 3.10.1

From the given information the equations for u(x) and w(x) are as follows

$u(x)={\alpha _{0}}+{\alpha _{1}}(x-3)+{\alpha _{2}}{(x-3)^{2}}$

(10.5)

$w(x)={\beta _{0}}+{\beta _{1}}(x-3)+{\beta _{2}}{(x-3)^{2}}$

(10.6)

Using the initial condition $w(3)=0$

$w(3)={\beta _{0}}+{\beta _{1}}(3-3)+{\beta _{2}}{(3-3)^{2}}\Rightarrow w(3)=0={\beta _{0}}$

(10.7)

Using the condition that $u(3)=0.001$

$w(3)={\alpha _{0}}+{\alpha _{1}}(3-3)+{\alpha _{2}}{(3-3)^{2}}\Rightarrow u(3)=0.001={\alpha _{0}}$

(10.8)

The partial derivatives of u(x) and w(x) are

${\frac {du(x)}{dx}}={\alpha _{1}}+2{\alpha _{2}}(x-3)$

(10.8)

${\frac {dw(x)}{dx}}={\beta _{1}}+2{\beta _{2}}(x-3)$

(10.9)

Replacing, the weak form becomes

$AE\int \limits _{1}^{3}{\left[{{\beta _{1}}+2{\beta _{2}}(x-3)}\right]}\left[{{\alpha _{1}}+2{\alpha _{2}}(x-3)}\right]dx=-0.1({\beta _{1}}-4{\beta _{2}})A+\int \limits _{1}^{3}{2x\left[{{\beta _{1}}+2{\beta _{2}}(x-3)}\right]}dx$

(10.10)

After integration the equation becomes

$AE[{\beta _{1}}(2{\alpha _{1}}-4{\alpha _{2}})+{\beta _{2}}(-4{\alpha _{1}}+{\frac {10}{3}}{\alpha _{2}})]=-0.1({\beta _{1}}-4{\beta _{2}})A+8{\beta _{1}}-{\frac {62}{3}}{\beta _{2}}$

(10.11)

Rearranging we get

${\beta _{1}}\left[{AE(2{\alpha _{1}}-4{\alpha _{2}})+0.1A-8}\right]+{\beta _{2}}\left[{AE(-4{\alpha _{1}}+{\frac {10}{3}}{\alpha _{2}})-.4A+{\frac {62}{3}}}\right]=0$

(10.12)

The $\beta _{0}$ and $\beta _{1}$ values must be arbitrary, the solutions to this equation are dependent on the bracketed terms equaling zero. The solution for this, in matrix form, is

$AE\left[{\begin{array}{*{20}{c}}2&{-4}\\{-4}&{\frac {10}{3}}\end{array}}\right]\left[{\begin{array}{*{20}{c}}{\alpha _{1}}\\{\alpha _{2}}\end{array}}\right]=\left[{\begin{array}{*{20}{c}}{0.1A-8}\\{.4A-20.666}\end{array}}\right]$

(10.13)

The solution to the weak form becomes

$\left[{\begin{array}{*{20}{c}}{\alpha _{1}}\\{\alpha _{2}}\end{array}}\right]={\left[{\begin{array}{*{20}{c}}2&{-4}\\{-4}&{\frac {10}{3}}\end{array}}\right]^{-1}}\left[{\begin{array}{*{20}{c}}{0.1A-8}\\{0.4A-20.666}\end{array}}\right]{(AE)^{-1}}$

(10.14)

Giving the solution of

$\left[{\begin{array}{*{20}{c}}{\alpha _{1}}\\{\alpha _{2}}\end{array}}\right]=\left[{\begin{array}{*{20}{c}}{-.2071+{\frac {11.7132}{a}}}\\{-{\rm {0}}{\rm {.1286}}+{\frac {{\rm {7}}{\rm {.8566}}}{a}}}\end{array}}\right]{E^{-1}}$

(10.15)

Problem 3.7: Rework Problem 3.1 without using a shift in the given polynomial basis.

Given

$\displaystyle 2{\frac {d^{2}u}{dx^{2}}}+3=0$

(7.1)

Boundary Conditions:

$\displaystyle u(x=1)=0$

(7.2)

$\displaystyle -{\frac {du}{dx}}(x=0)=4$

(7.3)

$\displaystyle b_{j}=(x)^{j}$

(7.4)

The exact solution to Eq.3.1:

$\displaystyle u(x)=-{\frac {3}{4}}x^{2}-4x+{\frac {19}{4}}$

(7.5)

Find

We want to find an approximation to Eq.3.5 in the form of

$\displaystyle u^{h}(x)=\sum _{i=0}^{n}a_{i}b_{i}$

(7.6)

as well as investigate how the accuracy of our approximation changes with n going from 2 to 4 and then 6.

Case:n=2

First let's expand Eq.3.6 out to n=2 terms.

$\displaystyle u^{h}(x)=a_{0}b_{0}+a_{1}b_{1}+a_{2}b_{2}$

(7.7)

Now consider Eq.7.4, and substitute the appropriate terms into Eq.7.7.

$\displaystyle u^{h}(x)=a_{0}(1)+a_{1}(x)+a_{2}(x)^{2}$

(7.8)

Let's rewrite this equation in terms of a product between 2 vectors.

$\displaystyle u^{h}(x)=\left[1\quad x\quad (x)^{2}\right]{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}$

(7.9)

In order for us to solve for the unknown $a_{i}'s$ in Eq.7.9, we must have 3 linearly independent equations describing linear combinations of the $a_{i}'s.$Two of the equations are easily obtained directly from the given boundary conditions in Eq.7.2 and Eq.7.3. It only makes practical sense that if we want our approximation to closely resemble our exact solution, we must require identical behavior between Eq.7.5 and Eq.7.8 at the boundaries.

$\displaystyle {\begin{aligned}u^{h}(x=1)=a_{0}(1)+a_{1}(1)+a_{2}(1)^{2}&=0\\a_{0}(1)+a_{1}(1)+a_{2}(1)^{2}&=0\end{aligned}}$

(7.10)

$\displaystyle a_{0}+a_{1}+a_{2}=0$

(7.11)

$\displaystyle {\frac {du}{dx}}(x=0)=a_{1}+a_{2}2(0)=-4$

(7.12)

$\displaystyle a_{1}=-4$

(7.13)

Now we just need one more linearly independent equation and we can then readily solve for the $a_{i}'s.$ To develop another set of equations which are linearly independent of Eq.7.11 and Eq.7.13, we take the inner product between our approximation of u(x), Eq.7.9 and our basis set of functions Eq.7.4. Before getting to the inner product first, let's define our linear differential operator as the following

$\displaystyle P(\cdot )={\frac {d^{2}}{dx^{2}}}+{\frac {3}{2}}=0$

(7.14)

Now, using our operator in Eq.7.14, we can now operate on our approximation in Eq.7.9.

$\displaystyle {\begin{aligned}P(u^{h})={\frac {d^{2}}{dx^{2}}}\left[1\quad x\quad (x)^{2}\right]{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}+{\frac {3}{2}}&=0\\\left[0\quad 0\quad 2\right]{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}+{\frac {3}{2}}&=0\end{aligned}}$

(7.15)

Which we can then rewrite Eq.7.15 as

$\displaystyle \left[0\quad 0\quad 2\right]{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}=-{\frac {3}{2}}$

(7.16)

Therefor

$2a_{2}=-{\frac {3}{2}}$

(7.21)

Using Eq.7.11 and Eq.7.13, we can form the following matrix equation.

$\displaystyle {\begin{bmatrix}1&1&1\\0&1&0\\0&0&2\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}={\begin{bmatrix}0\\-4\\-3/2\end{bmatrix}}$

(7.21)

Let's solve for the inverse of the matrix in Eq.7.21 using the Gauss-Jordan Method. Letting the matrix above in Eq.7.21 equal $\mathbf {K}$, we can find $\mathbf {K^{-1}}$ by

$\displaystyle {\begin{aligned}\mathbf {K} \cdot \mathbf {I} &={\begin{bmatrix}1&1&1&1&0&0\\0&1&0&0&1&0\\0&0&2&0&0&1\end{bmatrix}}\\&={\begin{bmatrix}1&0&0&1&-1&-1/2\\0&1&0&0&1&0\\0&0&2&0&0&1\end{bmatrix}}\\&={\begin{bmatrix}1&0&0&1&-1&-1/2\\0&1&0&0&1&0\\0&0&1&0&0&1/2\end{bmatrix}}\\\end{aligned}}$

(7.22)

$\displaystyle \therefore \mathbf {K^{-1}} ={\begin{bmatrix}1&-1&-1/2\\0&1&0\\0&0&1/2\end{bmatrix}}$

(7.23)

Now using Eq.7.23, we can solve Eq.7.21 for our unknown $a_{i}'s.$

$\displaystyle \mathbf {K^{-1}} \cdot \mathbf {K} {\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}=\mathbf {K^{-1}} {\begin{bmatrix}0\\-4\\-3/2\end{bmatrix}}$

(7.24)

$\displaystyle {\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}={\begin{bmatrix}1&-1&-1/2\\0&1&0\\0&0&1/2\end{bmatrix}}{\begin{bmatrix}0\\-4\\-3/2\end{bmatrix}}$

(7.25)

$\displaystyle {\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\end{bmatrix}}={\begin{bmatrix}19/4\\-4\\-3/4\end{bmatrix}}$

(7.26)

Solution

$\displaystyle \therefore u^{h}(x)={\frac {19}{4}}-4(x)-{\frac {3}{4}}(x)^{2}$

(7.27)

Now. let's check the absolute error between our approximate solution above in Eq. 7.26 and our exact solution in Eq.7.5.

Comparing $u^{h}(x)$, 7.27, to $u(x)$, 7.5, we can see that using this basis gives an exact solution to $u(x)$

$\displaystyle u^{h}(x)={\frac {19}{4}}-4(x)-{\frac {3}{4}}(x)^{2}$

(7.27)

$\displaystyle u(x)={\frac {19}{4}}-4(x)-{\frac {3}{4}}(x)^{2}$

(7.5)

Graphical Comparison Between Exact and Numerical Approximate Solutions

Case:n=4

We will follow exactly the same methodology as we have above, leaving out the explanations for each step.

$\displaystyle u^{h}(x)=a_{0}b_{0}+a_{1}b_{1}+a_{2}b_{2}+a_{3}b_{3}+a_{4}b_{4}$

(7.32)

$\displaystyle u^{h}(x)=a_{0}(1)+a_{1}(x)+a_{2}(x)^{2}+a_{3}(x)^{3}+a_{4}(x)^{4}$

(7.33)

$\displaystyle u^{h}(x)={\begin{bmatrix}1&x&x^{2}&x^{3}&x^{4}\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{bmatrix}}$

(7.34)

Using our given boundary conditions to obtain our first 2 linearly independent equations.

$\displaystyle {\begin{aligned}u^{h}(x=1)&=a_{0}(1)+a_{1}(1)+a_{2}(1)^{2}+a_{3}(1)^{3}+a_{4}(1)^{4}\\&=a_{0}+a_{1}+a_{2}+a_{3}+a_{4}\end{aligned}}$

(7.35)

$\displaystyle a_{0}+a_{1}+a_{2}+a_{3}+a_{4}=0$

(7.36)

$\displaystyle {\frac {du^{h}}{dx}}(x=0)=0+a_{1}+0+0+0$

(7.37)

$\displaystyle a_{1}=-4$

(7.38)

Using our predefined linear differential operator in Eq.7.14.

$\displaystyle P(u^{h})={\frac {d^{2}}{dx^{2}}}{\begin{bmatrix}1&x&x^{2}&x^{3}&x^{4}\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{bmatrix}}+{\frac {3}{2}}$

(7.39)

Because this time we have 5 unknown $a_{i}'s,$ we will have to choose 3 equations from our inner product operation.

$\displaystyle \left\langle b_{i},P(u^{h})\right\rangle =\int _{0}^{1}b_{i}P(u^{h})=0$

(7.40)

$\displaystyle \int _{0}^{1}{\begin{bmatrix}1\\x\\x^{2}\\x^{3}\\x^{4}\end{bmatrix}}{\begin{bmatrix}0&0&2&6(x)&12(x)^{2}\end{bmatrix}}\,dx{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{bmatrix}}=-{\frac {3}{2}}\int _{0}^{1}{\begin{bmatrix}1\\x\\x^{2}\\x^{3}\\x^{4}\end{bmatrix}}\,dx$

(7.41)

$\displaystyle \int _{0}^{1}{\begin{bmatrix}0&0&2&6(x)&12(x)^{2}\\0&0&2(x)&6(x)^{2}&12(x)^{3}\\0&0&2(x)^{2}&6(x)^{3}&12(x)^{4}\\0&0&2(x)^{3}&6(x)^{4}&12(x)^{5}\\0&0&2(x)^{4}&6(x)^{5}&12(x)^{6}\end{bmatrix}}\,dx{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{bmatrix}}=-{\frac {3}{2}}\int _{0}^{1}{\begin{bmatrix}1\\x\\x^{2}\\x^{3}\\x^{4}\end{bmatrix}}\,dx$

(7.42)

$\displaystyle {\begin{bmatrix}0&0&2&3&4\\0&0&1&2&3\\0&0&2/3&3/2&12/5\\0&0&1/2&6/5&2\\0&0&2/5&1&12/7\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{bmatrix}}=-{\begin{bmatrix}3/2\\3/4\\3/6\\3/8\\3/10\end{bmatrix}}$

(7.43)

Now we can use the bottom three rows of Eq.7.43, along with Eq.7.36 and Eq.7.38 to develop a system of equations

$\displaystyle {\begin{bmatrix}1&1&1&1&1\\0&1&0&0&0\\0&0&2/3&3/2&12/5\\0&0&1/2&6/5&2\\0&0&2/5&1&12/7\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{bmatrix}}=-{\begin{bmatrix}0\\4\\3/6\\3/8\\3/10\end{bmatrix}}$

(7.44)

Solving this system using MATLAB, the solution becomes

$\displaystyle {\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{bmatrix}}={\begin{bmatrix}19/4\\-4\\-3/4\\0\\0\end{bmatrix}}$

(7.45)

Solution

$\displaystyle \therefore u^{h}(x)={\frac {19}{4}}-4(x)-{\frac {3}{4}}(x)^{2}$

(7.46)

We can see from above in Eq.7.46 that our solution did not improve at all by an increase in n from 2 to 4. This should have been intuitive from the start since we showed earlier that using n=2 gives us an equivalent $u(x)$

Graphical Comparison Between Exact and Numerical Approximate Solutions

Egm5526.s11.team-2.langpm 14:38, 14 February 2011 (UTC)

Case:n=6

$\displaystyle {\begin{aligned}u^{h}(x)&=a_{0}b_{0}+a_{1}b_{1}+a_{2}b_{2}+a_{3}b_{3}+a_{4}b_{4}+a_{5}b_{5}+a_{6}b_{6}\\&=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+a_{4}x^{4}+a_{5}x^{5}+a_{6}x^{6}\end{aligned}}$

(7.47)

Using the given boundary conditions in Eq.7.2 and Eq.7.3 we have the following to restrictions on our approximated solution.

$\displaystyle {\begin{aligned}u^{h}(x=1)=&a_{0}+a_{1}+a_{2}+a_{3}+a_{4}+a_{5}+a_{6}=0\\&a_{0}+a_{1}+a_{2}+a_{3}+a_{4}+a_{5}+a_{6}=0\end{aligned}}$

(7.48)

$\displaystyle a_{0}+a_{1}+a_{2}+a_{3}+a_{4}+a_{5}+a_{6}=0$

(7.49)

$\displaystyle {\frac {du^{h}}{dx}}(x=0)=0+a_{1}+0+0+0+0+0$

(7.50)

$\displaystyle a_{1}=-4$

(7.51)

Using the same linear differential operator as we have for the past 2 cases, Eq.7.14 we have the following.

$\displaystyle P(u^{h})={\frac {d^{2}}{dx^{2}}}{\begin{bmatrix}1&x&x^{2}&x^{3}&x^{4}&x^{5}&x^{6}\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\a_{5}\\a_{6}\end{bmatrix}}+{\frac {3}{2}}$

(7.52)

Forming the inner product.

$\displaystyle \left\langle b_{i},P(u^{h})\right\rangle =\int _{0}^{1}b_{i}P(u^{h})\,dx=0$

(7.53)

$\displaystyle \int _{0}^{1}{\begin{bmatrix}0&0&2&6(x)&12(x)^{2}&20(x)^{3}&30(x)^{4}\\0&0&2(x)&6(x)^{2}&12(x)^{3}&20(x)^{4}&30(x)^{5}\\0&0&2(x)^{2}&6(x)^{3}&12(x)^{4}&20(x)^{5}&30(x)^{6}\\0&0&2(x)^{3}&6(x)^{4}&12(x)^{5}&20(x)^{6}&30(x)^{7}\\0&0&2(x)^{4}&6(x)^{5}&12(x)^{6}&20(x)^{7}&30(x)^{8}\\0&0&2(x)^{5}&6(x)^{6}&12(x)^{7}&20(x)^{8}&30(x)^{9}\\0&0&2(x)^{6}&6(x)^{7}&12(x)^{8}&20(x)^{9}&30(x)^{10}\end{bmatrix}}\,dx{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\a_{5}\\a_{6}\end{bmatrix}}=-{\frac {3}{2}}\int _{0}^{1}{\begin{bmatrix}1\\x\\(x)^{2}\\(x)^{3}\\(x)^{4}\\(x)^{5}\\(x)^{6}\end{bmatrix}}$

(7.54)

$\displaystyle {\begin{bmatrix}0&0&2&3&4&5&6\\0&0&1&2&3&4&5\\0&0&2/3&3/2&12/5&10/3&30/7\\0&0&1/2&6/5&2&20/7&15/2\\0&0&2/5&1&12/7&5/2&10/3\\0&0&1/3&6/7&3/2&20/9&3\\0&0&2/7&3/4&4/3&1/2&30/11\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\a_{5}\\a_{6}\end{bmatrix}}={\begin{bmatrix}-3/2\\-3/4\\-3/6\\-3/8\\-3/10\\-3/12\\-3/14\end{bmatrix}}$

(7.55)

We take the bottom 5 rows from Eq.7.55 and build a system of equations with Eq.7.49 and Eq.7.51.

$\displaystyle {\begin{bmatrix}1&1&1&1&1&1&1\\0&1&0&0&0&0&0\\0&0&2/3&3/2&12/5&10/3&30/7\\0&0&1/2&6/5&2&20/7&15/2\\0&0&2/5&1&12/7&5/2&10/3\\0&0&1/3&6/7&3/2&20/9&3\\0&0&2/7&3/4&4/3&1/2&30/11\end{bmatrix}}{\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\a_{5}\\a_{6}\end{bmatrix}}={\begin{bmatrix}0\\-4\\-3/6\\-3/8\\-3/10\\-3/12\\-3/14\end{bmatrix}}$

(7.56)

As in the $n=4$ case, we use a MATLAB to determine the unknown $a_{i}'s$ adn to determine the final form of our approximation to the solution of the differential equation 7.1.

$\displaystyle {\begin{bmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\a_{5}\\a_{6}\end{bmatrix}}={\begin{bmatrix}19/4\\-4\\-3/4\\0\\0\\0\\0\end{bmatrix}}$

(7.57)

Solution

$\displaystyle \therefore u^{h}(x)={\frac {19}{4}}-4(x)-{\frac {3}{4}}(x)^{2}$

(7.58)

We can see from above in Eq.7.58 that our solution did not improve at all by an increase in n from 2 to 6 as was expected. Comparing this to the results from problem 3.1 we can see that the shift is not necessary for this basis function. It will return an exact solution with or without the shift.

Graphical Comparison Between the Exact and Numerical Solutions