# Problem 5.1: Deriving the Weak Form From the Strong Form (Data 1b)

## Given: Strong Form, Boundary Conditions, and Basis Functions

The Strong Form from lecture 9-1 and Fish and Belytschko p. 73:

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (1.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (1.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (1.3)

 ${\Gamma _g} = \left\{ 0 \right\} \to u(x = 0) = 4$ (1.4)

 ${\Gamma _g} = \left\{ 1 \right\} \to \frac{{du}}{{dx}}\left( {x = 1} \right) = - 6$ (1.5)

## Find: The Weak Form in Order to Solve for $u^h(x)$ with Following Basis Functions

Solve for $\displaystyle u(x)$ and $\displaystyle u^h(x)$ using the following basis functions until the error $\displaystyle e(n) =|u^h(0.5)-u(0.5)| \leq 10^{-6}$.

1. Polynomial Basis Function $\displaystyle \mathcal{F}_P :=\{ x^j,j=0,1,2,...\}$
2. Exponential Basis Function $\displaystyle \mathcal{F}_e :=\{ e^{jx},j=0,1,2,...\}$
3. Fourier Basis Function $\qquad \displaystyle \mathcal{F}_F :=\{cos(jx),j=0,1,2,...,sin(kx),k=1,2,...\}$

On the domain defined by,

 $\Omega = \left] {\alpha ,\beta } \right[ = \left] { 0,1} \right[$

## Solution: Using Boundary Conditions and Constraint Breaking Solutions

### 1. Approximating Solution Using Polynomial Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

 $\frac{\partial}{\partial x} \left[ {\left( {{\text{2}} + {\text{3x}}} \right) \frac{{\partial u}}{{\partial x}}} \right] + 5x = 0$ (1.21)

Multiplying Equ 2.20 by w(x) and integrating within the domain:

 $\displaystyle \int_0^{1} {w(x) \cdot \left[\frac{\partial}{\partial x}\left [ (2+3x)\frac{\partial u}{\partial x} \right ]+5x\right] } dx = 0$ (1.22)

Next we will integrate by parts to obtain,

 $\displaystyle \left. {w(x)(2+3x)\frac{\partial u}{\partial x}} \right|_0^{1} - \int_0^{1} {\frac{\partial w}{\partial x}(2+3x)\frac{\partial u}{\partial x}} dx + \int_0^{1} 5wx dx= 0$ (1.23)

Since $\displaystyle w(x=0) = 0$ we can simplify (2.22) in the following way to obtain the weak form,

 $-30w(1) - \int_0^1 {\left[ {\frac{{dw}}{{dx}}\left( {2 + 3x} \right)\frac{{du}}{{dx}}} \right]} dx + \int_0^1 {5xw(x)} dx = 0$ (1.24)

Now we can use the polynomial basis functions,

 $w(x) = {b_i}{(x - 1)^i} \to \frac{{dw}}{{dx}} = i{b_i}{(x - 1)^{i - 1}}$ (1.25)
 $\displaystyle u(x) = {b_j}{(x - 1)^j} \to \frac{{du}}{{dx}}j{b_j}{(x - 1)^{j - 1}}$ (1.26)
 $\displaystyle such \quad that \quad w(0) = 0$ (1.27)

Which allows us to write the weak form (2.25) in terms of our matrices $\displaystyle {\color{red}\bold{K}_{FF}}$, $\displaystyle {\color{red}\bold{d}_{F}}$ and $\displaystyle {\color{red}\bold{F}_{F}}$.

 $\underbrace{\left[ \int_0^1 {{\frac{{\partial b_i(x)}}{{\partial x}} (2 + 3x)\frac{{\partial b_j(x)}}{{\partial x}}} } dx \right]}_{{\color{red}\bold{K}_{FF}}} \underbrace{d_i}_{{\color{red}\bold{d}_{F}}} = \underbrace{\left[ \int_0^1 {5xb_i(x)} dx - 30b_i(x)|_{x=1} \right]}_{{\color{red}\bold{F}_{F}}}$ (1.28)

Next we wish to obtain the solution to $\displaystyle \bold{d_F}$ by solving the following system of equations,

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}} \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$ (1.29)

Solving for the approximating solution, the d matrix can transposed and multiplied by the basis function matrix:

 $U_p^h = {d_0}{b_0} + {d_1}{b_1} + {d_2}{b_2}$ (1.30)

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=8

We found that for $\displaystyle n=8$ the error was of the order $\displaystyle 10^{-6}$. Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\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)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} 41/3 \\ 53/4 \\ 13 \\ 77/6 \\ 89/7 \\ 101/8 \\ 113/9 \\ 25/2\\ \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.0000000 \\ & 7.2497811 \\ & -5.4294546 \\ & 4.9208409 \\ & -5.0302374 \\ & 4.5191594 \\ & -3.0192656 \\ & 1.2522298 \\ & -0.23479307 \\ \end{matrix} \right)}}$

 $U_p^h(x) = 4 - 7.2497811(x - 1) + -5.4294546{(x - 1)^2} +4.9208409{(x - 1)^3} + -5.0302374{(x - 1)^4} + 4.5191594{(x - 1)^5} + -3.0192656{(x - 1)^6} 1.2522298 (x - 1)^7-0.23479307(x - 1)^8$

(1.31)

### Matlab Code

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

#### Plot: Showing the error $e(n) = |u(x)-u_h(x)|$

For the polynomial approximate solution, with the increasing number of nodes, the error decreases. For this particular instance, the error XXXXXXX

### 2. Approximating Soultion Using Fourier Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the fourier solution is obtained the same way. W(x) is now equation 5.35 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\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)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} 4 & -6.1075 & 11.6035 & -18.9907 & 13.0886 & -27.2503 & 3.4218 & -23.8065 & -8.5011 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} 4.0000 & 47.4091 & 48.5755 & -3.4993 & -30.6212 & -4.1185 & 7.2053 & 0.9062 & -0.4248 \\ \end{matrix} \right)^T}}$

### Matlab Code

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

