# Nonlinear finite elements/Axial bar finite element solution
## Axially loaded bar: The Finite Element Solution

The finite element method is a type of Galerkin method that has the following advantages:
* The functions $\varphi _{i}$ are found in a systematic manner.
* The functions $\varphi _{i}$ are chosen such that they can be used for arbitrary domains.
* The functions $\varphi _{i}$ are piecewise polynomials.
* The functions $\varphi _{i}$ are non-zero only on a small part of the domain.
As a result, computations can be done in a modular manner that is suitable for computer implementation.

### Discretization

The first step in the finite element approach is to divide the domain into elements and nodes, i.e., to create the finite element mesh.
Let us consider a simple situation and divide the rod into 3 elements and 4 nodes as shown in Figure 6.
* Figure 6. Finite element mesh and basis functions for the bar.

### Shape functions

The functions $\varphi _{i}$ have special characteristics in finite element methods and are generally written as $N_{i}$ and are called basis functions, shape functions, or interpolation functions.
Therefore, we may write
    ${\begin{aligned}K_{ij}&=\int _{0}^{L}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\\f_{j}&=\int _{0}^{L}N_{j}\mathbf {q} ~dx+\left.N_{j}{\boldsymbol {R}}\right|_{x=L}~.\end{aligned}}$
The finite element basis functions are chosen such that they have the following properties:
* The functions $N_{i}$ are bounded and continuous.
* If there are $n$ nodes, then there are $n$ basis functions - one for each node. There are four basis functions for the mesh shown in Figure 6.
* Each function $N_{i}$ is nonzero only on elements connected to node $i$.
* $N_{i}$ is 1 at node $i$ and zero at all other nodes.

### Stiffness matrix

