User:Eml5526.s11.team2/hwk5

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

Exact Solution

Exact Solution Showing all Steps Find the exact solution of $\displaystyle u(x)$ As presented in the problem statement we will solve the following ODE $\displaystyle \frac{d}{{dx}}\left[ (2+3x)\frac{{du}}{{dx}}\right] + 5x = 0$ $\displaystyle \begin{align} s.t.& \quad u(x=0) & = 4 \\ and & \quad \frac{{du}}{{dx}}(x=1)&= -6 \end{align}$ (1.6) Rearranging the ODE provided in (1.1) and integrating once gives $\frac{{du}}{{dx}} = \frac{-\frac{5}{2}x^2 + a}{\quad 2+3x}$ (1.7) Next we must satisfy the essential boundary conditions $\frac{{du}}{{dx}}(x=1) = \frac{-\frac{5}{2}x^2 + a}{\quad 2+3x} = -6$ $\implies a = \frac{-55}{2}$ (1.8) Therefore the first order separable ODE has the form $\frac{du}{dx} = \frac{-\frac{5}{2}x^2 - \frac{55}{2}}{\quad 2+3x}$ ${du} = \frac{-\frac{5}{2}x^2 - \frac{55}{2}}{\quad 2+3x} dx$ (1.9) To solve for $u(x)$ in (1.9) we must evaluate the following integral $\int {du} = \int {\frac{-\frac{5}{2}x^2 - \frac{55}{2}}{\quad 2+3x} dx}$ (1.10) Recall that a polynomial P(x) with a numerator N(x), denominator D(x), quotient Q(x), and remainder R(x) such that P(x) has degree(D(x)) < degree(N(x)), then P(x) can be written as $\frac{N(x)}{D(x)} = Q(x) + \frac{R(x)}{D(x)}$ (1.11) Next we will use long division to simplify the integrand. $\begin{matrix} \;\, {\color{blue}-\frac{5}{6}x} \; {\color{blue} + \frac{5}{9}} \\ {\color{cyan}2 + 3x}\overline{) -\frac{5}{2}x^2 + 0 x - \frac{55}{2}}\\ \quad \;\; \underline{\;\; -\frac{5}{2}x^2 - \; \frac{5}{3}x}\\ \qquad\qquad\quad\quad\quad\; \frac{5}{3}x - \frac{55}{2}\\ \qquad\qquad\quad\quad\quad\; \underline{\frac{5}{3}x + \frac{10}{9}}\\ \qquad\qquad\qquad\qquad\qquad {\color{red}-\frac{515}{18}}\\ \end{matrix}$ (1.12) The polynomial in the integrand can therefore be written as $\underbrace{\frac{-\frac{5}{2}x^2 - \frac{55}{2}}{{\color{cyan} 2+3x}}}_{\frac{N(x)}{{\color{cyan}D(x)}}} = {\color{blue}\underbrace{-\frac{5}{6}x+\frac{5}{9}}_{Q(x)}} - \underbrace{ \frac{{\color{red}515}}{{\color{red}18}({\color{cyan}2+3x})}}_{{\color{red}R(x)}/{\color{cyan}D(x)}}$ (1.13) Now the integration becomes trivial since it can be written as a sum of three easily evaluated integrals $\displaystyle \int {\frac{-\frac{5}{2}x^2 - \frac{55}{2}}{2+3x} dx} = -\int{\frac{515}{18}\left(\frac{1}{2+3x} \right)dx} - \int{\frac{5}{6}xdx} + \int{\frac{5}{9}dx}$ (1.14) Then if we integrate using u-substitution for the first integral in (1.14) by letting $\displaystyle u = 2+3x \quad \implies \quad du=3dx$ (1.15) We get the following result $\displaystyle \begin{align} -\int{\frac{515}{18}\left( \frac{1}{2+3x} \right) dx} & = -\int{\frac{515}{18} \left( \frac{1}{3} \right) \left( \frac{1}{u}du \right) } \\ & = -\frac{515}{54} ln|2+3x| + b \end{align}$ (1.16) $\displaystyle u(x) = \int {\frac{-\frac{5}{2}x^2 - \frac{55}{2}}{2+3x} dx} = -\frac{515}{54} ln|2+3x| - \frac{5}{12}x^2 + \frac{5}{9}x + c_{5}$ (1.17) To determine the $\displaystyle c_{5}$, we must apply the necessary boundary conditions to the result obtained in (1.17) $\displaystyle u(x=0) = -\frac{515}{54} ln(5) - \frac{5}{12} + \frac{5}{9} + c_{5} = 4$ (1.18) Solving for $\displaystyle c_{5}$ gives the following algebraic expression $\displaystyle c_{5} = 4 + \frac{515}{54} ln(2)$ (1.19) Now $\displaystyle c_{5}$ in (1.19) can be substituted into (1.17) where $\displaystyle u(x)$ can explicitly be expressed as   $\displaystyle u(x) = -\frac{515}{54} ln|2+3x| - \frac{5}{12}x^2 + \frac{5}{9}x + 4 + \frac{515}{54} ln(2)$ (1.20)

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