### 3. Approximating Soultion Using Exponential Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the exponential solution is obtained the same way. W(x) is now equation 5.41 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}= 10^8\left( \begin{matrix} 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0001 & 0.0003 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0001 & 0.0002 & 0.0006 & 0.0017\\ 0.0000 & 0.0000 & 0.0000 & 0.0001 & 0.0003 & 0.0008 & 0.0022 & 0.0062\\ 0.0000 & 0.0000 & 0.0001 & 0.0003 & 0.0008 & 0.0025 & 0.0072 & 0.0206\\ 0.0000 & 0.0001 & 0.0003 & 0.0008 & 0.0026 & 0.0077 & 0.0225 & 0.0649\\ 0.0000 & 0.0002 & 0.0008 & 0.0025 & 0.0077 & 0.0232 & 0.0682 & 0.1973\\ 0.0001 & 0.0006 & 0.0022 & 0.0072 & 0.0225 & 0.0682 & 0.2014 & 0.5858\\ 0.0003 & 0.0017 & 0.0062 & 0.0206 & 0.0649 & 0.1973 & 0.5858 & 1.7106\\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{10^4\left( \begin{matrix} 4 & 0.0023 & 0.0085 & 0.0249 & 0.0692 & 0.1885 & 0.5107 & 1.3817 & 3.7387 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.000000000000 \\ & -416.206015822685 \\ & 1802.371470713180 \\ & -4682.097146605770 \\ & 7738.463302468000 \\ & -8239.102282958180 \\ & 5487.180599937650 \\ & -2083.316654302480 \\ & 344.581494129590 \\ \end{matrix} \right)}}$

### Matlab Code

#### Plots

The plots for solutions and error at x=.5 are as follows:

# Problem 5.2: Deriving the Weak Form From the Strong Form (Data 1)

## Given: Strong Form, Boundary Conditions, and Basis Functions

The Strong Form from lecture 9-1 and Fish and Belytschko p. 73:

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (2.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (2.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (2.3)

 ${\Gamma _g} = \left\{ 1 \right\} \to u(x = 1) = 4$ (2.4)

 ${\Gamma _g} = \left\{ 0 \right\} \to \frac{{du}}{{dx}}\left( {x = 0} \right) = - 6$ (2.5)

## Find: The Weak Form in Order to Solve for $u^h(x)$ with Following Basis Functions

Solve for $\displaystyle u(x)$ and $\displaystyle u^h(x)$ using the following basis functions until the error $\displaystyle e(n) =|u^h(0.5)-u(0.5)| \leq 10^{-6}$.

1. Polynomial Basis Function $\displaystyle \mathcal{F}_P :=\{ x^j,j=0,1,2,...\}$
2. Exponential Basis Function $\displaystyle \mathcal{F}_e :=\{ e^{jx},j=0,1,2,...\}$
3. Fourier Basis Function $\qquad \displaystyle \mathcal{F}_F :=\{cos(jx),j=0,1,2,...,sin(kx),k=1,2,...\}$

On the domain defined by,

 $\Omega = \left] {\alpha ,\beta } \right[ = \left] { 0,1} \right[$

## Solution: Using Boundary Conditions and Constraint Breaking Solutions

### 1. Approximating Solution Using Polynomial Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

 $\frac{\partial}{\partial x} \left[ {\left( {{\text{2}} + {\text{3x}}} \right) \frac{{\partial u}}{{\partial x}}} \right] + 5x = 0$ (2.21)

Multiplying Equ 2.20 by w(x) and integrating within the domain:

 $\displaystyle \int_0^{1} {w(x) \cdot \left[\frac{\partial}{\partial x}\left [ (2+3x)\frac{\partial u}{\partial x} \right ]+5x\right] } dx = 0$ (2.22)

Next we will integrate by parts to obtain,

 $\displaystyle \left. {w(x)(2+3x)\frac{\partial u}{\partial x}} \right|_0^{1} - \int_0^{1} {\frac{\partial w}{\partial x}(2+3x)\frac{\partial u}{\partial x}} dx + \int_0^{1} 5wx dx= 0$ (2.23)

Since $\displaystyle w(x=1) = 0$ we can simplify (2.22) in the following way to obtain the weak form,

 $12w(0) - \int_0^1 {\left[ {\frac{{dw}}{{dx}}\left( {2 + 3x} \right)\frac{{du}}{{dx}}} \right]} dx + \int_0^1 {5xw(x)} dx = 0$ (2.24)

Now we can use the polynomial basis functions,

 $w(x) = {b_i}{(x - 1)^i} \to \frac{{dw}}{{dx}} = i{b_i}{(x - 1)^{i - 1}}$ (2.25)
 $\displaystyle u(x) = {b_j}{(x - 1)^j} \to \frac{{du}}{{dx}}j{b_j}{(x - 1)^{j - 1}}$ (2.26)
 $\displaystyle such \quad that \quad w(1) = 0$ (2.27)

Which allows us to write the weak form (2.25) in terms of our matrices $\displaystyle {\color{red}\bold{K}_{FF}}$, $\displaystyle {\color{red}\bold{d}_{F}}$ and $\displaystyle {\color{red}\bold{F}_{F}}$.

 $\underbrace{\left[ \int_0^1 {{\frac{{\partial b_i(x)}}{{\partial x}} (2 + 3x)\frac{{\partial b_j(x)}}{{\partial x}}} } dx \right]}_{{\color{red}\bold{K}_{FF}}} \underbrace{d_i}_{{\color{red}\bold{d}_{F}}} = \underbrace{\left[ \int_0^1 {5xb_i(x)} dx + (12+18x)b_i(x)|_{x=0} \right]}_{{\color{red}\bold{F}_{F}}}$ (2.28)

Next we wish to obtain the solution to $\displaystyle \bold{d_F}$ by solving the following system of equations,

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}} \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$ (2.29)

Solving for the approximating solution, the d matrix can transposed and multiplied by the basis function matrix:

 $U_p^h = {d_0}{b_0} + {d_1}{b_1} + {d_2}{b_2}$ (2.30)

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=2

 $U_p^h(x) = 4 - 2.558(x - 1) + 1.293{(x - 1)^2}$ (2.31)

#### n=4

 $U_p^h(x) = 4 - 2.876(x - 1) + .602{(x - 1)^2} + .310{(x - 1)^3} + .698{(x - 1)^4}$ (2.32)

#### n=6

 $U_p^h(x) = 4 - 2.899(x - 1) + .400{(x - 1)^2} - .121{(x - 1)^3} + .672{(x - 1)^4} + .604{(x - 1)^5} + .378{(x - 1)^6}$ (2.33)

#### n=8

