User:Fem.tm7/HW3

Problem Description

 * Consider $$\left\{ b_j(x);j=0,1,...,n \right\}$$ when $$\displaystyle b_j(x) = (x+k)^j$$


 * Choose $$\displaystyle k=1$$, such that $$\displaystyle b_j(x=0) \ne 0$$

1. Let n=2 $$\rightarrow$$ ndof = n+1 = 2+1 = 3 with  $$\mathbf {d} = \left\{d_j;j=0,...,n\right\}_{(n+1)\times1}$$.

2. Find 2 equations that enforce the boundary conditions for
 * {| style="width:100%" border="0"

$$u^h(x) = \sum_{j=0}^n{d_jb_j(x)}$$. (3.1.1)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

3. Find 1 more equation (i.e. j = 0,1,2) to solve for $$\mathbf {d} = \left\{d_j\right\}_{3\times1}$$ by projecting the residue, $$\displaystyle P(u^h)$$, on a basis function, $$\displaystyle b_k(x)$$, with k = 0,1,2 such that the additional equation is linearly independent from the equations in part 2.


 * Notes regarding residue:
 * The general form of the residue is defined as $$\displaystyle P(u^h)$$. The exact solution is expressed as $$\displaystyle P(u) = 0$$.


 * $$\displaystyle u^h$$ approximates $$\displaystyle u$$. In other words $$\displaystyle u^h \cong u $$. Thus, $$ P(u^h) \ne 0 \quad \forall x \in \Omega $$ (for all values of x in the domain).
 * Notes regarding projection:
 * $$\displaystyle\int_{\Omega } b_i(x) P\left(u^h(x)\right)dx = 0 \quad\Rightarrow\quad\int_0^1 b_k(x) P(u^h)dx = 0$$
 * $$\displaystyle\mathbf b_i \cdot \mathbf P(\mathbf v) = 0 \quad\Rightarrow\quad \mathbf b_i \cdot \mathbf P(\mathbf v) \ = \  <\mathbf b_i,\mathbf P(\mathbf v)>$$


 * Combining the previous two equations yields
 * {| style="width:100%" border="0"

$$\displaystyle  = \int_0^1 b_k(x) P(u^h)dx$$. (3.1.2) 4. Display 3 equations in matrix form, $$\mathbf K \cdot \mathbf d = \mathbf F$$.
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

5. Solve for $$\mathbf d$$.

6. Construct $$\displaystyle u_n^h(x)$$ and plot versus the exact solution, $$\displaystyle u(x)$$.

7. Repeat steps 1 through 6 for:
 * n = 4
 * n = 6

Given
The data set for the general 1-D model with "simple" boundary conditions (G1DM1.0/D2) is given in Vu-Quoc, lecture 12, page 1.


 * $$\displaystyle \Omega = ]{0,1}[$$


 * $$\displaystyle a_2=2$$


 * $$\displaystyle f=3$$

Static case since $$\displaystyle \frac{\partial^su}{\partial t^s}(x,t)=0$$.

The natural (eq 3.1.3) and essential (ed 3.1.4) boundary conditions are defined as


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

$$-\frac{d{u}^{h}}{dx}(0)=4$$ (3.1.3)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }
 * {| style="width:100%" border="0"

$$\displaystyle u(1)=0$$ (3.1.4)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Solution
For the case when $$\displaystyle n=2$$

Parts 1 & 2

 * $$\left\{ b_j(x);j=0,1,...,n \right\}$$ where $$\displaystyle b_j(x) = (x+k)^j$$


 * $$\Rightarrow \quad b_0(x)=(x+k)^0=1, b_1(x)=(x+k)^1=(x+k),\ and\ b_2(x)=(x+k)^2$$.

Choosing $$k=1$$ to avoid $$b_j'(x)=0$$ yields


 * $$\quad b_0(x)=1, b_1(x)=(x+1),\ and\ b_2(x)=(x+1)^2$$.

The natural boundary condition can be implemented by differentiating $$\displaystyle u^h(x)$$ with respect to $$\displaystyle x$$ and taking its value at x=0.


 * $$\frac{d{u}^{h}}{dx}(x) = \sum_{j=0}^2{d_jb_j'(x)}=d_0b_0'(x)+d_1b_1'(x)+d_2b_2'(x)=d_1+2d_2(x+1)$$
 * $$\frac{d{u}^{h}}{dx}(0) = \sum_{j=0}^2{d_jb_j'(0)}=d_0b_0'(0)+d_1b_1'(0)+d_2b_2'(0)=d_1+2d_2=-4$$

So the relevant equation is


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

$$ \displaystyle d_1+2d_2=-4$$ (3.1.5)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * 
 * }

The essential boundary condition, $$\displaystyle u^h(1)=0$$, is implemented as follows:


 * $$u^h(x)= \sum_{j=0}^2{d_jb_j(x)}=d_0b_0(x)+d_1b_1(x)+d_2b_2(x)=d_0+d_1(x+1)+d_2(x+1)^2$$
 * $$u^h(1)= \sum_{j=0}^2{d_jb_j(1)}=d_0b_0(1)+d_1b_1(1)+d_2b_2(1)=d_0+2d_1+4d_2=0$$

So the relevant equation is


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

$$ \displaystyle d_0+2d_1+4d_2=0$$ (3.1.6)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * 
 * }

Part 3
Equation 3.1.2 can be used to project the residue onto the basis function. First, the partial differential equation $$\displaystyle P(u^h)$$ is defined.


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

$$ \displaystyle P(u^h) = \frac{\partial }{\partial x}\left[a_2(x)\frac{\partial u}{\partial x}\right]+f(x,t) - \overline {m} (x) \frac{\partial ^s u}{\partial t^s}$$ (3.1.7)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Substituting the given conditions into this equation and recalling equation 3.1.1,


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

$$ \displaystyle P(u^h) = \frac{\partial }{\partial x}\left\{ 2 \left[d_1+2d_2(x+1)\right] \right\}+3$$ (3.1.8)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Differentiating with respect to $$\displaystyle x$$ and substituting known values,


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

$$\displaystyle P(u^h) = \left[4d_2\right]+3$$ (3.1.9)
 * style="width:10%; |
 * style="width:10%; |
 * 
 * }

Substituting into 3.1.2


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

$$ \displaystyle  = \int_0^1 b_k(x) P(u^h)dx = \int_0^1 b_k(x) \left[4d_2 +3\right]dx=0$$ (3.1.10)
 * style="width:10%;|
 * style="width:10%;|
 * 
 * }

Since the equation is valid for all $$\displaystyle b_k(x), \ b_0(x)=(x+1)$$ is utilized.


 * $$\displaystyle  = \int_0^1 (x+1) \left[4d_2 +3\right]dx= \int_0^1 [(3+4d_2)x+4d_2+3]dx=(3+4d_2)\frac{1^2}{2}+(4d_2+3)1 = 6d_2+\frac{9}{2}=0$$


 * $$\displaystyle $$

Simplifying, the relevant equation becomes


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

$$ \displaystyle 4d_2=-3$$ (3.1.11)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * 
 * }

Part 4
The coefficient ($$K$$) and constant ($$f$$) matrices are constructed from the equations determined aboved.


 * Notes:
 * The first row of the matrix $$\mathbf K$$ is based on equation 3.1.1 and the essential boundary condition, 3.1.4, when $$\displaystyle x=1$$.


 * The second row of the matrix $$\mathbf K$$ is from equations 3.1.5 and 3.1.6 (for the general case).


 * The third row of the matrix $$\mathbf K$$ is from equation 3.1.11.


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

$$ \begin{bmatrix} 0 & 1 & 2\\ 1& 2 & 4\\  0& 0 & 4 \end{bmatrix} \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |

Note $$\displaystyle\mathbf K$$ is not Symmetric. (3.1.12)
 * 
 * }

Part 5
To solve for the matrix $$\mathbf d$$, first recognize that $$\mathbf K \cdot \mathbf d = \mathbf F$$ can be rewritten as $$\mathbf d = \mathbf K^{-1} \cdot \mathbf F$$.


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

$$ \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 2\\ 1& 2 & 4\\  0& 0 & 4 \end{bmatrix}^{-1} \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix}$$
 * style="width:10%; |
 * style="width:10%; |

$$ \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} -2 & 1 & 0\\ 1& 0 & -0.5\\  0& 0 & 0.25 \end{bmatrix} \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix} $$


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

$$ \displaystyle \underline {d}_{n=2} = \begin{bmatrix} 8\\ -2.5\\ -0.75 \end{bmatrix} $$     (3.1.13)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * 
 * }
 * }

The inverse matrix and d coefficients were found using the following matlab script:

Parts 6 & 7
The approximate function, $$\displaystyle u_n^h(x)$$, can now be written using the determined coefficients in equation 3.1.13. The exact solution, was obtained and is displayed here:



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

The same procedure is repeated for $$n = 4$$ and $$n = 6$$ in equation 3.1.1.

As the procedure is the same shown above, a Matlab script was written to enforce the boundary conditions, perform the inner products and integration, solve the system of equations, and plot the results.

Figure 3.1.1 shows the approximate solutions for $$n=2,4,6$$ plotted against the exact solution in the domain $$ \Omega $$. It is clear that $$n=2$$ provides a good approximation. This is likely due to the fact that the approximate and exact solutions are both second order polynomial functions.



Part 8
The previous Matlab code also calculates the error in the estimate for u or the following equation:


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

$$\displaystyle e_n(x)=u(x)-u^h(x)$$ (3.1.14)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

This equation was evaluated at .5 and plotted vs n giving the following figure:



Author
Johnathan Whittaker Bullard