MATLAB Code for Problem 5.1

clear all; clc; format long syms x Uh; U = 241/54*log((2+3*x)/2)+(20*x-15*x^2+144)/36; % Exact Solution m = 9; % Global Equidistant Nodes d_0 = 4; % Basis Function Vector b_i = sym(ones(m,1)); for i=1:m if rem(i,2)==0 b_i(i)=sin(i/2*(x-1)); else b_i(i)=cos((floor(i/2)+1)*(x-1))-1; end end % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) - subs((30)*b_i(i),x,0); end d = zeros(m+1,1); Kfff = zeros(m+1,m+1); Kfff(1,1) = 1; Kfff(2:m+1,2:m+1) = Kff % Displaying Kff Fff = d_0*ones(m+1,1) Fff(2:m+1,1) =Ff % Displaying Ff d = Kfff\Fff; % Solve for d i=1,2,...,m Uh = dot(d,b_i) % This come out simplified with fractions. % You can just look at b_i and d to write out Uh U05=vpa(subs(U,0.5)) Uh05=vpa(subs(Uh,0.5)) error = abs(U05 - Uh05)

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

MATLAB Code for Problem 5.1

clear all; clc; format long syms x Uh; U = 241/54*log((2+3*x)/2)+(20*x-15*x^2+144)/36; % Exact Solution m = 9; % Global Equidistant Nodes d_0 = 4; % Basis Function Vector b_i = sym(ones(m,1)); for i=1:m if rem(i,2)==0 b_i(i)=sin(i/2*(x-1)); else b_i(i)=cos((floor(i/2)+1)*(x-1))-1; end end % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) - subs((30)*b_i(i),x,0); end d = zeros(m+1,1); Kfff = zeros(m+1,m+1); Kfff(1,1) = 1; Kfff(2:m+1,2:m+1) = Kff % Displaying Kff Fff = d_0*ones(m+1,1) Fff(2:m+1,1) =Ff % Displaying Ff d = Kfff\Fff; % Solve for d i=1,2,...,m Uh = dot(d,b_i) % This come out simplified with fractions. % You can just look at b_i and d to write out Uh U05=vpa(subs(U,0.5)) Uh05=vpa(subs(Uh,0.5)) error = abs(U05 - Uh05)

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)|$

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

MATLAB Code for Problem 5.1