We found that for $\displaystyle n=8$ the error was of the order $\displaystyle 10^{-6}$. Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\left( \begin{matrix} 3.50 & -3.00 & 2.75 & -2.60 & 2.50 & -2.49 & 2.38 & -2.33 \\ -3.00 & 3.67 & -3.90 & 4.00 & -4.05 & 4.07 & -4.08 & 4.09 \\ 2.75 & -3.90 & 4.50 & -4.86 & 5.09 & -5.25 & 5.37 & -5.45 \\ -2.60 & 4.00 & -4.86 & 5.43 & -5.83 & 6.13 & -6.36 & 6.54 \\ 2.50 & -4.05 & 5.09 & -5.83 & 6.39 & -6.82 & 7.16 & -7.44 \\ -2.43 & 4.07 & -5.25 & 6.13 & -6.82 & 7.36 & -7.81 & 8.18 \\ 2.38 & -4.08 & 5.37 & -6.36 & 7.16 & -7.81 & 8.35 & -8.80 \\ -2.33 & 4.09 & -5.45 & 6.54 & -7.44 & 8.18 & -8.80 & 9.33 \\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} -12.83 & 12.42 & -12.15 & 12.17 & -12.12 & 12.09 & -12.07 & 12.06 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.000000000000 \\ & -2.899914248271 \\ & 0.373061337131 \\ & -0.279875467206 \\ & 0.324030513018 \\ & 0.438725044951 \\ & 0.810687654673 \\ & 0.613125135156 \\ & 0.229921935886 \\ \end{matrix} \right)}}$

 $U_p^h(x) = 4 - 2.899(x - 1) + .400{(x - 1)^2} - .121{(x - 1)^3} + .672{(x - 1)^4} + .604{(x - 1)^5} + .378{(x - 1)^6}$

(2.34)

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

#### Plot: Showing the error $e(n) = |u(x)-u_h(x)|$

For the polynomial approximate solution, with the increasing number of nodes, the error decreases. For this particular instance, the error XXXXXXX

### 2. Approximating Soultion Using Fourier Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the fourier solution is obtained the same way. W(x) is now equation 5.35 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=2

 $U_f^h(x) = 4 - 4.125(\cos (x - 1) - 1) - 2.297(\sin (x - 1))$ (2.34)

#### n=4

 $U_f^h(x) = 4 - 11.309(\cos (x - 1) - 1) - 2.066(\sin (x - 1))- 2.422(\cos (2x - 2) - 1) - 0.394(\sin (2x - 2))$ (2.35)

etc...

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

### 3. Approximating Soultion Using Exponential Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the exponential solution is obtained the same way. W(x) is now equation 5.41 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=2

 $U_e^h(x) = 4 - 14.678({\operatorname{e} ^{(x - 1)}} - 1) + 6.434({e^{(2x - 2)}} - 1)$ (2.36)

#### n=4

 $U_e^h(x) = 4 - 58.869({\operatorname{e} ^{(x - 1)}} - 1) + 93.9({e^{(2x - 2)}} - 1) - 72.659({e^{(3x - 3)}} - 1) + 21.596({e^{(4x - 4)}} - 1)$ (2.37)

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\left( \begin{matrix} 1.72 & 2.63 & 3.17 & 3.51 & 3.75 & 3.92 & 4.05 & 4.15 \\ 2.63 & 4.23 & 5.27 & 6.00 & 6.53 & 6.94 & 7.26 & 7.52 \\ 3.17 & 5.27 & 6.74 & 7.83 & 8.67 & 9.33 & 9.87 & 10.31 \\ 3.51 & 6.00 & 7.83 & 9.25 & 10.37 & 11.28 & 12.03 & 12.67 \\ 3.75 & 6.53 & 8.67 & 10.37 & 11.75 & 12.89 & 13.85 & 14.67 \\ 3.92 & 6.94 & 9.33 & 11.28 & 12.89 & 14.25 & 15.41 & 16.41 \\ 4.05 & 7.26 & 9.87 & 12.03 & 13.85 & 15.41 & 16.75 & 17.92 \\ 4.15 & 7.52 & 10.31 & 12.67 & 14.67 & 16.41 & 17.92 & 19.25 \\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} -8.25 & -11.46 & -12.76 & -13.34 & -13.62 & -13.78 & -13.88 & -13.95 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.000000000000 \\ & -416.206015822685 \\ & 1802.371470713180 \\ & -4682.097146605770 \\ & 7738.463302468000 \\ & -8239.102282958180 \\ & 5487.180599937650 \\ & -2083.316654302480 \\ & 344.581494129590 \\ \end{matrix} \right)}}$

The plots for solutions and error at x=.5 are as follows:

# Problem 5.3: Torsion of a Bar

Modified version of HW 4.4 defined on page 2 of Lecture 21 from Fish & Belytschko Page 73 Problem 3.7

## Given

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (5.3.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (5.3.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (5.3.3)
 ${\Gamma _g} = \left\{ 0 \right\} \to u(x = 0) = 4$ (5.3.4)

 ${\Gamma _g} = \left\{ 1 \right\} \to \frac{{du}}{{dx}}\left( {x = 1} \right) = - 6$ (5.3.5)

Where the Weak Form for a bar in torsion can be defined as

 $\int\limits_{0}^{1} \frac{\partial w}{\partial x} a_2(x) \frac{\partial \phi}{{\partial x}}dx = w(1) h + \int\limits_{0}^{1} w(x)f(x)dx ~~~~~~~\forall w ~ with ~ w(0)=0 ~ and ~ u(x=0) = 4$ (5.3.2)

Make use of the Lagrange Interpolation Basis Function (LIBF) with uniform discretization (equidistant nodes) using $\displaystyle m = 4,6,8$ where the LIBF can defined as follows

 $L_{i,m}(x) = \prod_{\begin{smallmatrix}{k=1}\\ k\neq i\end{smallmatrix}}^{m} \frac{x-x_k}{x_i-x_k} = \frac{(x-x_1)}{(x_i-x_1)} \cdots \frac{(x-x_{k-1})}{(x_i-x_{k-1})} \frac{(x-x_{k})}{(x_i-x_{k})}$. (5.3.3)

## Find

a) Explain how Langrange Interpolation Basis Functions (LIBF) are used as Constraint Breaking Solutions
b) Plot all LIBF used
c) Use matlab quad function or WolframAlpha.com to integrate
d) Plot approximate vs exact solution and convergence error versus equidistant nodes (m)

## Solution

### (a)

Langrange Interpolation Basis Functions satisfy the condition of the constraint breaking solution by assuming a value of zero at every node except for the node where the LIBF is defined there it assumes a value of 1. It acts like the Kronecker Delta

 $L_{i,m}(x) = \begin{cases} 1, & \text{if } x = x_i \\ 0, & \text{if } x = x_k \cdots x_m \end{cases}$ (5.3.4)

### (b)

The following is the LIBF for the $\displaystyle m = 4$

 Figure 1: LIBF m = 4([1])

The following is the LIBF for the $\displaystyle m = 6$

 Figure 1: LIBF m = 6([2])

The following is the LIBF for the $\displaystyle m = 8$

 Figure 1: LIBF m = 8([3])

### (c)

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

# Problem 5.4: Torsion of a Bar

Modified version of HW 4.4 defined on page 2 of Lecture 21 from Fish & Belytschko Page 73 Problem 3.7