Problem Description
Consider figure 2.16 from "A First Course in Finite Elements" by Fish and Belytschko p. 37 pb.2.1




 * Part a: Number the elements and nodes.


 * Part b: Assemble the global stiffness and force matrix.


 * Part c: Partition the system and solve for the nodal displacements.


 * Part d: Compute the reaction forces.

Part a



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


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

The elements are numbered in $$ \displaystyle {\color{red} red} $$ and the nodes are labled in $$ \displaystyle {\color{blue} blue} $$
 * style = |
 * }
 * }

Part b
Assemble the global stiffness and force matrix.

$$ \displaystyle {\color{red} Element 1} $$,  $$ \displaystyle {\color{blue} I=1} $$,  $$ \displaystyle {\color{blue} J=4}$$

$$ \displaystyle \mathbf{K}^{(1)}=3k \begin{vmatrix} 1 & -1\\ -1 & 1 \end{vmatrix} \Rightarrow \mathbf{\tilde{K}}^{(1)}= \begin{vmatrix} 3k & 0 & 0 & -3k\\ 0 & 0 & 0 &0 \\ 0 & 0  & 0 &0 \\ -3k &0  &0  & 3k \end{vmatrix}$$

$$ \displaystyle {\color{red} Element 2} $$,  $$ \displaystyle {\color{blue} I=1} $$,  $$ \displaystyle {\color{blue} J=3}$$

$$ \displaystyle \mathbf{K}^{(2)}=k \begin{vmatrix} 1 & -1\\ -1 & 1 \end{vmatrix} \Rightarrow \mathbf{\tilde{K}}^{(2)}= \begin{vmatrix} k & 0 & -k & 0\\ 0 & 0 & 0 &0 \\ -k &0  & k &0 \\ 0 &0 &0  & 0 \end{vmatrix}$$

$$ \displaystyle {\color{red} Element 3} $$,  $$ \displaystyle {\color{blue} I=3} $$,  $$ \displaystyle {\color{blue} J=4}$$

$$ \displaystyle \mathbf{K}^{(3)}=2k \begin{vmatrix} 1 & -1\\ -1 & 1 \end{vmatrix} \Rightarrow \mathbf{\tilde{K}}^{(3)}= \begin{vmatrix} 0 & 0 & 0 & 0\\ 0& 0 & 0  &0 \\ 0 & 0  & 2k &-2k \\ 0 &0 & -2k  & 2k \end{vmatrix}$$

$$ \displaystyle {\color{red} Element 4} $$,  $$ \displaystyle {\color{blue} I=4} $$,  $$ \displaystyle {\color{blue} J=2}$$

$$ \displaystyle \mathbf{K}^{(4)}=k \begin{vmatrix} 1 & -1\\ -1 & 1 \end{vmatrix} \Rightarrow \mathbf{\tilde{K}}^{(4)}= \begin{vmatrix} 0 & 0 & 0 & 0\\ 0 & k & 0 & -k \\ 0 & 0 & 0 & 0 \\ 0 & -k & 0  & k \end{vmatrix}$$

Assembled system matrix:


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


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$ \displaystyle \mathbf{K}=\sum_{e=1}^{4}\mathbf{\tilde{K}}^{e}=\begin{vmatrix} 4k & 0 & -k & -3k\\ 0 & k & 0 & -k\\ -k & 0 & 3k & -2k\\ -3k & -k & -2k & 6k \end{vmatrix} $$


 * style = |
 * }
 * }

The displacement and force matrices for the system are:

$$ \displaystyle \mathbf{d}=\begin{vmatrix} 0\\ 0\\ u_{3}\\ u_{4}\\ \end{vmatrix},\;

\mathbf{f}=\begin{vmatrix} 0\\ 0\\ 0\\ 50\\ \end{vmatrix}N,\;

\mathbf{r}=\begin{vmatrix} r_{1}\\ r_{2}\\ 0\\ 0\\ \end{vmatrix}\; $$

The global system of equations is given by:


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


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$ \displaystyle \mathbf{K}=\sum_{e=1}^{4}\mathbf{\tilde{K}}^{e}=\begin{vmatrix} 4k & 0 & -k & -3k\\ 0 & k & 0 & -k\\ -k & 0 & 3k & -2k\\ -3k & -k & -2k & 6k \end{vmatrix} \displaystyle \begin{vmatrix} 0\\ 0\\ u_{3}\\ u_{4}\\ \end{vmatrix}= \begin{vmatrix} r_{1}\\ r_{2}\\ 0\\ 50\\ \end{vmatrix} $$


 * style = |
 * }
 * }

Part c
Partition the system and solve for the nodal displacements.

As the first two displacements are prescribed, we partition after two rows and columns:

$$ \displaystyle \begin{vmatrix} \mathbf{K}_{E} & \mathbf{K}_{EF}\\ \mathbf{K}^{T}_{EF} & \mathbf{K}_{F} \end{vmatrix}

\begin{vmatrix} \mathbf{\overline{d}}_{E}\\ \mathbf{d}_{F} \end{vmatrix}=

\begin{vmatrix} \mathbf{r}_{E}\\ \mathbf{f}_{F} \end{vmatrix} $$

Where

$$ \displaystyle \mathbf{K}_{E}= \begin{vmatrix} 4k & 0 \\ 0 & k \end{vmatrix} $$, $$ \displaystyle \mathbf{K}_{F}= \begin{vmatrix} 3k & 2k \\ -2k & 6k \end{vmatrix} $$,  $$ \displaystyle \mathbf{K}_{EF}= \begin{vmatrix} -k & -3k \\ 0 & -k \end{vmatrix} $$,  $$ \displaystyle \mathbf{\overline{d}}_{E}= \begin{vmatrix} 0\\ 0 \end{vmatrix} $$,  $$ \displaystyle \mathbf{d}_{f}= \begin{vmatrix} u_{3}\\ u_{4} \end{vmatrix} $$,  $$ \displaystyle \mathbf{f}_{F}= \begin{vmatrix} 0\\ 50 \end{vmatrix} $$,  $$ \displaystyle \mathbf{e}_{E}= \begin{vmatrix} r_{1}\\ r_{2} \end{vmatrix} $$

Solving for the nodal displacements: $$ \displaystyle \begin{vmatrix} \mathbf{K}_{F}\\ \end{vmatrix}

\begin{vmatrix} \mathbf{d}_{F}\\ \end{vmatrix}=

\begin{vmatrix} \mathbf{f}_{F}\\ \end{vmatrix} $$

$$ \displaystyle \begin{vmatrix} 3k & -2k \\ -2k & 6k \end{vmatrix}

\begin{vmatrix} u_{3}\\ u_{4} \end{vmatrix}=

\begin{vmatrix} 0 \\ 50 \end{vmatrix} $$

The systems of equations are: $$ \displaystyle 3k\cdot u_{3}-2k\cdot u_{4}=0 $$ Muliply by -3 $$ \displaystyle -3(3k\cdot u_{3}-2k\cdot u_{4})=0 $$ $$ \displaystyle -9k\cdot u_{3}-6k\cdot u_{4}=0 $$

and $$ \displaystyle -2k\cdot u_{3}-6k\cdot u_{4}=50 $$ Summing the two equations lets us solve for the displacements as:


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


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$ \displaystyle u_{3}=\frac{-50}{7k},\; u_{4}=\frac{75}{7k} $$


 * style = |
 * }
 * }

Part d
Compute the reaction forces: The global system now becomes:

$$ \displaystyle \mathbf{K}=\sum_{e=1}^{4}\mathbf{\tilde{K}}^{e}=\begin{vmatrix} 4k & 0 & -k & -3k\\ 0 & k & 0 & -k\\ -k & 0 & 3k & -2k\\ -3k & -k & -2k & 6k \end{vmatrix} \displaystyle \begin{vmatrix} 0\\ 0\\ \frac{-50}{11k} \\ \frac{75}{11k} \\ \end{vmatrix}= \begin{vmatrix} r_{1}\\ r_{2}\\ 0\\ 50\\ \end{vmatrix} $$

Solving for the reation forces:

$$ \displaystyle \begin{vmatrix} \mathbf{K}_{EF}\\ \end{vmatrix} \begin{vmatrix} \mathbf{d}_{F}\\ \end{vmatrix}= \begin{vmatrix} \mathbf{r}_{E}\\ \end{vmatrix} $$

The equations for the reaction forces are: $$ \displaystyle -k\left (\frac{-50}{11k} \right )-3k\left (\frac{75}{11k}  \right )=r_{1} $$

$$ \displaystyle -k\left (\frac{75}{11k} \right )=r_{2} $$

The reaction forces equal:


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


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$ \displaystyle r_{1}=\frac{-175}{11}N $$ and $$ r_{2}=\frac{75}{11}N $$
 * style = |
 * }
 * }

Author
Brandonhua 04:48, 15 February 2011 (UTC)

Given
Looking at the truss structure below where nodes A and B are fixed. A 10 N force in the positive x-direction acts at node C. The joint locations are in meters. The corresponding Young ‘ s modulus is $$ E = {10^{11}}{\text{Pa}} $$ and cross-sectional area for all bars are $$ A = 2 \cdot {10^{ - 2}}{{\text{m}}^2} $$.



Find

 * 1) Number the elements and nodes.
 * 2) Assemble the global stiffness and force matrix.
 * 3) Partition the system and solve for the nodal displacements.
 * 4) Compute the stresses and reactions.

Solution
1. The elements and nodes numbers are shown in the figure below:



2. Dividing the structure into 4 elements and 4 nodes, and deals with each element starting with element 1: Element 1 is numbered with global nodes 1 and 4. It is positioned at an angle $$ {\phi ^{\left( 1 \right)}} = {90^ \circ } $$with respect to positive x-axis. Then,
 * {| style="width:100%" border="0"

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

\cos {90^ \circ } = 0

$$     (Eq 3.1)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\sin {90^ \circ } = 1

$$     (Eq 3.2)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{l^{\left( 1 \right)}} = 1{\text{m}}

$$     (Eq 3.3)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{k^{\left( 1 \right)}} = \frac = 2 \times {10^9}{\text{N/m}}

$$     (Eq 3.4)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{K}}^{\left( 1 \right)}} = 2 \times {10^9}\left[ {\begin{array}{ccccccccccccccc} 0&0&0&0 \\  0&1&0&{ - 1} \\   0&0&0&0 \\   0&{ - 1}&0&1 \end{array}} \right]

$$     (Eq 3.5) Element 2 is numbered with global nodes 2 and 4. It is positioned at an angle $${\phi ^{\left( 2 \right)}} = {135^ \circ }  $$with respect to positive x-axis.
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\cos {135^ \circ } = - \frac{1}

$$     (Eq 3.6)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\sin {135^ \circ } = \frac{1}

$$     (Eq 3.7)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{l^{\left( 2 \right)}} = \sqrt 2 {\text{m}}

$$     (Eq 3.8)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{k^{\left( 2 \right)}} = \frac = \frac{2} \times {10^9}{\text{N/m}}

$$     (Eq 3.9)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{K}}^{\left( 2 \right)}} = \frac{2} \times {10^9}\left[ {\begin{array}{ccccccccccccccc} {\frac{1}{2}}&{ - \frac{1}{2}}&{ - \frac{1}{2}}&{\frac{1}{2}} \\ { - \frac{1}{2}}&{\frac{1}{2}}&{\frac{1}{2}}&{ - \frac{1}{2}} \\ { - \frac{1}{2}}&{\frac{1}{2}}&{\frac{1}{2}}&{ - \frac{1}{2}} \\ {\frac{1}{2}}&{ - \frac{1}{2}}&{ - \frac{1}{2}}&{\frac{1}{2}} \end{array}} \right]

$$     (Eq 3.10) Element 3 is numbered with global nodes 3 and 4. It is positioned at an angle $$ {\phi ^{\left( 3 \right)}} = {180^ \circ } $$with respect to positive x-axis.
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\cos {180^ \circ } = - 1

$$     (Eq 3.11)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\sin {180^ \circ } = 0

$$     (Eq 3.12)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{l^{\left( 3 \right)}} = 1{\text{m}}

$$     (Eq 3.13)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{k^{\left( 3 \right)}} = \frac = 2 \times {10^9}{\text{N/m}}

$$     (Eq 3.14)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{K}}^{\left( 2 \right)}} = 2 \times {10^9}\left[ {\begin{array}{ccccccccccccccc} 1&0&{ - 1}&0 \\  0&0&0&0 \\   { - 1}&0&1&0 \\   0&0&0&0 \end{array}} \right]

$$     (Eq 3.15) Element 4 is numbered with global nodes 2 and 3. It is positioned at an angle $$ {\phi ^{\left( 4 \right)}} = {90^ \circ } $$with respect to positive x-axis.
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\cos {90^ \circ } = 0

$$     (Eq 3.16)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\sin {90^ \circ } = 1

$$     (Eq 3.17)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{l^{\left( 4 \right)}} = 1{\text{m}}

$$     (Eq 3.18)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{k^{\left( 4 \right)}} = \frac = 2 \times {10^9}{\text{N/m}}

$$     (Eq 3.19)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{K}}^{\left( 4 \right)}} = 2 \times {10^9}\left[ {\begin{array}{ccccccccccccccc} 0&0&0&0 \\  0&1&0&{ - 1} \\   0&0&0&0 \\   0&{ - 1}&0&1 \end{array}} \right]

$$     (Eq 3.20) Assemble the global stiffness matrix
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\mathbf{K}} = 2 \times {10^9}\left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0&0&0&0 \\  0&1&0&0&0&0&0&{ - 1} \\   0&0&{\frac{1}}&{ - \frac{1}}&0&0&{ - \frac{1}}&{\frac{1}} \\ 0&0&{ - \frac{1}}&{1 + \frac{1}}&0&{ - 1}&{\frac{1}}&{ - \frac{1}} \\ 0&0&0&0&1&0&{ - 1}&0 \\  0&0&0&{ - 1}&0&1&0&0 \\   0&0&{ - \frac{1}}&{\frac{1}}&{ - 1}&0&{1 + \frac{1}}&{ - \frac{1}} \\ 0&{ - 1}&{\frac{1}}&{ - \frac{1}}&0&0&{ - \frac{1}}&{1 + \frac{1}} \end{array}} \right]

$$     (Eq 3.21) The external force and reaction matrix
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\mathbf{f}} = \left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   0 \\   0 \\   {10} \\   0 \\   0 \\   0 \end{array}} \right]

$$     (Eq 3.22)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\mathbf{r}} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\    \\   0 \\   0 \\   0 \\   0 \end{array}} \right]

$$     (Eq 3.23) And the displacement matrix
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\mathbf{d}} = \left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   0 \\   0 \\    \\    \\    \\ \end{array}} \right]

$$     (Eq 3.24) 3. Global system of equation
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

2 \times {10^9}\left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0&0&0&0 \\  0&1&0&0&0&0&0&{ - 1} \\   0&0&{\frac{1}}&{ - \frac{1}}&0&0&{ - \frac{1}}&{\frac{1}} \\ 0&0&{ - \frac{1}}&{1 + \frac{1}}&0&{ - 1}&{\frac{1}}&{ - \frac{1}} \\ 0&0&0&0&1&0&{ - 1}&0 \\  0&0&0&{ - 1}&0&1&0&0 \\   0&0&{ - \frac{1}}&{\frac{1}}&{ - 1}&0&{1 + \frac{1}}&{ - \frac{1}} \\ 0&{ - 1}&{\frac{1}}&{ - \frac{1}}&0&0&{ - \frac{1}}&{1 + \frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   0 \\   0 \\    \\    \\    \\ \end{array}} \right] = {\mathbf{r}} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\    \\   {10} \\   0 \\   0 \\   0 \end{array}} \right]

$$     (Eq 3.25) Partition the system
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\overline {\mathbf{d}} _E} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   0 \\   0 \end{array}} \right],{{\mathbf{d}}_F} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right],{{\mathbf{f}}_F} = \left[ {\begin{array}{ccccccccccccccc} {10} \\  0 \\   0 \\   0 \end{array}} \right],{{\mathbf{K}}_F} = \left[ {\begin{array}{ccccccccccccccc} 1&0&{ - 1}&0 \\  0&1&0&0 \\   { - 1}&0&{1 + \frac{1}}&{ - \frac{1}} \\ 0&0&{ - \frac{1}}&{1 + \frac{1}} \end{array}} \right]

$$
 * <p style="text-align:right">
 * }

$$ \displaystyle

{{\mathbf{r}}_E} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right],{{\mathbf{K}}_{EF}} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0 \\  0&0&0&{ - 1} \\   0&0&{ - \frac{2}}&{\frac{2}} \\ 0&{ - 1}&{\frac{2}}&{ - \frac{2}} \end{array}} \right]

$$ <p style="text-align:right"> Solve for the nodal displacements
 * {| style="width:100%" border="0"

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

2 \times {10^9}\left[ {\begin{array}{ccccccccccccccc} 1&0&{ - 1}&0 \\  0&1&0&0 \\   { - 1}&0&{1 + \frac{1}}&{ - \frac{1}} \\ 0&0&{ - \frac{1}}&{1 + \frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} {10} \\  0 \\   0 \\   0 \end{array}} \right]

$$     (Eq 3.26)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \frac{1}{2} \times {10^{ - 9}}\left[ {\begin{array}{ccccccccccccccc} {48.2843} \\  0 \\   {38.2843} \\   {10} \end{array}} \right]

$$     (Eq 3.27) 4. The reaction matrix is
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{r}}_E} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = {{\mathbf{K}}_E}{{\mathbf{\bar d}}_E} + {{\mathbf{K}}_{EF}}{{\mathbf{d}}_F} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0 \\  0&0&0&{ - 1} \\   0&0&{ - \frac{1}}&{\frac{1}} \\ 0&{ - 1}&{\frac{1}}&{ - \frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {48.2843} \\  0 \\   {38.2843} \\   {10} \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} 0 \\  { - 10} \\   { - 10} \\   {10} \end{array}} \right]

$$     (Eq 3.28) The stresses in the two elements are
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma ^e} = \frac\left[ {\begin{array}{ccccccccccccccc} { - \cos {\phi ^e}}&{ - \sin {\phi ^e}}&{\cos {\phi ^e}}&{\sin {\phi ^e}} \end{array}} \right]{{\mathbf{d}}^e}

$$     (Eq 3.29) For element 1:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{d}}^{\left( 1 \right)}} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \frac{1}{2} \times {10^{ - 9}}\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   {38.2843} \\   {10} \end{array}} \right]

$$     (Eq 3.30)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma ^{\left( 1 \right)}} = \left[ {\begin{array}{ccccccccccccccc} 0&{ - 1}&0&1 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   {38.2843} \\   {10} \end{array}} \right]\frac{1} = 500{\text{Pa}}

$$     (Eq 3.31) For element 2:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{d}}^{\left( 2 \right)}} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \frac{1}{2} \times {10^{ - 9}}\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   {38.2843} \\   {10} \end{array}} \right]