clear all; clc; format long syms x Uh; U = 241/54*log((2+3*x)/2)+(20*x-15*x^2+144)/36; % Exact Solution m = 9; % Global Equidistant Nodes d_0 = 4; % Basis Function Vector b_i = sym(ones(m,1)); for i=1:m if rem(i,2)==0 b_i(i)=sin(i/2*(x-1)); else b_i(i)=cos((floor(i/2)+1)*(x-1))-1; end end % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) - subs((30)*b_i(i),x,0); end d = zeros(m+1,1); Kfff = zeros(m+1,m+1); Kfff(1,1) = 1; Kfff(2:m+1,2:m+1) = Kff % Displaying Kff Fff = d_0*ones(m+1,1) Fff(2:m+1,1) =Ff % Displaying Ff d = Kfff\Fff; % Solve for d i=1,2,...,m Uh = dot(d,b_i) % This come out simplified with fractions. % You can just look at b_i and d to write out Uh U05=vpa(subs(U,0.5)) Uh05=vpa(subs(Uh,0.5)) error = abs(U05 - Uh05)

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

Exact Solution

Exact Solution Integration Find the exact solution of $\displaystyle u(x)$ As presented in the problem statement we will solve the following ODE $\displaystyle \frac{d}{{dx}}\left[ (2+3x)\frac{{du}}{{dx}}\right] + 5x = 0$ $\displaystyle \begin{align} s.t.& \quad u(x=1) & = 4 \\ and & \quad \frac{{du}}{{dx}}(x=0)&= -6 \end{align}$ (2.6) Rearranging the ODE provided in (2.1) and integrating once gives $\frac{{du}}{{dx}} = \frac{-\frac{5}{2}x^2 + a}{\quad 2+3x}$ (2.7) Next we must satisfy the essential boundary conditions $\frac{{du}}{{dx}}(x=0) = \frac{-\frac{5}{2}x^2 + a}{\quad 2+3x} = -6$ $\implies \frac{a}{2} = -6$ $\implies a = -12$ (2.8) Therefore the first order separable ODE has the form $\frac{du}{dx} = \frac{-\frac{5}{2}x^2 - 12}{\quad 2+3x}$ ${du} = \frac{-\frac{5}{2}x^2 - 12}{\quad 2+3x} dx$ (2.9) To solve for $u(x)$ in (2.9) we must evaluate the following integral $\int {du} = \int {\frac{-\frac{5}{2}x^2 - 12}{\quad 2+3x} dx}$ (2.10) Recall that a polynomial P(x) with a numerator N(x), denominator D(x), quotient Q(x), and remainder R(x) such that P(x) has degree(D(x)) < degree(N(x)), then P(x) can be written as $\frac{N(x)}{D(x)} = Q(x) + \frac{R(x)}{D(x)}$ (2.11) Next we will use long division to simplify the integrand. $\begin{matrix} \;\, {\color{blue}-\frac{5}{6}x} \; {\color{blue} + \frac{5}{9}} \\ {\color{cyan}2 + 3x}\overline{) -\frac{5}{2}x^2 + 0 x - 12}\\ \quad \;\; \underline{\;\; -\frac{5}{2}x^2 - \; \frac{5}{3}x}\\ \qquad\qquad\quad\quad\quad\; \frac{5}{3}x - 12\\ \qquad\qquad\quad\quad\quad\; \underline{\frac{5}{3}x + \frac{10}{9}}\\ \qquad\qquad\qquad\qquad\qquad {\color{red}-\frac{118}{9}}\\ \end{matrix}$ (2.12) The polynomial in the integrand can therefore be written as $\underbrace{\frac{-\frac{5}{2}x^2 - 12}{{\color{cyan} 2+3x}}}_{\frac{N(x)}{{\color{cyan}D(x)}}} = {\color{blue}\underbrace{-\frac{5}{6}x+\frac{5}{9}}_{Q(x)}} - \underbrace{ \frac{{\color{red}118}}{{\color{red}9}({\color{cyan}2+3x})}}_{{\color{red}R(x)}/{\color{cyan}D(x)}}$ (2.13) Now the integration becomes trivial since it can be written as a sum of three easily evaluated integrals $\displaystyle \int {\frac{-\frac{5}{2}x^2 - 12}{2+3x} dx} = -\int{\frac{118}{9}\left(\frac{1}{2+3x} \right)dx} - \int{\frac{5}{6}xdx} + \int{\frac{5}{9}dx}$ (2.14) Then if we integrate using u-substitution for the first integral in (2.14) by letting $\displaystyle u = 2+3x \quad \implies \quad du=3dx$ (2.15) We get the following result $\displaystyle \begin{align} -\int{\frac{118}{9}\left( \frac{1}{2+3x} \right) dx} & = -\int{\frac{118}{9} \left( \frac{1}{3} \right) \left( \frac{1}{u}du \right) } \\ & = -\frac{118}{27} ln|2+3x| + b \end{align}$ (2.16) $\displaystyle u(x) = \int {\frac{-\frac{5}{2}x^2 - 12}{2+3x} dx} = -\frac{118}{27} ln|2+3x| - \frac{5}{12}x^2 + \frac{5}{9}x + c_{5}$ (2.17) To determine the $\displaystyle c_{5}$, we must apply the necessary boundary conditions to the result obtained in (2.22) $\displaystyle u(x=1) = -\frac{118}{27} ln(5) - \frac{5}{12} + \frac{5}{9} + c_{5} = 4$ (2.18) Solving for $\displaystyle c_{5}$ gives the following algebraic expression $\displaystyle c_{5} = \frac{139}{36} + \frac{118}{27} ln(5)$ (2.19) Now $\displaystyle c_{5}$ in (2.24) can be substituted into (2.22) where $\displaystyle u(x)$ can explicitly be expressed as   $\displaystyle u(x) = -\frac{118}{27} ln|2+3x| - \frac{5}{12}x^2 + \frac{5}{9}x + \frac{139}{36} + \frac{118}{27} ln(5)$ (2.20)

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)