## Given

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (5.3.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (5.3.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (5.3.3)
 ${\Gamma _g} = \left\{ 0 \right\} \to u(x = 0) = 4$ (5.3.4)

 ${\Gamma _g} = \left\{ 1 \right\} \to \frac{{du}}{{dx}}\left( {x = 1} \right) = - 6$ (5.3.5)

Where the Weak Form for a bar in torsion can be defined as

 $\int\limits_{0}^{1} \frac{\partial w}{\partial x} a_2(x) \frac{\partial \phi}{{\partial x}}dx = w(1) h + \int\limits_{0}^{1} w(x)f(x)dx ~~~~~~~\forall w ~ with ~ w(0)=0 ~ and ~ u(x=0) = 4$ (5.3.6)

Make use of the Lagrange Interpolation Basis Function (LIBF) with uniform discretization (equidistant nodes) using $\displaystyle m = 4,6,8$ where the LIBF can defined as follows

 $L_{i,m}(x) = \prod_{\begin{smallmatrix}{k=1}\\ k\neq i\end{smallmatrix}}^{m} \frac{x-x_k}{x_i-x_k} = \frac{(x-x_1)}{(x_i-x_1)} \cdots \frac{(x-x_{k-1})}{(x_i-x_{k-1})} \frac{(x-x_{k})}{(x_i-x_{k})}$. (5.3.6)

## Find

a) Explain how Langrange Interpolation Basis Functions (LIBF) are used as Constraint Breaking Solutions
b) Plot all LIBF used
c) Use matlab quad function or WolframAlpha.com to integrate
d) Plot approximate vs exact solution and convergence error versus equidistant nodes (m)

## Solution

### (a)

Langrange Interpolation Basis Functions satisfy the condition of the constraint breaking solution by assuming a value of zero at every node except for the node where the LIBF is defined there it assumes a value of 1. It acts like the Kronecker Delta

 $L_{i,m}(x) = \begin{cases} 1, & \text{if } x = x_i \\ 0, & \text{if } x = x_k \cdots x_m \end{cases}$ (5.4.4)

### (b)

The following is the LIBF for the $\displaystyle m = 4$

 Figure 1: LIBF m = 4([4])

The following is the LIBF for the $\displaystyle m = 6$

 Figure 1: LIBF m = 6([5])

The following is the LIBF for the $\displaystyle m = 8$

 Figure 1: LIBF m = 8([6])

### (c)

Our solution $u \left ( x \right )$ and our weight functions can be approximated using the following definitions,

 \displaystyle \begin{align} & u^{h}(x)=\sum_{j} d_{j} b_{j} \\ &\\ & w^{h}(x)=\sum_{i} c_{i} b_{i} \\ \end{align} $(5.3.3)$

To approximate our solution u, we divide our x domain into m nodal subdivisions. Additionally, we will use m LIBF to approximate our solution. These LIBF will have the following form,

 \displaystyle \begin{align} & L_{i,m}(x)=\prod_{k=1}^{m}\frac{(x-x_{k})}{(x_{i}-x_{k})} \\ \end{align} $(5.4.5)$

Therefore, for m=4, we take the first 4 Lagrange Interpolation Basis Functions, or more specifically,

 \displaystyle \begin{align} & L_{1,4}=\frac{(x-x_{2})(x-x_{3})(x-x_{4})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})} \\ \end{align} $(5.4.6)$
 \displaystyle \begin{align} & L_{2,4}=\frac{(x-x_{1})(x-x_{3})(x-x_{4})}{(x_{2}-x_{1})(x_{2}-x_{3})(x_{2}-x_{4})} \\ \end{align} $(5.4.7)$
 \displaystyle \begin{align} & L_{3,4}=\frac{(x-x_{1})(x-x_{2})(x-x_{4})}{(x_{3}-x_{1})(x_{3}-x_{2})(x_{3}-x_{4})} \\ \end{align} $(5.4.8)$
 \displaystyle \begin{align} & L_{4,4}=\frac{(x-x_{1})(x-x_{2})(x-x_{3})}{(x_{4}-x_{1})(x_{4}-x_{2})(x_{4}-x_{3})}\\ \end{align} $(5.4.9)$

Here we want to develop equations for $\tilde{K_{KK}}$ and $\tilde{F_K}$ and solve for $\tilde{d_{K}}$, similar to the previous problems.

### (d)

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

# Problem 5.5 : Calculix-Using CCX module

## Disk Problem

After building disc in cgx module then we might want to know that properties of nodes,lines, surfaces, and bodies so now lets look at how we can get more information what we have done for our model.

First we have to open our program. To do this run the following ;

 $\displaystyle [All Programs]\Rightarrow [bConverged]\Rightarrow [Calculix]\Rightarrow [Calculix Command]$

Open calculix and go to directory where the basic examples stays in your harddrive.

 $\displaystyle ...\Rightarrow [Program Files]\Rightarrow [bConverged] \Rightarrow [Calculix] \Rightarrow [Cgx] \Rightarrow [Examples] \Rightarrow [Basic]$

Then type to open disc example:

 $\displaystyle [cgx]\Rightarrow [disc]$

Now we want to learn about point properties. To see point properties in our graph secren type

 $\displaystyle \begin{matrix} plot & pa & all \end{matrix}$

We need to mesh disc to see the mesh properties to mesh the disc please type following steps

 $\displaystyle \begin{matrix} elty & all & qu8 \end{matrix}$
 $\displaystyle \begin{matrix} mesh & all \end{matrix}$

What we have done so far is that we opened disc and meshed it with quadratic elements. So we are ready to see element properties on our graph screen. To see properties of elements type,

 $\displaystyle \begin{matrix} plot & ma & all \end{matrix}$

Also we can get more specific properties form the special editor of calculix called 'SciTE.' We can go to examples directory and just basically right click on disc example and open with by editing. SciTE automatically will be opened.

From below we can easly see the properties of all components of disc. Especially point properties are boxed.

## Generating Different Meshes

We will do the same procedures as we have done before. Our model is prepared. So we will mesh our model and take properties of different kinds of meshes for 2-D. Our meshes for 2-D are:

 $\displaystyle \begin{matrix} tr3 & 3 node & 2D element \\ gu4 & 4 node & 2D element\\ qu8 & 8 node & 2D element \end{matrix}$ (N)
 $\displaystyle \begin{matrix} Mesh & with & tr3 \end{matrix}$
 $\displaystyle \begin{matrix} Mesh & with & qu4 \end{matrix}$
 $\displaystyle \begin{matrix} Mesh & with & qu8 \end{matrix}$

## CCX module

In this module we will analyze for unknowns that we are interested in. This unknowns can be displacements or temperatures. Also we can find fluxes between these nodes.

For the first step go to directory where beam example stays in your harddrive.

 $\displaystyle ....[Program Files]\Rightarrow [bConverged]\Rightarrow [Calculix]\Rightarrow [ccx]\Rightarrow [test]$