$$     (Eq 3.32)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma ^{\left( 2 \right)}} = \left[ {\begin{array}{ccccccccccccccc} {\frac{1}}&{ - \frac{1}}&{ - \frac{1}}&{\frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   {38.2843} \\   {10} \end{array}} \right]\frac{1} =  - 1000{\text{Pa}}

$$     (Eq 3.33) For element 3:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{d}}^{\left( 3 \right)}} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \frac{1}{2} \times {10^{ - 9}}\left[ {\begin{array}{ccccccccccccccc} {48.2843} \\  0 \\   {38.2843} \\   {10} \end{array}} \right]

$$     (Eq 3.34)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma ^{\left( 3 \right)}} = \left[ {\begin{array}{ccccccccccccccc} 1&0&{ - 1}&0 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {48.2843} \\  0 \\   {38.2843} \\   {10} \end{array}} \right]\frac{1} = 500{\text{Pa}}

$$     (Eq 3.35) For element 4:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{{\mathbf{d}}^{\left( 4 \right)}} = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \frac{1}{2} \times {10^{ - 9}}\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   {48.2843} \\   0 \end{array}} \right]

$$     (Eq 3.36)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma ^{\left( 4 \right)}} = \left[ {\begin{array}{ccccccccccccccc} 0&{ - 1}&0&1 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   {48.2843} \\   0 \end{array}} \right]\frac{1} = 0

$$     (Eq 3.37)
 * <p style="text-align:right">
 * }

Author
Jiang Jin

Problem Description

 * Prove the asymmetry of the weighted residual stiffness matrix by showing:


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

$$ \displaystyle \int \mathbf b_j(x)\big [\frac{d}{dx}a_2(x)\frac{d}{dx}\mathbf b_i(x)]dx\neq\int \mathbf b_i(x)\big [\frac{d}{dx}a_2(x)\frac{d}{dx}\mathbf b_j(x)]dx$$ (3.4.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Given
The stiffness matrix is defined as:
 * {| style="width:100%" border="0"

$$ \displaystyle K_{i,j}=<\mathbf b_i,\mathbf P(u^h)>$$ (3.4.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Where < > is for the inner product. Here uh is an approximation of u from a summation of basis functions given as:


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

$$ \displaystyle u(x)\approx u^h(x)=\sum_j \mathbf b_j(x)$$ (3.4.3)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Also, P(u) is given by:


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

$$ \displaystyle \mathbf P(\mathbf u(x))=\frac{d}{dx}a_2(x)\frac{d}{dx}\mathbf u(x)$$ (3.4.4)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Plugging Eq. 3.4.3 into Eq. 3.4.4 and looking only at one term of the summation gives:


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

$$ \displaystyle \frac{d}{dx} a_2(x)\frac{d}{dx} \mathbf b_j(x)=(a_2*\mathbf b_j')'$$ (3.4.5) Which for a different numbered basis function than bj(x) is written as:
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \frac{d}{dx} a_2(x)\frac{d}{dx} \mathbf b_i(x)=(a_2*\mathbf b_i')'$$ (3.4.6)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Plugging 3.4.5 and 3.4.6 into 3.4.2 gives:


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

$$ \displaystyle \int \mathbf b_j(x)\big [\frac{d}{dx}a_2(x)\frac{d}{dx}\mathbf b_i(x)]dx$$ (3.4.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \int \mathbf b_i(x)\big [\frac{d}{dx}a_2(x)\frac{d}{dx}\mathbf b_j(x)]dx$$ (3.4.8)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Which our are equations for Eq.3.4.1, so we now must prove:


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

$$ \displaystyle K_{i,j}\neq K_{j,i}$$ (3.4.9)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Solution
First, expanding out the integrands of Eq. 3.4.1:


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

$$ \displaystyle \alpha = \mathbf b_i\big [a_2'*\mathbf b_j'+a_2*\mathbf b_j''\big ]$$ (3.4.10)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \beta = \mathbf b_j\big [a_2'*\mathbf b_i'+a_2*\mathbf b_i''\big ]$$ (3.4.11)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

So now we're simply proving:


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

$$ \displaystyle \alpha \neq \beta$$ (3.4.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

This will be proven using the following basis functions:


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

$$ \displaystyle \mathbf b_i(x)=cos(i*x)$$ (3.4.13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \mathbf b_j(x)=cos(j*x)$$ (3.4.14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Where $$\displaystyle i\neq j$$.

Plugging Eq. 3.4.13 and 3.4.14 into Eq. 3.4.10 gives:


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

$$ \displaystyle \alpha = cos(i*x)\big [a_2'(-jsin(j*x))+a_2(-j^2cos(j*x))\big]$$ (3.4.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Repeating for Eq. 3.4.11:


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

$$ \displaystyle \beta = cos(j*x)\big [a_2'(-isin(i*x))+a_2(-i^2cos(i*x))\big]$$ (3.4.16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Expanding out Eq. 3.4.15 and 3.4.16:


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

$$ \displaystyle \alpha=-ja_2'cos(i*x)sin(j*x)-j^2a_2cos(i*x)cos(j*x)$$ (3.4.17)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \beta=-ia_2'cos(j*x)sin(i*x)-i^2a_2cos(j*x)cos(i*x)$$ (3.4.18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Inspecting Eq. 3.4.17 and Eq. 3.4.18 shows:


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

$$ \displaystyle \alpha\neq\beta\text{ when } i\neq j $$ (3.4.19)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Author
Braden Snook

Problem Statement
Jacob Fish and Ted Belytschko, " A first Course in Finite Elements",John Wiley & Sons, Ltd, Chapter 2, P37 Problem 2.2. " Show that the Equivalent stiffness of a spring aligned in the x-direction for the bar of thickness t with a centered square hole in figure is

Where E is the Young's modulus and t is the width of the bar." Reference figure in textbook: 2.17

Solution
The bar can be sub-divided into 3 elements. The node numbers 1,2,3,4 are assgined. The element have stiffness $${\color{red} k^{\left( 1\right)}}, {\color{green} k^{\left( 2\right)}}, {\color{blue} k^{\left(3\right)}} $$. We know that stiffness $$ \mathbf{k=\frac{AE}{L}} $$ Where A is cross sectional area, E is the Elastic modulus and L is the length of the element

Now assembling the stiffness matrices we get,

Solving the above we get,

Author
Srilalithkumar 15:41, 15 February 2011 (UTC)

Problem Description
"A First Course in Finite Elements" by Fish and Belytschko p. 38 Problem 2.4

Given the three-bar structure subjected to the prescribed load ad point C equal to $${10^3}N$$ as shown. The Young’s modulus is $$E={{10}^{11}}$$ Pa. The cross-sectional area of the bar BC is $${{10}^{-2}}{{m}^{2}}$$ and that of BD and BF is $${{10}^{-2}}{{m}^{2}}$$. Note that point D is free to move in the x-direction. Coordinates of joins are given in meters.


 * 1) Construct the global stiffness matrix and load matrix.
 * 2) Partition the matrices and solve for the unknown displacement at Point B and displacement in the x-direction at point D.
 * 3) Find the stressed in the three bars.
 * 4) Find the reaction at the nodes C,D and F.

Solution
a) First subdivide the structure into elements. There are three elements and two nodes. Shown in Fig. Element 1:
 * {| style="width:100%" border="0"

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

\varphi = 135^\circ ,\cos \varphi  =  - \frac{1},\sin \varphi  = \frac{1}

$$
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{K^{(1)}} = \frac\left[ {\begin{array}{ccccccccccccccc} {\frac{1}{2}}&{ - \frac{1}{2}}&{ - \frac{1}{2}}&{\frac{1}{2}} \\ { - \frac{1}{2}}&{\frac{1}{2}}&{\frac{1}{2}}&{ - \frac{1}{2}} \\ { - \frac{1}{2}}&{\frac{1}{2}}&{\frac{1}{2}}&{ - \frac{1}{2}} \\ {\frac{1}{2}}&{ - \frac{1}{2}}&{ - \frac{1}{2}}&{\frac{1}{2}} \end{array}} \right]

$$     (3.6.1) Element 2:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\varphi = 90^\circ ,\cos \varphi  = 0,\sin \varphi  = 1

$$
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{K^{(2)}} = \frac{1}\left[ {\begin{array}{ccccccccccccccc} 0&0&0&0 \\  0&1&0&{ - 1} \\   0&0&0&0 \\   0&{ - 1}&0&1 \end{array}} \right]

$$     (3.6.2) Element 3:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\varphi = 45^\circ ,\cos \varphi  = \frac{1},\sin \varphi  = \frac{1}

$$
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{K^{(3)}} = \frac\left[ {\begin{array}{ccccccccccccccc} {\frac{1}{2}}&{\frac{1}{2}}&{ - \frac{1}{2}}&{ - \frac{1}{2}} \\ {\frac{1}{2}}&{\frac{1}{2}}&{ - \frac{1}{2}}&{ - \frac{1}{2}} \\ { - \frac{1}{2}}&{ - \frac{1}{2}}&{\frac{1}{2}}&{\frac{1}{2}} \\ { - \frac{1}{2}}&{ - \frac{1}{2}}&{\frac{1}{2}}&{\frac{1}{2}} \end{array}} \right]

$$     (3.6.3) The global stiffness matrix is:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

K = \frac{1}\left[ {\begin{array}{ccccccccccccccc} {\frac{1}}&{ - \frac{1}}&0&0&0&0&{ - \frac{1}}&{\frac{1}} \\ { - \frac{1}}&{\frac{1}}&0&0&0&0&{\frac{1}}&{ - \frac{1}} \\ 0&0&0&0&0&0&0&0 \\  0&0&0&2&0&0&0&{ - 2} \\   0&0&0&0&{\frac{1}}&{\frac{1}}&{ - \frac{1}}&{ - \frac{1}} \\ 0&0&0&0&{\frac{1}}&{\frac{1}}&{ - \frac{1}}&{ - \frac{1}} \\ { - \frac{1}}&{\frac{1}}&0&0&{ - \frac{1}}&{ - \frac{1}}&{\frac{1}}&0 \\ {\frac{1}}&{ - \frac{1}}&0&{ - 2}&{ - \frac{1}}&{ - \frac{1}}&0&{2 + \frac{1}} \end{array}} \right]

$$     (3.6.4) And
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

d = {\left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&&0&& \end{array}} \right]^T}

$$     (3.6.5)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

F = {\left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0&0&&0 \end{array}} \right]^T}

$$     (3.6.6)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

r = {\left[ {\begin{array}{ccccccccccccccc} &&&&0&&0&0 \end{array}} \right]^T}

$$     (3.6.7) b) Global system of equations:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\frac{1}\left[ {\begin{array}{ccccccccccccccc} {\frac{1}}&{ - \frac{1}}&0&0&0&0&{ - \frac{1}}&{\frac{1}} \\ { - \frac{1}}&{\frac{1}}&0&0&0&0&{\frac{1}}&{ - \frac{1}} \\ 0&0&0&0&0&0&0&0 \\  0&0&0&2&0&0&0&{ - 2} \\   0&0&0&0&{\frac{1}}&{\frac{1}}&{ - \frac{1}}&{ - \frac{1}} \\ 0&0&0&0&{\frac{1}}&{\frac{1}}&{ - \frac{1}}&{ - \frac{1}} \\ { - \frac{1}}&{\frac{1}}&0&0&{ - \frac{1}}&{ - \frac{1}}&{\frac{1}}&0 \\ {\frac{1}}&{ - \frac{1}}&0&{ - 2}&{ - \frac{1}}&{ - \frac{1}}&0&{2 + \frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   0 \\   0 \\    \\   0 \\    \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\    \\   0 \\    \\    \\   0 \end{array}} \right]

$$     (3.6.8) Then reduce the global system of equations: The global system is partitioned four rows and four columns: And:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\frac{1}\left[ {\begin{array}{ccccccccccccccc} {\frac{1}}&{\frac{1}}&{ - \frac{1}}&{ - \frac{1}} \\ {\frac{1}}&{\frac{1}}&{ - \frac{1}}&{ - \frac{1}} \\ { - \frac{1}}&{ - \frac{1}}&{\frac{1}}&0 \\ { - \frac{1}}&{ - \frac{1}}&0&{2 + \frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} \\  0 \\    \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} 0 \\   \\    \\   0 \end{array}} \right]

$$     (3.6.9) So:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\left[ {\begin{array}{ccccccccccccccc} \\  0 \\    \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} {(1 + 2\sqrt 2 ) \times {{10}^{ - 6}}} \\ 0 \\  {(\frac{1}{2} + 2\sqrt 2 ) \times {{10}^{ - 6}}} \\ {0.5 \times {{10}^{ - 6}}} \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} {3.828 \times {{10}^{ - 6}}} \\ 0 \\  {3.328 \times {{10}^{ - 6}}} \\ {0.5 \times {{10}^{ - 6}}} \end{array}} \right]

$$     (3.6.10)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{r_{3y}} = 0

$$     (3.6.11) To sum up, the displacement in X-direction of B is $$ {{10}^{-6}}$$ mm, in Y-direction is $$ {{10}^{-6}}$$ mm. The displacement in X-direction of C is $$ {{10}^{-6}}$$ mm.
 * <p style="text-align:right">
 * }