Matlab Code:

MATLAB Code for Problem 5.2

clear all; clc; format long syms x Uh; U = -118/27*log((2+3*x)/5) + (20*x-15*x^2+139)/36; % Exact Solution m = 2; % Global Equidistant Nodes d_0 = 4; x_m = zeros(m); for i=1:m x_m(i)=(i-1)/(m-1); end % Basis Function Vector b_i = sym(ones(m,1)); for i=1:m for j=1:m b_i(i)=(x-1)^(i); end end % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) + subs((12+18*x)*b_i(i),x,0); end d = zeros(m,1); Kff % displaying Kff Ff % displaying Ff d = Kff\Ff; % Solve for d i=1,2,...,m d = [d_0;d] % Prepend d_0

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}}$

Matlab

MATLAB Code for Problem 5.2

clear all; clc; format long syms x Uh; U = -118/27*log((2+3*x)/5) + (20*x-15*x^2+139)/36; % Exact Solution m = 9; % Global Equidistant Nodes d_0 = 4; % Basis Function Vector b_i = sym(ones(m,1)); for i=1:m if rem(i,2)==0 b_i(i)=sin(i/2*(x-1)); else b_i(i)=cos((floor(i/2)+1)*(x-1))-1; end end % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) + subs((12+18*x)*b_i(i),x,0); end d = zeros(m+1,1); Kfff = zeros(m+1,m+1); Kfff(1,1) = 1; Kfff(2:m+1,2:m+1) = Kff% Lazy way of displaying Kff Fff = d_0*ones(m+1,1) Fff(2:m+1,1) =Ff % Lazy way of displaying Ff d = Kfff\Fff; % Solve for d i=1,2,...,m Uh = dot(d,b_i) % This come out simplified with fractions. % You can just look at b_i and d to write out Uh U05=vpa(subs(U,0.5)) Uh05=vpa(subs(Uh,0.5)) error = abs(U05 - Uh05)

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)|$

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:

Matlab

MATLAB Code for Problem 5.2