Let us compute the values of $K_{ij}$ for the three element mesh. We have
    ${\begin{aligned}K_{ij}&=\int _{0}^{L}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\\&=\int _{0}^{x_{2}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{x_{2}}^{x_{3}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{x_{3}}^{L}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\\&\equiv \int _{\Omega _{1}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\end{aligned}}$
The components of $\mathbf {K}$ are
    ${\begin{aligned}K_{11}&=\int _{\Omega _{1}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{12}&=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{13}&=\int _{\Omega _{1}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{14}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{22}&=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx\\K_{23}&=\int _{\Omega _{1}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx\\K_{24}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx\\K_{33}&=\int _{\Omega _{1}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx\\K_{34}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx\\K_{44}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx~.\end{aligned}}$
The matrix $\mathbf {K}$ is symmetric, so we don't need to explicitly compute the other terms.
From Figure 6, we see that $N_{1}$ is zero in elements 3 and 4, $N_{2}$ is zero in element 4, $N_{3}$ is zero in element 1, and $N_{4}$ is zero in elements 1 and 2. The same holds for $dN_{i}/dx$.
Therefore, the coefficients of the $\mathbf {K}$ matrix become
    ${\begin{aligned}K_{11}&=\int _{\Omega _{1}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx~;~K_{12}=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx~;~K_{13}=0~;~K_{14}=0\\K_{22}&=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx~;~K_{23}=\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx~;~K_{24}=0\\K_{33}&=\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx~;~K_{34}=\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx\\K_{44}&=\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx~.\end{aligned}}$
We can simplify our calculation further by letting $N_{k}^{e}$ be the shape functions over an element $e$. For example, the shape functions over element $2$ are $N_{1}^{2}$ and $N_{2}^{2}$ where the local nodes $1$ and $2$ correspond to global nodes $2$ and $3$, respectively. Then we can write,
    ${\begin{aligned}K_{11}&=\int _{\Omega _{1}}{\cfrac {dN_{1}^{1}}{dx}}AE{\cfrac {dN_{1}^{1}}{dx}}~dx~;~K_{12}=\int _{\Omega _{1}}{\cfrac {dN_{2}^{1}}{dx}}AE{\cfrac {dN_{1}^{1}}{dx}}~dx~;~K_{13}=0~;~K_{14}=0\\K_{22}&=\int _{\Omega _{1}}{\cfrac {dN_{2}^{1}}{dx}}AE{\cfrac {dN_{2}^{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{1}^{2}}{dx}}AE{\cfrac {dN_{1}^{2}}{dx}}~dx~;~K_{23}=\int _{\Omega _{2}}{\cfrac {dN_{2}^{2}}{dx}}AE{\cfrac {dN_{1}^{2}}{dx}}~dx~;~K_{24}=0\\K_{33}&=\int _{\Omega _{2}}{\cfrac {dN_{2}^{2}}{dx}}AE{\cfrac {dN_{2}^{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{1}^{3}}{dx}}AE{\cfrac {dN_{1}^{3}}{dx}}~dx~;~K_{34}=\int _{\Omega _{3}}{\cfrac {dN_{2}^{3}}{dx}}AE{\cfrac {dN_{1}^{3}}{dx}}~dx\\K_{44}&=\int _{\Omega _{3}}{\cfrac {dN_{2}^{3}}{dx}}AE{\cfrac {dN_{2}^{3}}{dx}}~dx~.\end{aligned}}$
Let $K_{kl}^{e}$ be the part of the value of $K_{ij}$ that is contributed by element $e$. The indices $kl$ are local and the indices $ij$ are global. Then,
    ${\begin{aligned}K_{11}&=K_{11}^{1}~;~K_{12}=K_{12}^{1}~;~K_{13}=0~;~K_{14}=0\\K_{22}&=K_{22}^{1}+K_{11}^{2}~;~~K_{23}=K_{12}^{2}~;~~K_{24}=0\\K_{33}&=K_{22}^{2}+K_{11}^{3}~;~K_{34}=K_{12}^{3}\\K_{44}&=K_{22}^{3}\\\end{aligned}}$
We can therefore see that if we compute the stiffness matrices over each element and assemble them in an appropriate manner, we can get the global stiffness matrix $\mathbf {K}$.

#### Stiffness matrix for two-noded elements

For our problem, if we consider an element $e$ with two nodes, the local hat shape functions have the form
    $N_{1}^{e}(\mathbf {x} )={\cfrac {\mathbf {x} _{2}-\mathbf {x} }{h_{e}}}~;~~N_{2}^{e}(\mathbf {x} )={\cfrac {\mathbf {x} -\mathbf {x} _{1}}{h_{e}}}$
where $h_{e}$ is the length of the element.
Then, the components of the element stiffness matrix are
    ${\begin{aligned}K_{11}^{e}&=\int _{x_{1}}^{x_{2}}{\cfrac {dN_{1}^{e}}{dx}}AE{\cfrac {dN_{1}^{e}}{dx}}~dx=\int _{x_{1}}^{x_{2}}\left(-{\cfrac {1}{h_{e}}}\right)AE\left(-{\cfrac {1}{h_{e}}}\right)~dx={\cfrac {AE}{h_{e}}}\\K_{12}^{e}&=\int _{x_{1}}^{x_{2}}{\cfrac {dN_{2}^{e}}{dx}}AE{\cfrac {dN_{1}^{e}}{dx}}~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {1}{h_{e}}}\right)AE\left(-{\cfrac {1}{h_{e}}}\right)~dx=-{\cfrac {AE}{h_{e}}}\\K_{22}^{e}&=\int _{x_{1}}^{x_{2}}{\cfrac {dN_{2}^{e}}{dx}}AE{\cfrac {dN_{2}^{e}}{dx}}~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {1}{h_{e}}}\right)AE\left({\cfrac {1}{h_{e}}}\right)~dx={\cfrac {AE}{h_{e}}}\end{aligned}}$
In matrix form,
    ${\mathbf {K} ^{e}={\cfrac {AE}{h_{e}}}{\begin{bmatrix}1&-1\\-1&1\end{bmatrix}}}$
The components of the global stiffness matrix are
    ${\begin{aligned}K_{11}&={\cfrac {AE}{h_{1}}}~;~~K_{12}=-{\cfrac {AE}{h_{1}}}~;~~K_{13}=0~;~~K_{14}=0\\K_{22}&={\cfrac {AE}{h_{1}}}+{\cfrac {AE}{h_{2}}}~;~~K_{23}=-{\cfrac {AE}{h_{2}}}~;~~K_{24}=0\\K_{33}&={\cfrac {AE}{h_{2}}}+{\cfrac {AE}{h_{3}}}~;~~K_{34}=-{\cfrac {AE}{h_{3}}}\\K_{44}&={\cfrac {AE}{h_{3}}}\\\end{aligned}}$
In matrix form,
    ${\mathbf {K} =AE{\begin{bmatrix}{\cfrac {1}{h_{1}}}&-{\cfrac {1}{h_{1}}}&0&0\\&{\cfrac {1}{h_{1}}}+{\cfrac {1}{h_{2}}}&-{\cfrac {1}{h_{2}}}&0\\&&{\cfrac {1}{h_{2}}}+{\cfrac {1}{h_{3}}}&-{\cfrac {1}{h_{3}}}\\{\text{Symm.}}&&&{\cfrac {1}{h_{3}}}\\\end{bmatrix}}}$

### Load vector

Similarly, for the load vector $\mathbf {f}$, we have
    ${\begin{aligned}f_{j}&=\int _{0}^{L}N_{j}\mathbf {q} ~dx+\left.N_{j}{\boldsymbol {R}}\right|_{x=L}\\&=\int _{\Omega _{1}}N_{j}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{j}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{j}\mathbf {q} ~dx+\left.N_{j}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}$
The components of the load vector are
    ${\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{1}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{1}\mathbf {q} ~dx+\left.N_{1}{\boldsymbol {R}}\right|_{x=L}\\f_{2}&=\int _{\Omega _{1}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{2}\mathbf {q} ~dx+\left.N_{2}{\boldsymbol {R}}\right|_{x=L}\\f_{3}&=\int _{\Omega _{1}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{3}\mathbf {q} ~dx+\left.N_{3}{\boldsymbol {R}}\right|_{x=L}\\f_{4}&=\int _{\Omega _{1}}N_{4}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{4}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{4}\mathbf {q} ~dx+\left.N_{4}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}$
Once again, since $N_{1}$ is zero in elements 2 and 3, $N_{2}$ is zero in element 3, $N_{3}$ is zero in element 1, and $N_{4}$ is zero in elements 1 and 2, we have
    ${\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}\mathbf {q} ~dx+\left.N_{1}{\boldsymbol {R}}\right|_{x=L}\\f_{2}&=\int _{\Omega _{1}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{2}\mathbf {q} ~dx+\left.N_{2}{\boldsymbol {R}}\right|_{x=L}\\f_{3}&=\int _{\Omega _{2}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{3}\mathbf {q} ~dx+\left.N_{3}{\boldsymbol {R}}\right|_{x=L}\\f_{4}&=\int _{\Omega _{3}}N_{4}\mathbf {q} ~dx+\left.N_{4}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}$
Now, the boundary $\mathbf {x} =L$ is at node 4 which is attached to element 3. The only non-zero shape function at this node is $N_{4}$. Therefore, we have
    ${\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}\mathbf {q} ~dx\\f_{2}&=\int _{\Omega _{1}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{2}\mathbf {q} ~dx\\f_{3}&=\int _{\Omega _{2}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{3}\mathbf {q} ~dx\\f_{4}&=\int _{\Omega _{3}}N_{4}\mathbf {q} ~dx+\left.N_{4}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}$
In terms of element shape functions, the above equations can be written as
    ${\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}^{1}\mathbf {q} ~dx=f_{1}^{1}\\f_{2}&=\int _{\Omega _{1}}N_{2}^{1}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{1}^{2}\mathbf {q} ~dx=f_{2}^{1}+f_{1}^{2}\\f_{3}&=\int _{\Omega _{2}}N_{2}^{2}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{1}^{3}\mathbf {q} ~dx=f_{2}^{2}+f_{1}^{3}\\f_{4}&=\int _{\Omega _{3}}N_{2}^{3}\mathbf {q} ~dx+\left.N_{2}^{3}{\boldsymbol {R}}\right|_{x=L}=f_{2}^{3}+{\boldsymbol {R}}\end{aligned}}$
The above shows that the global load vector can also be assembled from the element load vectors if we use finite element shape functions.

#### Load vector for two-noded elements

Using the linear shape functions discussed earlier and replacing $\mathbf {q}$ with $a\mathbf {x}$, the components of the element load vector $\mathbf {f} ^{e}$ are
    ${\begin{aligned}f_{1}^{e}&=\int _{x_{1}}^{x_{2}}N_{1}^{e}a\mathbf {x} ~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {x_{2}-x}{h_{e}}}\right)a\mathbf {x} ~dx={\cfrac {a}{h_{e}}}\left({\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\\f_{2}^{e}&=\int _{x_{1}}^{x_{2}}N_{2}^{e}a\mathbf {x} ~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {x-x_{1}}{h_{e}}}\right)a\mathbf {x} ~dx=-{\cfrac {a}{h_{e}}}\left({\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\end{aligned}}$
In matrix form, the element load vector is written
    ${\mathbf {f} ^{e}={\cfrac {a}{h_{e}}}{\begin{bmatrix}{\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\\{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\end{bmatrix}}}$
Therefore, the components of the global load vector are
    ${\begin{aligned}f_{1}&={\cfrac {a}{h_{1}}}\left({\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\\f_{2}&={\cfrac {a}{h_{1}}}\left({\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\right)+{\cfrac {a}{h_{2}}}\left({\cfrac {x_{3}(x_{3}^{2}-x_{2}^{2})}{2}}-{\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}\right)\\f_{3}&={\cfrac {a}{h_{2}}}\left({\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}-{\cfrac {x_{2}(x_{3}^{2}-x_{2}^{2})}{2}}\right)+{\cfrac {a}{h_{3}}}\left({\cfrac {x_{4}(x_{4}^{2}-x_{3}^{2})}{2}}-{\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}\right)\\f_{4}&={\cfrac {a}{h_{3}}}\left({\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}-{\cfrac {x_{3}(x_{4}^{2}-x_{3}^{2})}{2}}\right)+{\boldsymbol {R}}\end{aligned}}$

### Displacement trial function

Recall that we assumed that the displacement can be written as
    $\mathbf {u} _{h}(\mathbf {x} )=a_{1}\varphi _{1}(\mathbf {x} )+a_{2}\varphi _{2}(\mathbf {x} )+\dots +a_{n}\varphi _{n}(\mathbf {x} )=\sum _{i=1}^{n}a_{i}\varphi _{i}(\mathbf {x} )~.$
If we use finite element shape functions, we can write the above as
    $\mathbf {u} _{h}(\mathbf {x} )=a_{1}~N_{1}(\mathbf {x} )+a_{2}~N_{2}(\mathbf {x} )+\dots +a_{n}~N_{n}(\mathbf {x} )=\sum _{i=1}^{n}a_{i}~N_{i}(\mathbf {x} )$
where $n$ is the total number of nodes in the domain. Also, recall that the value of $N_{i}$ is 1 at node $i$ and zero elsewhere. Therefore, we have
    ${\begin{aligned}u_{1}&:=\mathbf {u} _{h}(\mathbf {x} _{1})=a_{1}~N_{1}(\mathbf {x} _{1})+a_{2}~N_{2}(\mathbf {x} _{1})+\dots +a_{n}~N_{n}(\mathbf {x} _{1})=a_{1}\\u_{2}&:=\mathbf {u} _{h}(\mathbf {x} _{2})=a_{1}~N_{1}(\mathbf {x} _{2})+a_{2}~N_{2}(\mathbf {x} _{2})+\dots +a_{n}~N_{n}(\mathbf {x} _{2})=a_{2}\\\vdots \\u_{n}&:=\mathbf {u} _{h}(\mathbf {x} _{n})=a_{1}~N_{1}(\mathbf {x} _{n})+a_{2}~N_{2}(\mathbf {x} _{n})+\dots +a_{n}~N_{n}(\mathbf {x} _{n})=a_{n}\\\end{aligned}}$
Therefore, the trial function can be written as
    $\mathbf {u} _{h}(\mathbf {x} )=u_{1}~N_{1}(\mathbf {x} )+u_{2}~N_{2}(\mathbf {x} )+\dots +u_{n}~N_{n}(\mathbf {x} )=\sum _{i=1}^{n}u_{i}~N_{i}(\mathbf {x} )$
where $u_{i}$ are the nodal displacements.

### Finite element system of equations

If all the elements are assumed to be of the same length $h$, the finite element system of equations ($\mathbf {K} \mathbf {a} =\mathbf {f}$) can then be written as
    ${\cfrac {AE}{h}}{\begin{bmatrix}1&-1&0&0\\-1&2&-1&0\\0&-1&2&-1\\0&0&-1&1\end{bmatrix}}{\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\cfrac {a}{h}}{\begin{bmatrix}\left({\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\\\left({\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\right)+\left({\cfrac {x_{3}(x_{3}^{2}-x_{2}^{2})}{2}}-{\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}\right)\\\left({\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}-{\cfrac {x_{2}(x_{3}^{2}-x_{2}^{2})}{2}}\right)+\left({\cfrac {x_{4}(x_{4}^{2}-x_{3}^{2})}{2}}-{\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}\right)\\\left({\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}-{\cfrac {x_{3}(x_{4}^{2}-x_{3}^{2})}{2}}\right)+{\boldsymbol {R}}{\cfrac {h}{a}}\end{bmatrix}}$

### Essential boundary conditions

To solve this system of equations we have to apply the essential boundary condition $\mathbf {u} =0$ at $\mathbf {x} =0$. This is equivalent to setting $u_{1}=0$. The reduced system of equations is
    ${\cfrac {AE}{h}}{\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&1\end{bmatrix}}{\begin{bmatrix}u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\cfrac {a}{h}}{\begin{bmatrix}\left({\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\right)+\left({\cfrac {x_{3}(x_{3}^{2}-x_{2}^{2})}{2}}-{\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}\right)\\\left({\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}-{\cfrac {x_{2}(x_{3}^{2}-x_{2}^{2})}{2}}\right)+\left({\cfrac {x_{4}(x_{4}^{2}-x_{3}^{2})}{2}}-{\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}\right)\\\left({\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}-{\cfrac {x_{3}(x_{4}^{2}-x_{3}^{2})}{2}}\right)+{\boldsymbol {R}}{\cfrac {h}{a}}\end{bmatrix}}$
This system of equations can be solved for $u_{2}$, $u_{3}$, and $u_{4}$. Let us do that.
Assume that $A$, $E$, $L$, $a$, and $R$ are all equal to 1. Then $x_{1}=0$, $x_{2}=1/3$, $x_{3}=2/3$, $x_{4}=1$, and $h=1/3$. The system of equations becomes
    ${\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&1\end{bmatrix}}{\begin{bmatrix}u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\begin{bmatrix}0.037\\0.074\\0.383\end{bmatrix}}\qquad \implies \qquad {\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\begin{bmatrix}0.0\\0.494\\0.951\\1.333\end{bmatrix}}$

### Computing element strains and stresses

From the above, it is clear that the displacement field within an element $e$ is given by
    $\mathbf {u} ^{e}=u_{1}^{e}N_{1}^{e}(\mathbf {x} )+u_{2}^{e}N_{2}^{e}(\mathbf {x} )~.$
Therefore, the strain within an element is
    ${\boldsymbol {\varepsilon }}^{e}={\frac {\partial \mathbf {u} ^{e}}{\partial x}}=u_{1}^{e}{\frac {\partial N_{1}^{e}}{\partial x}}+u_{2}^{e}{\frac {\partial N_{2}^{e}}{\partial x}}~.$
In matrix notation,
    ${\boldsymbol {\varepsilon }}^{e}=\mathbf {B} ^{e}\mathbf {u} ^{e}=\left[{\frac {\partial N_{1}^{e}}{\partial x}}~~{\frac {\partial N_{2}^{e}}{\partial x}}\right]{\begin{bmatrix}u_{1}^{e}\\u_{2}^{e}\end{bmatrix}}=\left[-{\cfrac {1}{h}}~~{\cfrac {1}{h}}\right]{\begin{bmatrix}u_{1}^{e}\\u_{2}^{e}\end{bmatrix}}$
The stress in the element is given by
    ${\boldsymbol {\sigma }}^{e}=E{\boldsymbol {\varepsilon }}^{e}~.$
For our discretization, the element stresses are
    ${\begin{aligned}{\boldsymbol {\sigma }}^{1}&=1.48\\{\boldsymbol {\sigma }}^{2}&=1.37\\{\boldsymbol {\sigma }}^{3}&=1.15\end{aligned}}$
A plot of this solution is shown in Figure 7.
* Figure 7(a). FEM vs exact solutions for displacements of an axially loaded bar.
* Figure 7(b). FEM vs exact solutions for stresses in an axially loaded bar.

#### Matlab code

The finite element code (Matlab) used to compute this solution is given below.