c) Element 1:
 * {| style="width:100%" border="0"

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

\varphi = 135^\circ ,\cos \varphi  =  - \frac{1},\sin \varphi  = \frac{1}

$$     (3.6.12)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma _1} = \frac{E}\left[ {\begin{array}{ccccccccccccccc} {\frac{1}}&{ - \frac{1}}&{ - \frac{1}}&{\frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {3.328 \times {{10}^{ - 6}}} \\ {0.5 \times {{10}^{ - 6}}} \\ 0 \\  0 \end{array}} \right] =  1.414 \times {10^5}Pa

$$     (3.6.13) Element 2:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\varphi = 90^\circ ,\cos \varphi  = 0,\sin \varphi  = 1

$$     (3.6.14)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma _2} = \frac{E}{1}\left[ {\begin{array}{ccccccccccccccc} 0&{ - 1}&0&1 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {3.328 \times {{10}^{ - 6}}} \\ {0.5 \times {{10}^{ - 6}}} \\ 0 \\  0 \end{array}} \right] = -  0.5 \times {10^5}Pa

$$     (3.6.15) Element 3:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\varphi = 45^\circ ,\cos \varphi  = \frac{1},\sin \varphi  = \frac{1}

$$     (3.6.16)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\sigma _3} = \frac{E}\left[ {\begin{array}{ccccccccccccccc} { - \frac{1}}&{ - \frac{1}}&{\frac{1}}&{\frac{1}} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {3.328 \times {{10}^{ - 6}}} \\ {0.5 \times {{10}^{ - 6}}} \\ {3.828 \times {{10}^{ - 6}}} \\ 0 \end{array}} \right] = 0Pa

$$     (3.6.17)
 * <p style="text-align:right">
 * }

d) Considering the Global system of equations, then:
 * {| style="width:100%" border="0"

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

\frac{1}\left[ {\begin{array}{ccccccccccccccc} 0&0&{ - \frac{1}}&{\frac{1}} \\ 0&0&{\frac{1}}&{ - \frac{1}} \\ 0&0&0&0 \\  0&0&0&{ - 2} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} \\  0 \\    \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right]

$$     (3.6.18) So
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\left[ {\begin{array}{ccccccccccccccc} \\   \\    \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} { - 1 \times {{10}^3}} \\ {1 \times {{10}^3}} \\ 0 \\  { - 1 \times {{10}^3}} \end{array}} \right]

$$     (3.6.19)
 * <p style="text-align:right">
 * }


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

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

{F_{Dx}} = {F_{Dy}} = 0

$$     (3.6.20)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{F_{Cx}} = {r_{2x}} = 0,{F_{Cy}} = {r_{2y}} = - 1 \times {10^3}N

$$     (3.6.21)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{F_{Fx}} = {r_{1x}} - 1 \times {10^3},{F_{Fy}} = {r_{1y}} = 1 \times {10^3}N

$$     (3.6.22)
 * <p style="text-align:right">
 * }

Author
Zongyi Yang

Problem Description
Repeat problem 3.1 with $$k=0$$

Solution
For the case when $$\displaystyle n=2$$

Parts 1 & 2

 * $$\left\{ b_j(x);j=0,1,...,n \right\}$$ where $$\displaystyle b_j(x) = (x+k)^j$$


 * $$\Rightarrow \quad b_0(x)=(x+k)^0=1, b_1(x)=(x+k)^1=(x+k),\ and\ b_2(x)=(x+k)^2$$.

The natural boundary condition can be implemented by differentiating $$\displaystyle u^h(x)$$ with respect to $$\displaystyle x$$ and taking its value at x=0.


 * $$\frac{d{u}^{h}}{dx}(x) = \sum_{j=0}^2{d_jb_j'(x)}=d_0b_0'(x)+d_1b_1'(x)+d_2b_2'(x)=d_1+2d_2(x)$$
 * $$\frac{d{u}^{h}}{dx}(0) = \sum_{j=0}^2{d_jb_j'(0)}=d_0b_0'(0)+d_1b_1'(0)+d_2b_2'(0)=d_1=-4$$

So the relevant equation is


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

$$ \displaystyle d_1=-4$$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

The essential boundary condition, $$\displaystyle u^h(1)=0$$, is implemented as follows:


 * $$u^h(x)= \sum_{j=0}^2{d_jb_j(x)}=d_0b_0(x)+d_1b_1(x)+d_2b_2(x)=d_0+d_1(x)+d_2(x)^2$$
 * $$u^h(1)= \sum_{j=0}^2{d_jb_j(1)}=d_0b_0(1)+d_1b_1(1)+d_2b_2(1)=d_0+d_1+d_2=0$$

So the relevant equation is


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

$$ \displaystyle d_0+d_1+d_2=0$$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

Part 3
Equation 3.1.2 can be used to project the residue onto the basis function. First, the partial differential equation $$\displaystyle P(u^h)$$ is defined.


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

$$ \displaystyle P(u^h) = \frac{\partial }{\partial x}\left[a_2(x)\frac{\partial u}{\partial x}\right]+f(x,t) - \overline {m} (x) \frac{\partial ^s u}{\partial t^s}$$
 * style="width:10%; |
 * style="width:10%; |
 * }
 * }

Substituting the given conditions into this equation and recalling equation 3.1.1,


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

$$ \displaystyle P(u^h) = \frac{\partial }{\partial x}\left\{ 2 \left[d_1+2d_2(x)\right] \right\}+3$$
 * style="width:10%; |
 * style="width:10%; |
 * }
 * }

Differentiating with respect to $$\displaystyle x$$ and substituting known values,


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

$$\displaystyle P(u^h) = \left[4d_2\right]+3$$
 * style="width:10%; |
 * style="width:10%; |
 * }
 * }

Substituting into 3.1.2


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

$$ \displaystyle <b_k,P(u^h)> = \int_0^1 b_k(x) P(u^h)dx = \int_0^1 b_k(x) \left[4d_2 +3\right]dx=0$$
 * style="width:10%;|
 * style="width:10%;|
 * }
 * }