clear all; clc; format long syms x Uh; U = -118/27*log((2+3*x)/5) + (20*x-15*x^2+139)/36; % Exact Solution m = 30; % Global Equidistant Nodes d_0 = 4; % Basis Function Vector % Basis Function Vector b_i = sym(ones(m,1)); for i=1:m for j=1:m b_i(i)=exp(i*(x-1))-1; end end % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) + subs((12+18*x)*b_i(i),x,0); end d = zeros(m,1); Kff % displaying Kff Ff % displaying Ff d = Kff\Ff; % Solve for d i=1,2,...,m d = [d_0;d] % Prepend d_0 b_i = [1;b_i] % Prepend b_0 Uh = d'*b_i % This come out simplified with fractions. % You can just look at b_i and d to write out Uh U05=vpa(subs(U,0.5)) Uh05=vpa(subs(Uh,0.5)) error = abs(U05 - Uh05)

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:

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

Matlab Code

MATLAB Code for Problem 5.3

clear all; clc; syms x b; m = 4; % Uniform Discretization x_m = zeros(m,1); for i=1:m x_m(i)=(i-1)/(m-1); end % Construct the Lagrange interpolating basis function L = sym(ones(m,1)); for i=1:m for j=1:m if(j~=i) L(i)=L(i)*(x-x_m(j))/(x_m(i)-x_m(j)); L(i) = expand(L(i)); end end end L % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) - 6*subs((2+3*x)*b_i(i),x,1); end d = zeros(m,1); Kff % Lazy way of displaying Kff Ff % Lazy way of displaying Ff d = Kff\Ff; % Solve for d i=1,2,...,m d = [d_0;d] % Prepend d_0 b_i = [1;b_i] % Prepend b_0 Uh = d'*b_i % This come out simplified with fractions. % You can just look at b_i and d to write out Uh U05=vpa(subs(U,0.5)) Uh05=vpa(subs(Uh,0.5)) error = abs(U05 - Uh05) p1 = ezplot(libf(1,m), 1,8); set( p1,'Color','blue', 'LineStyle', '--', 'LineWidth', 2); title('LIBF m = 4') xlabel('Number of Nodes') ylabel('Value of LIBF') hold on p2 = ezplot(libf(2,m), 1,8); set( p2,'Color','red', 'LineStyle', '--', 'LineWidth', 2); p3 = ezplot(libf(3,m), 1,8); set( p3,'Color','green', 'LineStyle', '--', 'LineWidth', 2); p4 = ezplot(libf(4,m), 1,8); set( p4,'Color','cyan', 'LineStyle', '--', 'LineWidth', 2); axis tight p5 = ezplot(libf(5,m), 1,8); set( p5,'Color','magenta', 'LineStyle', '--', 'LineWidth', 2); p6 = ezplot(libf(6,m), 1,8); set( p6,'Color','yellow', 'LineStyle', '--', 'LineWidth', 2); axis tight p7 = ezplot(libf(7,m), 1,8); set( p7,'Color','blue', 'LineStyle', '--', 'LineWidth', 2); p8 = ezplot(libf(8,m), 1,8); set( p8,'Color','red', 'LineStyle', '--', 'LineWidth', 2); axis tight title('LIBF m = 8') xlabel('Number of Nodes') ylabel('Value of LIBF') h = legend('Node1','Node2','Node3','Node4','Node5','Node6','Node7','Node8',-1); hold off

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:

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

Matlab Code

MATLAB Code for Problem 5.3