Find beampl.inp and right click on it. Then edit it by selecting open with. SciTe automatically will be opened. Every command inside is ready so we have assumed the we have completed by the mesh and ready for solving.So click tools and click 'solve'.

You have to see that $\displaystyle Job Finished$ Then click tools and click 'post prosess'.

Results

 $\displaystyle Result 1$
 $\displaystyle Result2$

# Problem 5.6: The Quadratic Lagrangian Element Basis Function

Lecture 30-6

## Given:

Given the Linear Lagrangian Element Basis Function,

from Lecture 30-6

## Find:

The Quadratic Lagrangian Element Basis Function

## Solution:

The above is the depiction of the basis function satisfying the CBS for element node i-1

The above is the depiction of the basis function satisfying the CBS for element node i

The above is the depiction of the basis function satisfying the CBS for element node i+1

Next we overlay the three functions to show the shape function for the quadratic Lagrangian

# Problem 5.7 Weak form solution using LLEIF

## Given

Solve the following discrete weak form with Linear Lagrange Element Interpolation Basis.

 $\displaystyle \sum {{c_i}\left[ {\sum {{{\tilde K}_{ij}}{d_j} - {{\tilde F}_i}} } \right]} = 0$ (3.7.1)

Where

 $\displaystyle {{\tilde K}_{ij}} = \int\limits_0^1 {{{b'}_i}(2 + 3x){{b'}_j}dx}$ (3.7.2)
 $\displaystyle {F_i} = {b_i}\left( 1 \right)\left( {12} \right) + \int_0^1 {{b_i}\left( x \right)5x} dx$ (3.7.3)

Such that

 $\displaystyle {u^h}\left( 0 \right) = \sum {{d_j}{b_j}\left( 0 \right) = 4}$ (3.7.4)
 $\displaystyle {w^h}\left( 0 \right) = \sum {{c_i}{b_i}\left( 0 \right) = 0}$ (3.7.5)

## Find

### d. Plot ${u^h}$ vs ${u}$