Since the equation is valid for all $$\displaystyle b_k(x), \ b_1(x)=(x+1)$$ is utilized.


 * $$\displaystyle <b_k,P(u^h)> = \int_0^1 (x+1) \left[4d_2 +3\right]dx= \int_0^1 [(3+4d_2)x+4d_2+3]dx=(3+4d_2)\frac{1^2}{2}+(4d_2+3)1 = 6d_2+\frac{9}{2}=0$$


 * $$\displaystyle $$

Simplifying, the relevant equation becomes


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

$$ \displaystyle 4d_2=-3$$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }

Part 4
The coefficient ($$K$$) and constant ($$f$$) matrices are constructed from the equations determined aboved.


 * Notes:
 * The first row of the matrix $$\mathbf K$$ is based on equation 3.1.1 and the essential boundary condition, 3.1.4, when $$\displaystyle x=1$$.


 * The second row of the matrix $$\mathbf K$$ is from equations 3.1.5 and 3.1.6 (for the general case).


 * The third row of the matrix $$\mathbf K$$ is from equation 3.1.11.


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

$$ \begin{bmatrix} 0 & 1 & 0\\ 1& 1 & 1\\  0& 0 & 4 \end{bmatrix} \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |

Note $$\displaystyle\mathbf K$$is not symmetric
 * }
 * }

Part 5
To solve for the matrix $$\mathbf d$$, first recognize that $$\mathbf K \cdot \mathbf d = \mathbf F$$ can be rewritten as $$\mathbf d = \mathbf K^{-1} \cdot \mathbf F$$.


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

$$ \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0\\ 1& 1 & 1\\  0& 0 & 4 \end{bmatrix}^{-1} \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix}$$
 * style="width:10%; |
 * style="width:10%; |

$$ \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} -1 & 1 & -0.25\\ 1& 0 & 0\\  0& 0 & 0.25 \end{bmatrix} \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix} $$


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

$$ \displaystyle \underline {d}_{n=2} = \begin{bmatrix} 4.75\\ -4\\ -0.75 \end{bmatrix} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }
 * }
 * }

Parts 6 & 7
The approximate function, $$\displaystyle u_n^h(x)$$, can now be written using the determined coefficients. The exact solution, was previously obtained and is displayed here:



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

The same procedure is repeated for $$n = 4$$ and $$n = 6$$ in equation 3.1.1.

As the procedure is the same shown above, a Matlab script was written to enforce the boundary conditions, perform the inner products and integration, solve the system of equations, and plot the results.

Figure 3.1.1 shows the approximate solutions for $$n=2,4,6$$ plotted against the exact solution in the domain $$ \Omega $$. It is clear that $$n=2$$ provides a good approximation. This is likely due to the fact that the approximate and exact solutions are both second order polynomial functions.



Part 8
The previous Matlab code also calculates the error in the estimate for u or the following equation:


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

$$\displaystyle e_n(x)=u(x)-u^h(x)$$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

This equation was evaluated at .5 and plotted vs n giving the following figure:



Author
Johnathan Whittaker Bullard

Reviewers / Team Members
Whitaker Bullard/Braden Snook

Problem Description

 * [From Fish & Belytschko, Chapter 3, Problem 3.1]

Show that the weak form of


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


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


 * $$\displaystyle u(3) = 0.001$$

is given by
 * $$\int_{1}^{3}\frac{dw}{dx}AE\frac{du}{dx}dx=-0.1(wA)_{x=1}+\int_{1}^{3}2xw\ dx \qquad \forall w \ \text{ with }\ w(3)=0.$$

Solution
The strong form of this problem can be directly compared to the series of equations for the case of one-dimensional stress analysis, and the derivation is therefore adapted from the text [F&B, Section 3.5] First, the governing equation and traction boundary condition are multiplied by an arbitrary (or weight) function, $$ w(x) $$, and integrated from 1 to 3.


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