clear all; clc; syms x b; m = 4; % Uniform Discretization x_m = zeros(m,1); for i=1:m x_m(i)=(i-1)/(m-1); end % Construct the Lagrange interpolating basis function L = sym(ones(m,1)); for i=1:m for j=1:m if(j~=i) L(i)=L(i)*(x-x_m(j))/(x_m(i)-x_m(j)); L(i) = expand(L(i)); end end end L % Construct K and F Kff=ones(m,m); Ff=ones(m,1); for i=1:m for j=1:m Kff(i,j)=int(diff(b_i(i))*(2+3*x)*diff(b_i(j)),0,1); end Ff(i)=int(b_i(i)*5*x,0,1) - 6*subs((2+3*x)*b_i(i),x,1); end d = zeros(m,1); Kff % Lazy way of displaying Kff Ff % Lazy way of displaying Ff d = Kff\Ff; % Solve for d i=1,2,...,m d = [d_0;d] % Prepend d_0 b_i = [1;b_i] % Prepend b_0 Uh = d'*b_i % This come out simplified with fractions. % You can just look at b_i and d to write out Uh U05=vpa(subs(U,0.5)) Uh05=vpa(subs(Uh,0.5)) error = abs(U05 - Uh05) p1 = ezplot(libf(1,m), 1,8); set( p1,'Color','blue', 'LineStyle', '--', 'LineWidth', 2); title('LIBF m = 4') xlabel('Number of Nodes') ylabel('Value of LIBF') hold on p2 = ezplot(libf(2,m), 1,8); set( p2,'Color','red', 'LineStyle', '--', 'LineWidth', 2); p3 = ezplot(libf(3,m), 1,8); set( p3,'Color','green', 'LineStyle', '--', 'LineWidth', 2); p4 = ezplot(libf(4,m), 1,8); set( p4,'Color','cyan', 'LineStyle', '--', 'LineWidth', 2); axis tight p5 = ezplot(libf(5,m), 1,8); set( p5,'Color','magenta', 'LineStyle', '--', 'LineWidth', 2); p6 = ezplot(libf(6,m), 1,8); set( p6,'Color','yellow', 'LineStyle', '--', 'LineWidth', 2); axis tight p7 = ezplot(libf(7,m), 1,8); set( p7,'Color','blue', 'LineStyle', '--', 'LineWidth', 2); p8 = ezplot(libf(8,m), 1,8); set( p8,'Color','red', 'LineStyle', '--', 'LineWidth', 2); axis tight title('LIBF m = 8') xlabel('Number of Nodes') ylabel('Value of LIBF') h = legend('Node1','Node2','Node3','Node4','Node5','Node6','Node7','Node8',-1); hold off

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$.

b.Plots of LLEBF for nel=4.

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

Matlab Code to find approximate solution

Matlab Code: clc clear all close all for y=1:1:4 N=2*y; N_dof(y)=2*y+1; intial=0; syms z digits(15) for i=1:N+1 x(i)=intial; intial=intial+(1/N); u(i)=vpa(241/54*log((2+3*x(i))/2)+(20*x(i)-15*x(i)^2+144)/36);% Exact Solution end k=zeros(N+1,N+1); a_2=2+3*z; F=zeros(N+1,1); for e=1:N for r=1:N+1 b(r)=0*z; end b_d=zeros(N+1,1); b(e)=b(e)+((x(e+1)-z)/(x(e+1)-x(e)));%calculation of basis matrix b(e+1)=b(e+1)+((z-x(e))/(x(e+1)-x(e))); b_d(e)=b_d(e)+vpa(-1/(x(e+1)-x(e)));% derivatives of basis matrix b_d(e+1)=b_d(e+1)+vpa(1/(x(e+1)-x(e))); s=x(e); t=x(e+1); p=b(e)*5*z; q1=b(e+1)*5*z; F(e)=F(e)+vpa(int(p,s,t));%calculations of force matrix F(e+1)=F(e+1)+vpa(int(q1,s,t)); for i=1:N+1 for j=1:N+1 q=b_d(i)*a_2*b_d(j); k(i,j)=k(i,j)+vpa(int(q,z,s,t));% calculation of stiffness matrix end end end %N_2(N/2)=N_2+1; F(N+1)=F(N+1)+12;%inclusion of natural boundry condition in weak form k(1,1)=1e10;% Using Penality method F(1)=k(1,1)*4; %imposing essential boundry condition u_h=vpa((k)\F);% apporimate displacement err(y)=u((N/2)+1)-u_h((N/2)+1);% calculation of error figure(y) plot(x,u,x,u_h,'*-') legend('exact solution','FEA Solution') xlabel('x cordinate') ylabel('displacement') end figure(y+1) plot(N_dof,err,'o-') xlabel('Degrees of Freedom') ylabel('error u(0)-u_h(0)')

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.

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$

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