Plot ${u^h}\left( {0.5} \right) - u\left( {0.5} \right)$ vs n (total no. of dof's)

## Solution

### a. Using LLEBF in constraint breaking solution(CBS)

We have following constraint to be satisfied to use weak form

 $\displaystyle {w^h}\left( 0 \right) = \sum {{c_i}{b_i}\left( 0 \right) = 0}$ (3.7.6)

At node 1 (i.e. at x=0) all the basis function are zero except 1.

 $\displaystyle {w^h} = {c_1}{b_1} + {c_2}{b_2} + {c_3}{b_3} + {c_4}{b_4} + \cdots$ (3.7.7)

 $\displaystyle {w^h} = {c_1}(1) + {c_2}(0) + {c_3}(0) + {c_4}(0) + \cdots$ (3.7.8)
 $\displaystyle {w^h} = {c_1}$ (3.7.9)

For

 $\displaystyle {w^h}\left( 0 \right) = 0$ (3.7.10)
 $\displaystyle {c_1} = 0$ (3.7.11)

To satisfy homogeneous boundary condition while using Linear Lagrange Element Interpolation Functions, we just have to select our $w^h$ such that ${c_1} = 0$.

### c Formulation of K and F on element level

We can divide the total domain into elemental domain.

 $\displaystyle \Omega = \bigcup\limits_{e = 1}^{nel} {{w^e}}$ (3.7.12)

Values of K and F will be summation of K and F over the elements

 $\displaystyle {{\tilde K}_{ij}} = \sum\limits_{e = 1}^{nel} {K_{ij}^e}$ (3.7.13)

 $\displaystyle K_{ij}^e = \int\limits_{{w_e}} {{{b'}^e}_i(2 + 3x){{b'}^e}_jdx}$ (3.7.14)

 $\displaystyle {{\tilde F}_i} = \sum\limits_{e = 1}^{nel} {F_i^e}$ (3.7.15)
 $\displaystyle {F_i}^e = \int_{{w_e}} {b_i^e5xdx}$ (3.7.16)

There will be two basis function for each element. Generalised values of these basis function will be as follow

 $\displaystyle {b_e} = \frac{{x - {x_{e + 1}}}}{{{x_e} - {x_{e + 1}}}}$ (3.7.17)
 $\displaystyle {b_{e + 1}} = \frac{{x - {x_e}}}{{{x_{e + 1}} - {x_e}}}$ (3.7.18)

### d. Solution using Matlab

Plot of exact solution and approximate solution with 2 elements

Figure 1: For elements = 2 [8]

Plot of exact solution and approximate solution with 4 elements

Figure 1: For elements = 4 [9]

Plot of exact solution and approximate solution with 6 elements

Figure 1: For elements = 6 [10]

Plot of exact solution and approximate solution with 8 elements

Figure 1: For elements = 8 [11]

Plot of error vs dof's

 $\displaystyle err = u\left( {0.5} \right) - {u^h}\left( {0.5} \right)$ (3.7.19)
 $\displaystyle dof's = nel + 1$ (3.7.20)

where nel= number of elements.

Figure 1: Error vs dof's [12]

for 16 elements

dof's=16+1=17

Figure 1: Error vs dof's [13]

## comment

Solution is converging with increasing no. of elements. For n=16 error is less than 0.0001.

# Problem 5.8:Using LLEBF to solve G1DM1.0/D1

## Given

G1DM1.0/D1 : Model for elastic bar problem was given in meeting 9.

 $\displaystyle \frac{\mathrm{d} }{\mathrm{d} x}\left [ (2+3x)\frac{\mathrm{d} u}{\mathrm{d} x} \right ]+5x=0$ (8.1)
 $\displaystyle u(1)=4$ (8.2)
 $\displaystyle \frac{\mathrm{d} u(x=0)}{\mathrm{d} x}=-3$ (8.3)

Weak Form

 $\displaystyle \int_{\Omega }\frac{\mathrm{d} w}{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx-\left ( wAt \right )_{\Gamma _{h}}-\int_{\Omega }wbdx=0 \Rightarrow \Omega =]0,1[$ (8.4)

We have to manipulate weak form in order to use for element matrices. As to element mesh we use weak form for each elements and then we can take summation for all elements to get global matrices.

 $\displaystyle \sum\limits_{e=1}^{nel}{{{w}^{eT}}\left\{ \int\limits_{{{x}_{1}}^{e}}^{{{x}_{2}}^{e}}{{{B}^{eT}}{{A}^{e}}{{E}^{e}}{{B}^{e}}dx{{d}^{e}}-\int\limits_{{{x}_{1}}^{e}}^{x_{2}^{e}}{{{N}^{eT}}bdx-{{({{N}^{eT}}{{A}^{e}}\overline{t})}_{x=0}}}} \right\}=0}$ (8.5)

Approximate functions gives as in F&B Book,

 $\displaystyle u^{h}=\sum_{e=1}^{nel}N^{e}d^{e}=\left ( \sum_{e=1}^{nel}N^{e}L^{e}\right )\mathbf{d}$ (8.6)
 $\displaystyle w^{h}=\sum_{e=1}^{nel}N^{e}w^{e}=\left ( \sum_{e=1}^{nel}N^{e}L^{e}\right )\mathbf{w}$ (8.7)

Where N is called as 'shape functions' and L is called 'gather matrices'.

 $\displaystyle N^{e}=\frac{1}{l^{e}}[\begin{matrix} x_{2}^{e}-x & x-x_{1}^{e} \end{matrix}]$ (8.8)
 $\displaystyle B^{e}=\frac{1}{l^{e}}[\begin{matrix} -1 & 1 \end{matrix}]$ (8.9)

## Solution

### Analytical Solution

 $\displaystyle \int \left ( (2+3x)\frac{\mathrm{d} u}{\mathrm{d} x} \right )dx=-\int (5x)dx+c_{1}$ (8.10)
 $\displaystyle 2\frac{\mathrm{d} u}{\mathrm{d} x}(x=0)=0+c_{1} \Rightarrow c_{1}=-6$ (8.11)
 $\displaystyle u=-\frac{5}{2}\int \frac{x^2}{2+3x}dx-\int \frac{6}{2+3x}dx+c_{2}$ (8.12)
 $\displaystyle -\frac{5}{2}\left [ \frac{1}{54}(9x^2-12x+8log(3x+2)-12) \right ]-\left [ \frac{1}{3}log(3x+2) \right ]+c_{2}$ (8.13)

Using other boundry condition ;

 $\displaystyle c_{2}=\frac{119}{36}+\frac{64}{27}log(5)$ (8.14)
 $\displaystyle u(x)=-\frac{5}{12}x^2+\frac{5}{9}x-\frac{64}{27}log(2+3x)+\frac{139}{36}+\frac{64}{27}log(5)$ (8.15)

### FEA solutions

When evaluating unknowns on distributed nodes our shape functions must meet some conditions.Our approximate functions must be complete in each element. As saying completeness we are saying that by using our nodal values we must be able to identify coefficients in our approximate solutions. After some manipulation we conclude with shape functions. For each nodes we have shape functions and for middle nodes shape functions come from adjacent elements.

For example two element mesh ;

 $\displaystyle N=N^{(1)}L^{(1)}+N^{(2)}L^{(2)}=[\begin{matrix} N_{1}^{(1)} & N_{2}^{(1)}+N_{1}^{(2)} & N_{2}^{(2)} \end{matrix}]$ (8.16)

For this example at x1 (first node) only N1 will be 1. Others will be zero. So we verify CBS.[1]

#### Using 4 elements

 $\displaystyle \Omega =]0,1[$

Lets find shape functions for 4 elements.

 \displaystyle\begin{align} &N^{(1)}=\frac{1}{0.25}[\begin{matrix} 0.25-x & x \end{matrix}]\\ &N^{(2)}=\frac{1}{0.25}[\begin{matrix} 0.5-x & x-0.25 \end{matrix}]\\ &N^{(3)}=\frac{1}{0.25}[\begin{matrix} 0.75-x & x-0.5 \end{matrix}]\\ &N^{(4)}=\frac{1}{0.25}[\begin{matrix} 1-x & x-0.75 \end{matrix}] \end{align} (8.17)

These shape functions are useful for us to make connection between nodes. After defining these tools now we can find stiffness matrices for all elements. Then after finding all of stiffness matrices for elements we can construct global stiffness matrix.

 $\displaystyle K^{e}=\int_{x_{1}^{e}}^{x_{2}^{e}}B^{eT}A^{e}E^{e}B^{e}dx$ (8.18)
 \displaystyle\begin{align} {{K}^{e}}&=\int\limits_{x_{1}^{e}}^{x_{2}^{e}}{\frac{1}{{{l}^{e}}}\left[ \begin{matrix} -1 \\ 1 \\ \end{matrix} \right]}(2+3x)\frac{1}{{{l}^{e}}}[\begin{matrix} -1 & 1 \\ \end{matrix}]dx \\ &=\frac{1}{{{({{l}^{e}})}^{2}}}\left[ \begin{matrix} 1 & -1 \\ -1 & 1 \\ \end{matrix} \right]\left[ 2(x_{2}^{e}-x_{1}^{e})+\frac{3}{2}(x_{2}^{e2}-x_{1}^{e2}) \right] \end{align} (8.19)

 $\displaystyle {{K}^{e}}=\frac{1}{{{l}^{e}}}\left[ \begin{matrix} 1 & -1 \\ -1 & 1 \\ \end{matrix} \right]\left[ 2+\frac{3}{2}(x_{2}^{e}+x_{1}^{e}) \right]$ (8.20)
 $\displaystyle K^{(1)}=\begin{bmatrix} 9.5 & -9.5\\ -9.5& 9.5 \end{bmatrix}$ (8.21)
 $\displaystyle K^{(2)}=\begin{bmatrix} 12.5 & -12.5\\ -12.5& 12.5 \end{bmatrix}$ (8.22)
 $\displaystyle K^{(3)}=\begin{bmatrix} 15.5 & -15.5\\ -15.5& 15.5 \end{bmatrix}$ (8.23)
 $\displaystyle K^{(4)}=\begin{bmatrix} 18.5 & -18.5\\ -18.5& 18.5 \end{bmatrix}$ (8.24)

We can construct global stiffness matrix now. This will be distribution of all four element stiffness matrices into appropriate places. For our convenience i will use direct method. But it is useful using below equation for computing purposes.

 $\displaystyle K=\sum_{e=1}^{nel}L^{eT}K^{e}L^{e}=L^{1T}K^{1}L^{1}+L^{2T}K^{2}L^{2}+L^{3T}K^{3}L^{3}+L^{4T}K^{4}L^{4}$ (8.25)
 $\displaystyle K=\left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} 9.5 \\ -9.5 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} & \begin{matrix} \begin{matrix} -9.5 \\ 22 \\ \end{matrix} \\ -12.5 \\ 0 \\ 0 \\ \end{matrix} \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ -12.5 \\ \end{matrix} \\ 28 \\ \begin{matrix} -15.5 \\ 0 \\ \end{matrix} \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ 0 \\ \end{matrix} \\ -15.5 \\ 34 \\ -18.5 \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ 0 \\ \end{matrix} \\ 0 \\ -18.5 \\ 18.5 \\ \end{matrix} \\ \end{matrix} \right]$ (8.26)