$$ \int_{1}^{3} w \left[ \frac{d}{dx} \left( AE\frac{du}{dx} \right) + 2x \right] dx=0 \qquad \forall w,$$ (3.8.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$\left(wA \left(E \frac{du}{dx} -0.1 \right ) \right )_{x=1}=0 \qquad \forall w.$$ (3.8.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

In its current form, equation 3.8.1 is twice differentiable and the stiffness matrix is not symmetric. Therefore it will need to be reduced to first derivatives only and hence a symmetric stiffness matrix. It can be rewritten as


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

$$ \int_{1}^{3} w \frac{d}{dx} \left( AE\frac{du}{dx} \right) dx + \int_{1}^{3} w(2x)dx=0 \qquad \forall w.$$ (3.8.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Using integration by parts we obtain


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

$$ \left( wAE \frac{du}{dx} \right )|_{1}^{3}-\int_{1}^{3}\frac{dw}{dx}AE \frac{du}{dx}dx + \int_{1}^{3} w(2x)dx = 0 \qquad \forall w \ \text{with} \ w(3)=0.$$ (3.8.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

It is important to note that in the derivation of the weak form, the weight functions and trial solutions are constructed so that $$w(3)=0 \text{ and } u=\bar{u} \text{ on } \Gamma_u$$, respectively. $$\displaystyle\Gamma_u$$ is defined as the boundary where the displacements are prescribed (x=3 in this case); $$\displaystyle\Gamma_t$$ is the boundary where the traction is prescribed (x=1).

Recalling Hooke's Law and the definition of strain equation 3.8.4 is rewritten as


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

$$ (wA \sigma )_{x=3} - (wA \sigma )_{x=1} -\int_{1}^{3}\frac{dw}{dx}AE \frac{du}{dx}dx + \int_{1}^{3} w(2x)dx = 0 \qquad \forall w \ \text{with} \ w(3)=0.$$ (3.8.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The original traction boundary can be applied to the second term. Recalling that $$\displaystyle w(3)=0$$ removes the first term. Thus we are left with the weak form in the case of one-dimensional stress analysis.


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

$$ \int_{1}^{3}\frac{dw}{dx}AE \frac{du}{dx}dx = -0.1(wA)_{x=1} + \int_{1}^{3} w(2x)dx = 0 \qquad \forall w \ \text{with} \ w(3)=0.$$ (3.8.6)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Author
Philip Flater

Problem Description
"A First Course in Finite Elements" by Fish and Belytschko p. 72, pb.3.3 Consider a trial (candidate) solution of the form $$ \displaystyle u(x)=\alpha _{0}+\alpha _{1}(x-3)$$ and a weight function of the same form. Obtain a solution to the weak form in pb.3.1 from Fish and Belytschko. Check the equilibrium equation in the strong from in pb.3.1 from Fish and Belytschko; is it satisfied? Check the natural boundary condition; is it satisfied?

Solution
The candidate solution:
 * {| style="width:100%" border="0"

$$ \displaystyle u(x)=\alpha _{0}+\alpha _{1}(x-3) $$       (3.9.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

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

$$ \displaystyle w(x)=\beta_{0}+\beta_{1}(x-3) $$       (3.9.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

From pb.3.1 from Fish and Belytschko (eqn 3.8.6 from this hw submission):


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

$$ \displaystyle \int_{1}^{3}\frac{dw}{dx}AE \frac{du}{dx}dx = -0.1(wA)_{x=1} + \int_{1}^{3} w(2x)dx = 0 \qquad \forall w \ \text{with} \ w(3)=0$$ (3.9.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The weight function $$ \displaystyle w(3)=0 $$

$$ \displaystyle w(3)=0=\beta_{0}+\beta _{1}(3-3) $$


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

$$ \displaystyle \beta_{0}=0 $$       (3.9.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

From pb.3.1 from Fish and Belytschko:


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

$$ \displaystyle \alpha_{0}=.001 $$       (3.9.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle u(3)=\alpha_{0}+\alpha_{1}(3-3)=.001 $$       (3.9.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \alpha_{0}=.001 $$       (3.9.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The candidate solution and the weight function is now:


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

$$ \displaystyle u(x)=.001+\alpha_{1}(x-3) $$       (3.9.7)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle w(x)=\beta_{1}(x-3) $$       (3.9.8)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Recall that the weak form is
 * {| style="width:100%" border="0"

$$ \displaystyle \int_{1}^{3}\frac{dw}{dx}AE\frac{du}{dx}dx=-0.1(wA)_{x=1}+\int_{1}^{3}2xw\ dx \qquad \forall w \ \text{ with }\ w(3)=0 $$       (3.9.9)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

From (3.9.7) and (3.9.8):
 * {| style="width:100%" border="0"

$$ \displaystyle \frac{du(x)}{dx}=\alpha_{1} $$       (3.9.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \frac{dw(x)}{dx}=\beta_{1} $$       (3.9.11)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Plug (3.9.7), (3.9.8), (3.9.10) and (3.9.11) into (3.9.9):


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

$$ \displaystyle \int_{1}^{3}\beta_{1}AE\alpha_{1}dx=-0.1(\beta_{1}(x-3)A)_{x=1}+\int_{1}^{3}2x\beta_{1}(x-3)\ dx \qquad \forall w \ \text{ with }\ w(3)=0 $$       (3.9.12)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Integrage (3.9.12)


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

$$ \displaystyle 2\beta_{1}AE\alpha_{1}=-0.1(\beta_{1}(1-3)A)-2\frac{10}{3}\beta_{1} $$       (3.9.13)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Bring everything to one side and factor out $$ \displaystyle \beta_{1} $$


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

$$ \displaystyle \beta_{1}(2AE\alpha_{1}-0.2A+\frac{20}{3}) =0 $$       (3.9.14)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \beta_{1}=0 $$       (3.9.15)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle 2AE\alpha_{1}-0.2A+\frac{20}{3} =0 $$       (3.9.16)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle \alpha_{1}=\left (0.2A-\frac{20}{3}  \right )\left ( \frac{1}{2AE} \right ) $$       (3.9.17)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Plug (3.9.17) into (3.9.7) and the candidate solution becomes:
 * {| style="width:100%" border="0"


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |


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

$$ \displaystyle u(x)=.001+\left (\frac{0.1}{E}-\frac{10}{3AE} \right )(x-3) $$       (3.9.18)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * style = |
 * }
 * }

From pb.3.1 from Fish and Belytschko, the equilibrium equation is:


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

$$ \displaystyle \frac{d}{dx}\left ( AE\frac{du}{dx} \right )+2x=0$$ (3.9.19)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Substitute (3.9.18) into (3.9.19)


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

$$ \displaystyle \frac{d}{dx}\left ( AE\left (\frac{0.1}{E}-\frac{10}{3AE} \right) \right )+2x=0 $$       (3.9.20)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |


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

$$ \displaystyle 0+2x \neq 0 $$       (3.9.21)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The equilibrium equation is not satisfied.
 * style = |
 * }
 * }

From pb.3.1 from Fish and Belytschko, the boundary condition is:


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

$$ \displaystyle \sigma (1)=\left ( E\frac{du}{dx} \right )_{x=1}=0.1 $$       (3.9.22)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Substitute (3.9.18) into (3.9.22):
 * {| style="width:100%" border="0"

$$ \displaystyle \sigma (1)=\left ( E \left (\frac{0.1}{E}-\frac{10}{3AE} \right) \right )_{x=1}=0.1 $$       (3.9.22)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |


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

$$ \displaystyle {0.1}-\frac{10}{3A} =0.1 $$       (3.9.23) Which is true when $$ \displaystyle  A\gg\frac{10}{3} $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * style = |
 * }
 * }

Author
Brandonhua 04:48, 15 February 2011 (UTC)

Problem Description
Adapted from "A First Course in Finite Elements" by Fish and Belytschko p. 72 pb.3.4

Consider a trial solution of the form $$ u\left( x \right) = {\alpha _0} + {\alpha _1}\left( {x - 3} \right) + {\alpha _2}{\left( {x - 3} \right)^2}$$ and a weight function of the same form. Obtain a solution to the weak form in Problem 3.1 in Fish and Belytschko; is it satisfied?

Check the boundary condition; is it satisfied?

The weak form is given by
 * {| style="width:100%" border="0"

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

\begin{array}{ccccccccccccccc} {\int\limits_1^3 {\fracAE\fracdx = - 0.1{{\left( {wA} \right)}_{x = 1}} + \int\limits_1^3 {2xwdx} } }&{\forall w{\text{ with }}w\left( 3 \right) = 0} \end{array}

$$     (Eq 10.1)
 * <p style="text-align:right">
 * }

Solution
Trial solution and weight function are in the form of
 * {| style="width:100%" border="0"

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

u\left( x \right) = {\alpha _0} + {\alpha _1}\left( {x - 3} \right) + {\alpha _2}{\left( {x - 3} \right)^2}

$$     (Eq 10.2)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

w\left( x \right) = {\beta _0} + {\beta _1}\left( {x - 3} \right) + {\beta _2}{\left( {x - 3} \right)^2}

$$     (Eq 10.3) Where $${\alpha _0}$$, $${\alpha _1}$$ and $${\alpha _2}$$ are unknown parameters and $${\beta _0}$$, $${\beta _1}$$ and $${\beta _2} $$are arbitrary parameters.
 * <p style="text-align:right">
 * }

Assume that the area A is a constant.

To be admissible the weight function must vanish at $$x = 3$$, so $${\beta _0} = 0$$, and the trial solution must satisfy the essential boundary condition $$u\left( 3 \right) = 0.001$$, so $${\alpha _0} = 0.001$$.

Then, we can get
 * {| style="width:100%" border="0"

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

\frac = {\alpha _1} + 2{\alpha _2}\left( {x - 3} \right)

$$     (Eq 10.4)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\frac = {\beta _1} + 2{\beta _2}\left( {x - 3} \right)

$$     (Eq 10.5)
 * <p style="text-align:right">
 * }

Substituting the above into the weak form gives
 * {| style="width:100%" border="0"

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

\int\limits_1^3 {\left( {{\beta _1} + 2{\beta _2}\left( {x - 3} \right)} \right)\left( {{\alpha _1} + 2{\alpha _2}\left( {x - 3} \right)} \right)AEdx = - 0.1A{{\left( {{\beta _1}\left( {x - 3} \right) + {\beta _2}{{\left( {x - 3} \right)}^2}} \right)}_{x = 1}} + \int\limits_1^3 {2x\left( {{\beta _1}\left( {x - 3} \right) + {\beta _2}{{\left( {x - 3} \right)}^2}} \right)dx} }

$$     (Eq 10.6) Integrating and rearranging the terms gives
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\beta _1}\left( {E\left( {{\text{2}}{\alpha _{\text{1}}}A - 4{\alpha _2}A} \right) - \left( {0.2A - \frac{3}} \right)} \right) + {\beta _2}\left( {E\left( { - 4{\alpha _1}A + \frac{3}{\alpha _2}A} \right) - \left( {8 - 0.4A} \right)} \right) = 0

$$     (Eq 10.7)
 * <p style="text-align:right">
 * }

As the above must hold for arbitrary weight functions, it must for arbitrary $${\beta _1}$$ and $${\beta _2}$$, which gives the following linear algebraic equation:
 * {| style="width:100%" border="0"

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

E\left[ {\begin{array}{ccccccccccccccc} {2A}&{ - 4A} \\ { - 4A}&{\frac{3}A} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} \\ \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} {0.2A - \frac{3}} \\ {8 - 0.4A} \end{array}} \right]

$$     (Eq 10.8)
 * <p style="text-align:right">
 * }

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

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

{\alpha _1} = - \frac{1}{E}\left( {\frac + 0.1} \right)

$$     (Eq 10.9)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

{\alpha _2} = - \frac{2}

$$     (Eq 10.10)
 * <p style="text-align:right">
 * }

Thus, the trial solution is:
 * {| style="width:100%" border="0"

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

u\left( x \right) = 0.001 - \frac{1}{E}\left( {\frac - 0.1} \right)\left( {x - 3} \right) - \frac{2}{\left( {x - 3} \right)^2}

$$     (Eq 10.11)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

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

\sigma \left( x \right) = - \frac{4}{A}x + \frac + 0.1

$$     (Eq 10.12)
 * <p style="text-align:right">
 * }

Substituting the trial solution into the equilibrium equation in the strong form gives
 * {| style="width:100%" border="0"

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

- 4 + 2x = 0

$$     (Eq 10.13)
 * <p style="text-align:right">
 * }

Only at $$x = 2$$ satisfied.

Substituting the trial solution into the natural boundary condition gives
 * {| style="width:100%" border="0"

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

{\sigma _1}\left( 1 \right) = {\left( {E\frac} \right)_{x = 1}} = 0.1 + \frac{2} = 0.1

$$     (Eq 10.14)
 * <p style="text-align:right">
 * }

Only $$A \to \infty$$ Satisfied.

Author
Jiang Jin

Problem Description

 * Adapted from 3.1 Problem Description.


 * Consider basis functions $$\left\{ b_j(x);j=0,1,...,n \right\}$$ where

$$\displaystyle \mathbf b_j(x)=\left\{1,cos(x),sin(x),cos(2x),sin(2x).....cos(mx)_{j-1},sin(mx)_j\right\};$$

$$\displaystyle m=1.....n/2$$

1. Let n=2 $$\rightarrow$$ ndof = n+1 = 2+1 = 3 with basis coefficients given by $$\mathbf {d} = \left\{d_j;j=0,...,n\right\}_{(n+1)\times1}$$.

2. Find 2 equations that enforce the boundary conditions for
 * {| style="width:100%" border="0"

$$u^h(x) = \sum_{j=0}^n{d_jb_j(x)}$$. (3.11.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

3. Find 1 more equation (i.e. j = 0,1,2) to solve for $$\mathbf {d} = \left\{d_j\right\}_{3\times1}$$ by projecting the residue, $$\displaystyle P(u^h)$$, on a basis function, $$\displaystyle b_k(x)$$, with k = 0,1,2 such that the additional equation is linearly independent from the equations in part 2.


 * Notes regarding residue:
 * The general form of the residue is defined as $$\displaystyle P(u^h)$$. The exact solution is expressed as $$\displaystyle P(u) = 0$$.


 * $$\displaystyle u^h$$ can approximate $$\displaystyle u$$. In other words $$\displaystyle u^h \cong u $$. Due to this approximation for $$\displaystyle u$$, $$ P(u^h) \ne 0 \quad \forall x \in \Omega $$ (for all values of x in the domain).
 * Notes regarding projection:
 * $$\displaystyle\int_{\Omega } b_i(x) P\left(u^h(x)\right)dx = 0 \quad\Rightarrow\quad\int b_k(x) P(u^h)dx = 0$$
 * $$\displaystyle\mathbf b_i \cdot \mathbf P(\underline v) = 0 \quad\Rightarrow\quad \mathbf b_i \cdot \mathbf P(\underline v) \ = \  <\underline b_i,\underline P(\underline v)>$$


 * Combining the previous two equations yields
 * {| style="width:100%" border="0"

$$\displaystyle <b_k,P(u^h)> = \int b_k(x) P(u^h)dx$$. (3.11.2) 4. Display 3 equations in matrix form, $$\underline K \cdot \underline d = \underline F$$.
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

5. Solve for $$\mathbf d$$.

6. Construct $$\displaystyle u_n^h(x)$$ and plot versus the exact solution, $$\displaystyle u(x)$$.

7. Repeat steps 1 through 6 for:
 * n = 4
 * n = 6


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

$$ \displaystyle \mathbf b_i(x)=\left\{1,cos(x),sin(x),cos(2x),sin(2x).....cos(mx),sin(mx)\right\};$$
 * style="width:95%" |
 * style="width:95%" |

$$\displaystyle i=0,1,....n,m=1.....n/2$$ (3.11.3)
 * <p style="text-align:right">
 * }

Given

 * See 3.1 Given.

Solution
Following a similar method as on 3.1 Solution

For the case when $$\displaystyle n=2$$

Boundary Condition Equations

 * $$ \quad b_0(x)=1, b_1(x)=cos(x),\ and\ b_2(x)=sin(x)$$.

The natural boundary condition (Eq. 3.1.3 in Given) can be implemented by differentiating $$\displaystyle u^h(x)$$ with respect to $$\displaystyle x$$ and taking its value at x=0.


 * $$\frac{d{u}^{h}}{dx}(x) = \sum_{j=0}^2{d_jb_j'(x)}=d_0b_0'(x)+d_1b_1'(x)+d_2b_2'(x)=$$

$$-d_1sin(x)+d_2cos(x)$$
 * $$\frac{d{u}^{h}}{dx}(0) = \sum_{j=0}^2{d_jb_j'(0)}=-d_1sin(0)+d_2cos(0)$$

So the relevant equation is


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

$$ \displaystyle d_2=-4$$ (3.11.4)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

The essential boundary condition, $$\displaystyle u^h(1)=0$$, is implemented as follows:


 * $$u^h(x)= \sum_{j=0}^2{d_jb_j(x)}=d_0b_0(x)+d_1b_1(x)+d_2b_2(x)=$$

$$d_0+d_1cos(x)+d_2sin(x)$$
 * $$u^h(1)= \sum_{j=0}^2{d_jb_j(1)}=d_0+d_1cos(1)+d_2sin(1)=0$$

So the relevant equation is


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

$$ \displaystyle d_0+d_1cos(1)+d_2sin(1)=0$$ (3.11.5)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Part 3
Equation 3.1.2 can be used to project the residue onto the basis function. First, the partial differential equation $$\displaystyle P(u^h)$$ is defined.


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

$$ \displaystyle P(u^h) = \frac{\partial }{\partial x}\left[a_2(x)\frac{\partial u}{\partial x}\right]+f(x,t) - \overline {m} (x) \frac{\partial ^s u}{\partial t^s}$$ (3.11.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Substituting the given conditions into this equation and recalling equation 3.1.1,


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

$$ \displaystyle P(u^h) = \frac{\partial }{\partial x}\left\{ 2 \left[-d_1sin(x)+d_2cos(x)\right] \right\}+3$$ (3.11.7)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Differentiating with respect to $$\displaystyle x$$:


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

$$\displaystyle P(u^h) = \left[-2d_1cos(x)-2d_2sin(x)\right]+3$$ (3.11.8)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Substituting into 3.11.2 and using given domain


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

$$ \displaystyle <b_k,P(u^h)> = \int_0^1 b_k(x) P(u^h)dx = \int_0^1 b_k(x) \bigg [\left[-2d_1cos(x)-2d_2sin(x)\right]+3\bigg ] dx=0$$ (3.11.9)
 * style="width:10%;|
 * style="width:10%;|
 * <p style="text-align:right">
 * }

Since the equation is valid for all $$\displaystyle b_k(x), \ b_0(x)=1$$ is utilized.


 * $$\displaystyle <b_k,P(u^h)> = \int_0^1 b_k(x) \big [\left[-2d_1cos(x)-2d_2sin(x)\right]+3\big ]=-2d_1sin(x)+2d_2cos(x)+3x|_0^1=0$$


 * $$\displaystyle $$

Simplifying, the relevant equation becomes


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

$$ \displaystyle -2d_1sin(1)+d_2(cos(1)-2)+3=0$$ (3.11.10)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Part 4
The coefficient ($$K$$) and constant ($$f$$) matrices are constructed from the equations determined above.


 * Notes:
 * The first row of the matrix $$\mathbf K$$ is based on equation 3.11.5 from the essential boundary condition.


 * The second row of the matrix $$\mathbf K$$ is from equations 3.11.4 from the natural boundary condition.


 * The third row of the matrix $$\mathbf K$$ is from equation 3.11.10 projecting 0 residue on the basis functions (minimizing error).


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

$$ \begin{bmatrix} 1 & cos(1) & sin(1)\\ 0 & 0 & 1\\ -2sin(1)& cos(1)-3 & 0 \end{bmatrix} \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} 0\\ -4\\ -3 \end{bmatrix} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |

Note $$\displaystyle\mathbf K$$ is not symmetric. (3.11.11)
 * <p style="text-align:right">
 * }

Part 5
To solve for the matrix $$\mathbf d$$, first recognize that $$\mathbf K \cdot \mathbf d = \mathbf F$$ can be rewritten as $$\mathbf d = \mathbf K^{-1} \cdot \mathbf F$$.


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

$$ \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} 1 & cos(1) & sin(1)\\ 0 & 2 & 0\\ -2sin(1)& cos(1)-3 & 0 \end{bmatrix}^{-1} \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix}$$
 * style="width:10%; |
 * style="width:10%; |

$$ \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} 1 & -.5463 & .32105\\ 0& -.54630 & -0.5942\\  0& 1 & 0 \end{bmatrix} \begin{bmatrix} -4\\ 0\\ -3 \end{bmatrix} $$


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

$$ \displaystyle \underline {d}_{n=2} = \begin{bmatrix} 1.2221\\ 3.9678\\ -4 \end{bmatrix} $$     (3.11.12)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }
 * }

Parts 6 & 7
The approximate function, $$\displaystyle u_n^h(x)$$, can now be written using the determined coefficients in equation 3.1.13. The exact solution, was obtained and is displayed here:



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

The same procedure is repeated for $$n = 4$$ and $$n = 6$$ in equation 3.1.1.

As the procedure is the same shown above, a Matlab script was written to enforce the boundary conditions, perform the inner products and integration, solve the system of equations, and plot the results.

In order to code this in Matlab, the stiffness matrix had to be defined as follows. If the stiffness matrix row was even (consequently multiplying by cos basis functions) the following was used (using Wolframalpha, integral of cos(nx)*cos(mx) and m equal to n, integral of cos(nx)*cos(mx) n not equal to m), integral of cos(nx)sin(mx) with n equal to m, and integral of cos(nx)sin(mx) with n not equal to m.


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

$$ \displaystyle K_{i,j}=2*-m^2\int cos(nx)cos(mx)dx $$
 * style="width:10%; |
 * style="width:10%; |

$$ = 2*-m^2\left\{ \begin{array}{l l}   \frac{2nx+sin(2nx)}{4n} & \quad \text{if }n=m\\ \frac{msin(mx)cos(nx)-ncos(mx)sin(nx)}{m^2-n^2} & \quad \text{if }n\neq m\\ \end{array} \right.$$ (3.11.13)
 * <p style="text-align:right">
 * }


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

$$ \displaystyle K_{i,j+1}=2*-m^2\int cos(nx)sin(mx) $$
 * style="width:10%; |
 * style="width:10%; |

$$ =2*-m^2 \left\{ \begin{array}{l l}   -\frac{cos(nx)^2}{2n} & \quad \text{if }n=m\\ -\frac{nsin(mx)sin(nx)+mcos(mx)cos(nx)}{m^2-n^2} & \quad \text{if }n\neq m\\ \end{array} \right.$$ (3.11.14)
 * <p style="text-align:right">
 * }

When on an odd row of the stiffness matrix the basis function used for the inner product was sin. Once again using WolframAlpha for a symbolic solution (integral of sin(nx)*cos(mx) n equal to m, integral of sin(nx)cos(mx) n not equal to m, integral of sin(nx)sin(mx) n equal to m, and integral of sin(nx)sin(mx) n not equal to m.


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

$$ \displaystyle K_{i,j}=2*-m^2\int sin(nx)cos(mx) $$
 * style="width:10%; |
 * style="width:10%; |

$$ =2*-m^2 \left\{ \begin{array}{l l}   -\frac{cos(nx)^2}{2n} & \quad \text{if }n=m\\ -\frac{nsin(mx)sin(nx)+mcos(mx)cos(nx)}{m^2-n^2} & \quad \text{if }n\neq m\\ \end{array} \right.$$ (3.11.15)
 * <p style="text-align:right">
 * }


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

$$ \displaystyle K_{i,j+1}=2*-m^2\int sin(nx)sin(mx) $$
 * style="width:10%; |
 * style="width:10%; |

$$ =2*-m^2 \left\{ \begin{array}{l l}   \frac{x}{2}-\frac{sin(2nx)}{4n} & \quad \text{if }n=m\\ \frac{nsin(mx)cos(nx)-mcos(mx)sin(nx)}{m^2-n^2} & \quad \text{if }n\neq m\\ \end{array} \right.$$ (3.11.16)
 * <p style="text-align:right">
 * }

This was implemented in the following Matlab code:

Figure 3.11.1 shows the approximate solutions for $$n=2,4,6$$ plotted against the exact solution in the domain $$ \Omega $$. It appears that n=4 and n=6 are both better approximations than n=2.



Part 8
The previous Matlab code also calculates the error in the estimate for u or the following equation:


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

$$\displaystyle e_n(x)=u(x)-u^h(x)$$ (3.11.17)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

This equation was evaluated at .5 and plotted vs n giving the following figure:



Author
Braden Snook