The next step we have to find force matrices. Which ha two components. Lets first find force matrix at natural boundary condition and second distributed force matrix on our domain.

 $\displaystyle f_{\Gamma _{h}}=-(N^{eT}A^{e}b)_{\Gamma _{h}^{e}}=-N^{eT}(x_{1})(2+3x)_{x=0}(-3)=6N^{eT}(x_{1})$ (8.27)
 $\displaystyle f_{{\Gamma }^{(1)}}=6\left[ \begin{matrix} N_{1}^{(1)}({{x}_{1}}) \\ N_{2}^{(1)}({{x}_{2}}) \\ \end{matrix} \right]=\left[ \begin{matrix} 6 \\ 0 \\ \end{matrix} \right]$ (8.28)
 \displaystyle\begin{align} &f_{{\Gamma }^{(2)}}=f_{{\Gamma }^{(3)}}=f_{{\Gamma }^{(4)}}=6\left[ \begin{matrix} N_{1}^{(e)}({{x}_{1}}) \\ N_{2}^{(e)}({{x}_{2}}) \\ \end{matrix} \right]=\left[ \begin{matrix} 0 \\ 0 \\ \end{matrix} \right] &for & e=2,3,4 \end{align} (8.29)
 \displaystyle\begin{align} f_{\Gamma _{h}} &=\sum_{e=1}^{nel}L^{eT}f_{\Gamma _{h}}^{e} \\ &=\left[ \begin{matrix} \begin{matrix} 6 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right] \end{align} (N)
 \displaystyle\begin{align} f_{\Omega }^{e} & =\int_{x_{1}^{e}}^{x_{2}^{e}}N^{eT}b(x)dx\\ &=\int\limits_{x_{1}^{e}}^{x_{2}^{e}}{\frac{1}{{{l}^{e}}}\left[ \begin{matrix} x_{2}^{e}-x \\ x-x_{1}^{e} \\ \end{matrix} \right]}5xdx \\ &=\frac{5}{{{l}^{e}}}\left[ \begin{matrix} \frac{x_{2}^{e}(x_{2}^{e2}-x_{1}^{e2})}{2}-\frac{(x_{2}^{e3}-x_{1}^{e3})}{3} \\ \frac{(x_{2}^{e3}-x_{1}^{e3})}{3}-\frac{x_{1}^{e}(x_{2}^{e2}-x_{1}^{e2})}{2} \\ \end{matrix} \right] \end{align} (8.30)
 $\displaystyle f_{\Omega }^{(1)}=\left[ \begin{matrix} 0.05208334 \\ 0.10416666 \\ \end{matrix} \right]$ (8.31)
 $\displaystyle f_{\Omega }^{(2)}=\left[ \begin{matrix} 0.2083334 \\ 0.2604166 \\ \end{matrix} \right]$ (8.32)
 $\displaystyle f_{\Omega }^{(3)}=\left[ \begin{matrix} 0.36458334 \\ 0.41666666 \\ \end{matrix} \right]$ (8.33)
 $\displaystyle f_{\Omega }^{(4)}=\left[ \begin{matrix} 0.52083334 \\ 0.57291666 \\ \end{matrix} \right]$ (8.34)
 \displaystyle\begin{align} f_{\Omega }&=\sum_{e=1}^{4}L^{eT}f_{\Omega }^{e} &\left[ \begin{matrix} \begin{matrix} 0.05208334 \\ 0.3125 \\ \end{matrix} \\ 0.625 \\ 0.9375 \\ 0.57291666 \\ \end{matrix} \right] \end{align} (8.35)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} 9.5 \\ -9.5 \\ \end{matrix} & \begin{matrix} -9.5 \\ 22 \\ \end{matrix} & \begin{matrix} 0 \\ -12.5 \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 & 0 \\ \end{matrix} \\ \begin{matrix} 0 & 0 \\ \end{matrix} \\ \end{matrix} \\ 0 & -12.5 & 28 & \begin{matrix} -15.5 & 0 \\ \end{matrix} \\ 0 & 0 & -15.5 & \begin{matrix} 34 & -18.5 \\ \end{matrix} \\ 0 & 0 & 0 & \begin{matrix} -18.5 & 18.5 \\ \end{matrix} \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ \end{matrix} \\ {{d}_{3}} \\ {{d}_{4}} \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 0.05208334 \\ 0.3125 \\ \end{matrix} \\ 0.625 \\ 0.9375 \\ 0.57291666 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 6 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ {{r}_{1}} \\ \end{matrix} \right]$ (8.36)
 $\displaystyle \left[ \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ {{d}_{4}} \\ \end{matrix} \right]=\left[ \begin{matrix} 2.0257 \\ 1.3886 \\ 0.8794 \\ 0.4285 \\ \end{matrix} \right]$ (8.37)

Comparing results with exact ones at the nodes;

 \displaystyle\begin{align} &\left[ \begin{matrix} 2.0257 \\ 1.3886 \\ 0.8794 \\ 0.4285 \\ \end{matrix} \right]+\left[ \begin{matrix} 4 \\ 4 \\ 4 \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} 6.0257 \\ 5.3886 \\ 4.8794 \\ 4.4285 \\ \end{matrix} \right] &{{d}_{exact}}=\left[ \begin{matrix} 6.0331 \\ 5.3911 \\ 4.8802 \\ 4.4286 \\ \end{matrix} \right] \end{align} (8.38)
 $\displaystyle 4 elements$

#### Using 6 elements

Shape functions for 6 elements

 \displaystyle\begin{align} & N^{(1)}=\frac{1}{0.1666}[\begin{matrix} 0.1666-x & x \end{matrix}]\\ &N^{(2)}=\frac{1}{0.1666}[\begin{matrix} 0.3333-x & x-0.16666 \end{matrix}]\\ & N^{(3)}=\frac{1}{0.1666}[\begin{matrix} 0.5-x & x-0.3333 \end{matrix}]\\ &N^{(4)}=\frac{1}{0.1666}[\begin{matrix} 0.6666-x & x-0.5 \end{matrix}]\\ & N^{(5)}=\frac{1}{0.1666}[\begin{matrix} 0.8332-x & x-0.6666 \end{matrix}]\\ & N^{(6)}=\frac{1}{0.1666}[\begin{matrix} 1-x & x-0.8332 \end{matrix}]\\ \end{align} (8.39)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} 13.5048 & -13.5048 & 0 & 0 \\ \end{matrix} & 0 & 0 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} -13.5048 & 30.0105 & -16.5057 & 0 \\ \end{matrix} & 0 & 0 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & -16.5057 & 36.0132 & -19.5075 \\ \end{matrix} & 0 & 0 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & -19.5075 & 42.0157 \\ \end{matrix} & -22.5084 & 0 & 0 \\ \end{matrix} \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & 0 & -22.5084 \\ \end{matrix} & 48.0164 & -25.508 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & 0 & 0 \\ \end{matrix} & -25.508 & 54.008 & -28.5 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & 0 & 0 \\ \end{matrix} & 0 & -28.5 & 28.5 \\ \end{matrix} \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ {{d}_{4}} \\ \end{matrix} \\ {{d}_{5}} \\ {{d}_{6}} \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 6.0232199942 \\ 0.13890282 \\ 0.277972087 \\ 0.416641681 \\ \end{matrix} \\ 0.555277743 \\ 0.694970317 \\ 0.394287575 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ {{r}_{1}} \\ \end{matrix} \right]$ (8.40)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ \end{matrix} \\ {{d}_{4}} \\ {{d}_{5}} \\ {{d}_{6}} \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 2.0292 \\ 1.5831 \\ 1.2098 \\ \end{matrix} \\ 0.8797 \\ 0.575 \\ 0.2845 \\ \end{matrix} \right]$ (8.41)
 \displaystyle\begin{align} & \left[ \begin{matrix} \begin{matrix} 2.0292 \\ 1.5831 \\ 1.2098 \\ \end{matrix} \\ 0.8797 \\ 0.575 \\ 0.2845 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 4 \\ 4 \\ 4 \\ \end{matrix} \\ 4 \\ 4 \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 6.0292 \\ 5.5831 \\ 5.2098 \\ \end{matrix} \\ 4.8797 \\ 4.575 \\ 4.2845 \\ \end{matrix} \right] & {{d}_{exact}}=\left[ \begin{matrix} \begin{matrix} 6.0331 \\ 5.5853 \\ 5.2109 \\ \end{matrix} \\ 4.8802 \\ 4.5753 \\ 4.2847 \\ \end{matrix} \right] \end{align} (8.42)
 $\displaystyle 6 elements$

#### Using 8 elements

Shape functions for 8 elements

 \displaystyle\begin{align} & N^{(1)}=\frac{1}{0.125}[\begin{matrix} 0.125-x & x \end{matrix}]\\ &N^{(2)}=\frac{1}{0.125}[\begin{matrix} 0.25-x & x-0.125 \end{matrix}]\\ & N^{(3)}=\frac{1}{0.125}[\begin{matrix} 0.375-x & x-0.25 \end{matrix}]\\ &N^{(4)}=\frac{1}{0.125}[\begin{matrix} 0.5-x & x-0.375 \end{matrix}]\\ & N^{(5)}=\frac{1}{0.125}[\begin{matrix} 0.625-x & x-0.5 \end{matrix}]\\ & N^{(6)}=\frac{1}{0.125}[\begin{matrix} 0.75-x & x-0.625 \end{matrix}]\\ & N^{(7)}=\frac{1}{0.125}[\begin{matrix} 0.875-x & x-0.75 \end{matrix}]\\ & N^{(8)}=\frac{1}{0.125}[\begin{matrix} 1-x & x-0.875 \end{matrix}]\\ \end{align} (8.43)
 $\displaystyle \left[ \begin{matrix} 17.5 & -17.5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -17.5 & 38 & -20.5 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -20.5 & 44 & -23.5 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -23.5 & 50 & -26.5 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -26.5 & 56 & -29.5 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -29.5 & 62 & -32.5 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -32.5 & 68 & -35.5 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -35.5 & 74 & -38.5 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -38.5 & 38.5 \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ \end{matrix} \\ {{d}_{4}} \\ {{d}_{5}} \\ {{d}_{6}} \\ \end{matrix} \\ {{d}_{7}} \\ {{d}_{8}} \\ {{d}_{9}} \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} 6.01302084 \\ 0.078125 \\ \end{matrix} \\ 0.15625 \\ \end{matrix} \\ 0.234375 \\ 0.31250004 \\ 0.390625 \\ \end{matrix} \\ 0.46874373 \\ 0.546875 \\ 0.29947916 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} \begin{matrix} 0 \\ 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ {{r}_{1}} \\ \end{matrix} \right]$ (8.44)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ \end{matrix} \\ {{d}_{3}} \\ \end{matrix} \\ {{d}_{4}} \\ {{d}_{5}} \\ \end{matrix} \\ {{d}_{6}} \\ {{d}_{7}} \\ {{d}_{8}} \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} \begin{matrix} 2.0312 \\ 1.6876 \\ \end{matrix} \\ 1.3904 \\ 1.1246 \\ 0.8800 \\ \end{matrix} \\ 0.6497 \\ 0.4286 \\ 0.213 \\ \end{matrix} \right]$ (8.45)
 \displaystyle\begin{align} &\left[ \begin{matrix} \begin{matrix} \begin{matrix} 2.0312 \\ 1.6876 \\ \end{matrix} \\ 1.3904 \\ 1.1246 \\ 0.8800 \\ \end{matrix} \\ 0.6497 \\ 0.4286 \\ 0.213 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} \begin{matrix} 4 \\ 4 \\ \end{matrix} \\ 4 \\ 4 \\ 4 \\ \end{matrix} \\ 4 \\ 4 \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} \begin{matrix} 6.0312 \\ 5.6876 \\ \end{matrix} \\ 5.3904 \\ 5.1246 \\ 4.88 \\ \end{matrix} \\ 4.6497 \\ 4.4286 \\ 4.213 \\ \end{matrix} \right] &{{d}_{exact}}=\left[ \begin{matrix} \begin{matrix} \begin{matrix} 6.0331 \\ 5.6866 \\ \end{matrix} \\ 5.3911 \\ 5.1249 \\ 4.8802 \\ \end{matrix} \\ 4.6498 \\ 4.4286 \\ 4.213 \\ \end{matrix} \right] \end{align} (8.46)
 $\displaystyle 8 elements$

# References

1. Fish and Belytschko FEA Anlaysis

# Contributing Members & Referenced Lecture

 Team Contribution Table Problem Number Lecture Assigned To Solved By Typed By Proofread By 5.1 Lecture 26-2 Chris,Tim Chris,Tim Chris,Tim XXX 5.2 Lecture 26-2 Chris,Tim Chris,Tim Chris,Tim XXX 5.3 Lecture 29-6 Danny,Phil XXX XXX XXX 5.4 Lecture 29-6 Danny,Phil XXX XXX XXX 5.5 Lecture 29-7 Chris,Tim Oztekin Oztekin XXX 5.6 Lecture 30-6 Danny,Phil XXX XXX xxx 5.7 Lecture 30-6 Erman,Taj Shahtaj Shahtaj Oztekin 5.8 Lecture 30-7 Erman,Taj Oztekin Oztekin XXX