PDF (letter size)

## Examples solving an ODE using ﬁnite elements method. Mathematica and Matlab implementation and animation

Sept 28,2006   Compiled on September 9, 2023 at 3:04pm

### Contents

This report is a basic review of weighted residual methods and FEM. Some basic diﬀerential equations are used to illustrate the method. Mathematica and Matlab code written to solve numerically a ﬁrst and second order ODE using FEM.

### 1 Introduction

This is a basic review of Finite Elements Methods from Mathematical point of view with examples of how it can be used to numerically solve ﬁrst and second order ODE’s. Currently I show how to use FEM to solve ﬁrst and second order ODE. I am also working on a detailed derivation and implementation using FEM to solve the 2D Poisson’s equation but this work is not yet completed.

FEM is a numerical method for solving diﬀerential equations (ordinary or partial). It can also be used to solve non-linear diﬀerential equations but I have not yet studied how this is done. FEM is a more versatile numerical method than the ﬁnite diﬀerence methods for solving diﬀerential equations as it supports more easily diﬀerent types of geometry and boundary conditions, in addition the solution of the diﬀerential equation found using FEM can be used at any point in the domain and not just on the grid points as the case is with ﬁnite diﬀerence methods. On the other hand FEM is more mathematically complex method, and its implementation is not as straight forward as with ﬁnite diﬀerence methods.

Considering only ordinary diﬀerential equations with constant coeﬃcients over the $$x$$ domain (real line).\begin {align*} a\frac {dy}{dx}+b\ y\left ( x\right ) & =f\left ( x\right ) \\ a\frac {dy}{dx}+b\ y\left ( x\right ) -f\left ( x\right ) & =0 \end {align*}

deﬁned over $$0\leq x\leq 1$$ with the boundary condition $$y\left ( 0\right ) =y_{0}$$. In the above, only when $$y\left ( x\right )$$ is the exact solution, call it $$y_{e}\left ( x\right )$$, do we have the above identity to be true.

In other words, only when $$y=\allowbreak y_{e}$$ we can write that $$a\frac {dy_{e}}{dx}+b\ y_{e}\left ( x\right ) -f\left ( x\right ) =0$$. Such a diﬀerential equations can be represented as an operator $$L$$\begin {equation} L:=\left ( y,x\right ) \rightarrow a\frac {dy}{dx}+by\left ( x\right ) -f\left ( x\right ) \tag {1} \end {equation} If we know the exact solution, call it $$y_{e}\left ( x\right )$$ then we write \begin {align*} L\left ( y_{e},x\right ) & \rightarrow a\frac {dy_{e}}{dx}+by_{e}\left ( x\right ) -f\left ( x\right ) \\ & \rightarrow 0 \end {align*}

FEM is based on the weighted residual methods (WRM) where we assume that the solution of an diﬀerential equation is the sum of weighted basis functions represented by the symbol $$\phi _{i}$$ in here. This is in essence is similar to Fourier series, where we represent a function as a weighted sums of series made up of the basis functions which happened to be in that case the sin and cosine functions.

So the ﬁrst step in solving the diﬀerential equation is to assume that the solution, called $$\tilde {y}\left ( x\right ) ,$$ can be written as$\tilde {y}\left ( x\right ) ={\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \phi _{j}\left ( x\right )$ Where $$q_{j}$$ are unknown coeﬃcients (the weights) to be determined. Hence the main computational part of FEM will be focused on determining these coeﬃcients.

When we substitute this assumed solution in the original ODE such as shown in the above example, equation (1) now becomes\begin {align*} L\left ( \tilde {y},x\right ) & =a\frac {d\tilde {y}}{dx}+b\tilde {y}\left ( x\right ) -f\left ( x\right ) \\ & =R\left ( x\right ) \end {align*}

Where $$R\left ( x\right )$$ is called the diﬀerential equation residual, which is a function over $$x$$ and in general will not be zero due to the approximate nature of our assumed solution. Our goal is to to determine the coeﬃcients $$q_{j}$$ which will make $$R\left ( x\right )$$ the minimum over the domain of the solution.

The optimal case is for $$R\left ( x\right )$$ to be zero over the domain. One method to be able to achieve this is by forcing $$R\left ( x\right )$$ to meet the following requirement ${\displaystyle \int \limits _{\Omega }} R\left ( x\right ) \ v_{i}\left ( x\right ) \ dx=0$ for all possible sets of function $$v_{i}\left ( x\right )$$ which are also deﬁned over the same domain. The functions $$v_{i}\left ( x\right )$$ are linearly independent from each others. If we can make $$R\left ( x\right )$$ satisfy the above for each one of these functions, then this implies that $$R\left ( x\right )$$ is zero. And the solution $$\tilde {y}\left ( x\right )$$ will be as close as possible to the exact solution. We will ﬁnd out in FEM that the more elements we use, the closer to the exact solution we get. This property of convergence when it comes to FEM is important, but not analyzed here.

Each one of these functions $$v_{i}\left ( x\right )$$ is called a test function (or a weight function), hence the name of this method.

In the Galerkin method of FEM, the test functions are chosen from the same class of functions as the trial functions as will be illustrated below.

By making $$R\left ( x\right )$$ satisfy the above integral equation (called the weak form of the original diﬀerential equation) for $$N$$ number of test functions, where $$N$$ is the number of the unknown coeﬃcients $$q_{i}$$, then we obtain a set of $$N$$ algebraic equations, which we can solve for $$q_{i}$$.

The above is the basic outline of all methods based on the weighted residual methods. The choice of the trial basis functions, and the choice of the test functions, determine the method used. Diﬀerent numerical schemes use diﬀerent types of trial and test functions.

In the above, the assumed solution $$\tilde {y}\left ( x\right )$$ is made up of a series of trial functions (the basis). This solution is assumed to be valid over the whole domain. This is called a global trial function. In methods such as Finite Elements and Finite volume, the domain itself is discretized, and the assumed solution is made up of a series of solutions, each of which is deﬁned over each element resulting from the discretization process.

In addition, in FEM, the unknown coeﬃcients, called $$q_{i}$$ above, have a physical meaning, they are taken as the solution values at each node. The trial functions themselves are generated by using polynomial interpolation between the nodal values. The polynomial can be linear, quadratic or cubic polynomial or higher order.

Lagrangian interpolation method is normally used for this step. The order of the polynomial is determined by the number of unknowns at the nodes. For example, if our goal is to determine the axial displacement at each node, then we have 2 unknowns, one at each end of the element. Hence a linear interpolation will be suﬃcient in this case, since a linear polynomial $$a_{0}+a_{1}x$$ contain 2 unknowns, the $$a_{1}$$ and $$a_{0}$$. If in addition to the axial displacement, we wish to also solve for the rotation at each end of the element, hence we have a total of 4 unknowns, 2 at each end of the element, which are the displacement and the rotation.

Hence in this case the minimum interpolating polynomial needed will be a cubic polynomial $$a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}$$. In the examples below, we assume that we are only solving for axial displacement, hence a linear polynomial will be suﬃcient.

At ﬁrst, we will work with global trial functions to illustrate how to use weighted residual method.

### 2 Weighted Residual method. Global trial functions.

The best way to learn how to use WRM is by working over and programming some examples.

We analyze the solution in terms of errors and the eﬀect of changing $$N$$ on the result.

#### 2.1 First example. First order ODE

Given the following ODE $\frac {dy}{dx}-y\left ( x\right ) =0$ deﬁned over $$0\leq x\leq 1$$ with the boundary condition $$y\left ( 0\right ) =1$$, we wish to solve this numerically using the WRM. This ODE has an exact solution of $$y=e^{x}$$.

The solution using WRM will always follow these steps.

STEP 1

Assume a solution that is valid over the domain $$0\leq x\leq 1$$ to be a series solution of trial (basis) functions. We start by selecting a trial functions. The assumed solution takes the form of $\tilde {y}\left ( x\right ) =\tilde {y}_{0}+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \phi _{j}\left ( x\right )$ Where $$\phi _{j}\left ( x\right )$$ is the trial function which we have to choose, and $$q_{j}$$ are the $$N$$ unknown coeﬃcients to be determined subject to a condition which will be shown below. $$\tilde {y}_{0}$$ is the assumed solution which needs to be valid only at the boundary conditions. Hence in this example, since we are given that the solution must be $$1$$ at the initial condition $$x=0$$, then $$\tilde {y}_{0}=1$$ will satisfy this boundary condition. Hence our trial solution is$\tilde {y}=1+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \phi _{j}\left ( x\right )$ STEP 2

Now we decide on what trial function $$\phi \left ( x\right )$$ to use. For this example, we can select the trial functions to be polynomials in $$x$$ or trigonometric functions. Let us choose a polynomial $$\phi _{j}\left ( x\right ) =x^{j}$$, hence our assumed solution becomes\begin {equation} \tilde {y}=1+{\sum \limits _{j=1}^{N}}q_{j}\ x^{j} \tag {1} \end {equation} Now we need to determine the coeﬃcients $$q_{j}$$, and then our solution will be complete. This is done in the following step.

STEP 3

Substituting the above assumed solution back into the original ODE, we obtain the residual $$R\left ( x\right )$$$\frac {d\tilde {y}}{dx}-\tilde {y}\left ( x\right ) =R\left ( x\right )$ $$R\left ( x\right )$$ is the ODE residual. This is the error which will result when the assumed solution is used in place of the exact solution.

Hence from (1), we ﬁnd the residual to be\begin {align} R\left ( x\right ) & =\frac {d}{dx}\left ( 1+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ x^{j}\right ) -\left ( 1+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ x^{j}\right ) \nonumber \\ & ={\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ jx^{j-1}-\left ( 1+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ x^{j}\right ) \nonumber \\ & =-1+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \left ( jx^{j-1}-x^{j}\right ) \tag {2} \end {align}

Our goal now is to reduce this residual to minimum. The way we achieve this is by requiring that the residual satisﬁes the following integral equation\begin {equation} {\displaystyle \int \limits _{0}^{1}} v_{i}\left ( x\right ) R\left ( x\right ) dx=0\ \qquad i=1\cdots N \tag {3} \end {equation} The above is a set of $$N$$ equations. The integration is carried over the whole domain, and $$v_{i}\left ( x\right )$$ is a weight (test) function, which we have to also select. Depending on the numerical scheme used, the test function will assume diﬀerent forms.

For the Galerkin method, we select the test function to be from the same family of functions as the trial (basis) functions. Hence in this example, let us select the test function to the following polynomial     \begin {equation} v_{i}\left ( x\right ) =x^{i-1} \tag {4} \end {equation} STEP 4

We now choose a value for $$N$$ and solve the set of equations generated from (3). Let us pick $$N=3,$$ hence $$R\left ( x\right )$$ becomes\begin {align*} R\left ( x\right ) & =-1+{\displaystyle \sum \limits _{j=1}^{3}} q_{j}\left ( jx^{j-1}-x^{j}\right ) \\ & =-1+q_{1}\left ( 1-x^{1}\right ) +q_{2}\ \left ( 2x-x^{2}\right ) +q_{3}\left ( 3x^{2}-x^{3}\right ) \end {align*}

Substituting the above in (3) gives\begin {align*} {\displaystyle \int \limits _{x=0}^{x=1}} x^{i-1}R\left ( x\right ) dx & =0\qquad i=1\cdots N\\{\displaystyle \int \limits _{x=0}^{x=1}} x^{i-1}\left ( -1+q\ \left ( 1-x^{1}\right ) +q_{2}\left ( 2x-x^{2}\right ) +q_{3}\ \left ( 3x^{2}-x^{3}\right ) \right ) dx & =0\qquad i=1\cdots 3 \end {align*}

The above generates $$N=3$$ equations to solve. They are\begin {align*} {\displaystyle \int \limits _{x=0}^{x=1}} \left ( -1+q_{1}\ \left ( 1-x^{1}\right ) +q_{2}\ \left ( 2x-x^{2}\right ) +q_{3}\ \left ( 3x^{2}-x^{3}\right ) \right ) dx & =0\qquad i=1\\{\displaystyle \int \limits _{x=0}^{x=1}} x\left ( -1+q_{1}\ \left ( 1-x^{1}\right ) +q_{2}\ \left ( 2x-x^{2}\right ) +q_{3}\ \left ( 3x^{2}-x^{3}\right ) \right ) dx & =0\qquad i=2\\{\displaystyle \int \limits _{x=0}^{x=1}} x^{2}\left ( -1+q_{1}\ \left ( 1-x^{1}\right ) +q_{2}\ \left ( 2x-x^{2}\right ) +q_{3}\ \left ( 3x^{2}-x^{3}\right ) \right ) dx & =0\qquad i=3 \end {align*}

Performing the integration above, gives the following 3 equations\begin {align*} -1-q_{1}\left ( -\frac {1}{2}\right ) -q_{2}\left ( -\frac {2}{3}\right ) -q_{3}\left ( -\frac {3}{4}\right ) & =0\\ -\frac {1}{2}-q_{1}\left ( -\frac {1}{6}\right ) -q_{2}\left ( -\frac {5}{12}\right ) -q_{3}\left ( -\frac {11}{20}\right ) & =0\\ -\frac {1}{3}-q_{1}\left ( -\frac {1}{12}\right ) -q_{2}\left ( -\frac {3}{10}\right ) -q_{3}\left ( -\frac {13}{30}\right ) & =0 \end {align*}

Which can be written in matrix form as$\begin {bmatrix} \frac {1}{2} & \frac {2}{3} & \frac {3}{4}\\ \frac {1}{6} & \frac {5}{12} & \frac {11}{20}\\ \frac {1}{12} & \frac {3}{10} & \frac {13}{30}\end {bmatrix}\begin {bmatrix} q_{1}\\ q_{2}\\ q_{3}\end {bmatrix} =\begin {bmatrix} 1\\ \frac {1}{2}\\ \frac {1}{3}\end {bmatrix}$ The solution is $q_{1}=\frac {72}{71},q_{2}=\frac {30}{71},q_{3}=\frac {20}{71}$ Hence our assumed series solution is now complete, using the above coeﬃcients, and from equation (1) we write\begin {align*} \tilde {y} & =1+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ x^{j}\\ \tilde {y} & =1+q_{1}\ x+q_{2}x^{2}+q_{3}x^{3} \end {align*}

Hence$\tilde {y}\left ( x\right ) =1+\frac {72}{71}x+\frac {30}{71}x^{2}+\frac {20}{71}x^{3}$ Let use compare the above solution to the exact solution $$y=e^{x}$$ by comparing the values of the solution over a number of points. This is done using the following small Mathematica code

To make this more useful, we can examine how the error changes as $$N$$ changes. The following Mathematica code determines the solution and calculates the same table as above for $$N=1\cdots 5$$

This code below plots the absolute error as $$N$$ changes. Notice that the number of peaks in the error plot is also $$N$$ which is the polynomial order (the trial solution) used to approximate the exact solution, which is to be expected.

#### 2.2 Second example. 4th order ODE

Now we will use a more complex example and repeat the above steps. We now want to numerically solve the following$\frac {d^{4}y}{dx^{4}}+y\left ( x\right ) =1$ deﬁned over $$0\leq x\leq 1\$$with the boundary conditions $$y\left ( 0\right ) =0,y\left ( 1\right ) =0,y^{\prime }\left ( 0\right ) ,y^{\prime }\left ( 1\right ) =0$$. This problem is taken from Professor S.N.Atluri text book ’Methods of computer modeling in engineering and the sciences’ Volume 1, page 47-50. Professor Atluri used a trigonometric functions for the trial function $\tilde {y}={\displaystyle \sum \limits _{i=1}^{N}} q_{i}\sin \left ( \left ( 2i-1\right ) \pi x\right )$ Which already satisﬁes the boundary conditions. For the test function, the same function as above is used hence the test (weight) function is $v_{j}\left ( x\right ) =\sin \left ( \left ( 2j-1\right ) \pi x\right )$ The book above then reduces the residual equation to a symmetric form by doing integration by parts before solving it for the coeﬃcients. In here, we will use the unsymmetrical weak form and compare the results with those shown in the above textbook. We now start again with the same steps as we did in the above example.

step 1

Selecting the trial solution.$\tilde {y}\left ( x\right ) =\tilde {y}_{0}+{\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \phi _{j}\left ( x\right )$ $$\tilde {y}_{0}=0$$ as this will satisfy the boundary conditions. Hence the trial solution is$\tilde {y}\left ( x\right ) ={\displaystyle \sum \limits _{j=1}^{N}} q_{j}\phi _{j}\left ( x\right )$ step 2

Selecting trial basis function $$\phi _{j}\left ( x\right )$$. As mentioned above, we select $$\phi _{j}\left ( x\right ) =\sin \left ( \left ( 2j-1\right ) \pi x\right )$$ , hence the trial solution is$\tilde {y}={\displaystyle \sum \limits _{j=1}^{N}} q_{j}\sin \left ( \left ( 2j-1\right ) \pi x\right )$ step 3

Substituting the above assumed solution into the original ODE, gives the diﬀerential equation residual $$R\left ( x\right )$$$\frac {d^{4}\tilde {y}}{dx^{4}}+\tilde {y}\left ( x\right ) -1=R\left ( x\right )$$\frac {d^{4}}{dx^{4}}\left ( {\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \phi _{j}\left ( x\right ) \right ) +\left ( {\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \phi _{j}\left ( x\right ) \right ) -1=R\left ( x\right )$ Notice the requirement above that the trial basis functions must be 4 times diﬀerentiable, which is the case here. From above we obtain$R\left ( x\right ) =\left ( {\displaystyle \sum \limits _{j=1}^{N}} \left ( 1+\left ( \left ( 2j-1\right ) \pi \right ) ^{4}\right ) q_{j}\sin \left ( \left ( 2j-1\right ) \pi x\right ) \right ) -1$ Our goal now is to reduce this residual to minimum. The way we achieve this is by requiring that the residual satisﬁes the following weak form integral equation\begin {equation} {\displaystyle \int \limits _{\Omega }} v_{i}\left ( x\right ) R\left ( x\right ) dx=0\qquad i=1\cdots N \tag {3} \end {equation} The above is a set of $$N$$ equations. The integration is carried over the whole domain, and $$v_{i}\left ( x\right )$$ is a weight (test) function, which we have to also select. As mentioned above, in this problem we select the test function to be \begin {equation} v_{i}\left ( x\right ) =\sin \left ( \left ( 2i-1\right ) \pi x\right ) \tag {4} \end {equation} step 4

Deciding on a value for $$N$$ and solving the set of equations generated from (3). Let us pick $$N=3,$$ hence $$R\left ( x\right )$$ becomes\begin {align*} R\left ( x\right ) & =\left ( {\displaystyle \sum \limits _{j=1}^{3}} \left ( 1+\left ( \left ( 2j-1\right ) \pi \right ) ^{4}\right ) \left ( q_{j}\sin \left ( 2j-1\right ) \pi x\right ) \right ) -1\\ & =\left ( \left ( 1+\pi ^{4}\right ) \left ( q_{1}\sin \pi x\right ) +\left ( 1+\left ( 3\pi \right ) ^{4}\right ) \left ( q_{2}\sin \left ( 3\pi x\right ) \right ) +\left ( 1+\left ( 5\pi \right ) ^{4}\right ) \left ( q_{3}\sin 5\pi x\right ) \right ) -1 \end {align*}

Hence (3) becomes\begin {align*} {\displaystyle \int \limits _{\Omega }} v_{i}\left ( x\right ) R\left ( x\right ) dx & =0\qquad i=1\cdots N\\{\displaystyle \int \limits _{0}^{1}} \left ( \sin \left ( 2i-1\right ) \pi x\right ) R\left ( x\right ) dx & =0\qquad i=1\cdots N\\{\displaystyle \int \limits _{0}^{1}} \left ( \sin \left ( 2i-1\right ) \pi x\right ) \left \{ \left ( \left ( 1+\pi ^{4}\right ) \left ( q_{1}\sin \pi x\right ) +\left ( 1+\left ( 3\pi \right ) ^{4}\right ) \left ( q_{2}\sin 3\pi x\right ) +\left ( 1+\left ( 5\pi \right ) ^{4}\right ) \left ( q_{3}\sin 5\pi x\right ) \right ) -1\right \} dx & =0\qquad i=1\cdots N \end {align*}

The above generates $$N$$ equations to solve for the coeﬃcients $$q_{i}$$\begin {align*} {\displaystyle \int \limits _{0}^{1}} \left ( \sin \pi x\right ) \ \left \{ \left ( \left ( 1+\pi ^{4}\right ) \left ( q_{1}\sin \pi x\right ) +\left ( 1+\left ( 3\pi \right ) ^{4}\right ) \left ( q_{2}\sin 3\pi x\right ) +\left ( 1+\left ( 5\pi \right ) ^{4}\right ) \left ( q_{3}\sin 5\pi x\right ) \right ) -1\right \} dx & =0\qquad i=1\\{\displaystyle \int \limits _{0}^{1}} \left ( \sin 3\pi x\right ) \ \left \{ \left ( \left ( 1+\pi ^{4}\right ) \left ( q_{1}\sin \pi x\right ) +\left ( 1+\ \ \left ( 3\pi \right ) ^{4}\right ) \left ( q_{2}\sin \left ( 3\pi x\right ) \right ) +\left ( 1+\left ( 5\pi \right ) ^{4}\right ) \left ( q_{3}\sin 5\pi x\right ) \right ) -1\right \} dx & =0\qquad i=2\\{\displaystyle \int \limits _{0}^{1}} \left ( \sin 5\pi x\right ) \ \left \{ \left ( \left ( 1+\pi ^{4}\right ) \left ( q_{1}\sin \pi x\right ) +\left ( 1+\ \ \left ( 3\pi \right ) ^{4}\right ) \left ( q_{2}\sin 3\pi x\right ) +\left ( 1+\left ( 5\pi \right ) ^{4}\right ) \left ( q_{3}\sin 5\pi x\right ) \right ) -1\right \} dx & =0\qquad i=3\ \end {align*}

Carrying the integration above and simplifying and solving for $$q_{i}$$ gives the numerical solution.

This below is a Mathematica code which solves this problem for diﬀerent $$N$$ values, and compares the error as $$N$$ changes. The error shown is the percentage error in the solution (approximate compared to exact) for up to $$N=10$$. The result below agrees well with the result in Professor’s Atluri textbook.

### 3 Finite element method

#### 3.1 Example one. First order ODE, linear interpolation

Let us ﬁrst summarize what we have done so far. Given a diﬀerential equation deﬁned over domain $$\Omega$$, we assume its solution to be of the form $$\tilde {u}\left ( x\right ) ={\displaystyle \sum \limits _{j=1}^{N}} q_{j}\ \phi _{j}\left ( x\right )$$.

The function $$\phi _{j}\left ( x\right )$$ is called the $$j^{th}$$ basis function. $$q_{j}\ \phi _{j}\left ( x\right )$$ is called a trial function.

The function $$\phi _{j}\left ( x\right )$$ is made up of functions called the shape functions$$\ N_{k\text { }}$$as they are normally called in structural mechanics books.

$$q_{j}\$$are the unknown coeﬃcients which are determined by solving $$N$$ set of equations generated by setting $$N$$ integrals of the form $${\displaystyle \int \limits _{\Omega }} R\left ( x\right ) \ v_{i}\left ( x\right ) \ dx$$ to zero. Where $$R\left ( x\right )$$ is the diﬀerential equation residual. In all what follows $$N$$ is the taken as the number of nodes.

In FEM, we also carry the same basic process as was described above, the diﬀerences are the following:

Now we divide the domain itself into a number of elements. Earlier we did not do this.

Next, the $$\phi _{j}\left ( x\right )$$ function is found by assuming the solution to be an interpolation between the nodes of the element. The solution values at the nodes are the $$q_{i}$$ and are of course unknown except at the boundaries as given by the problem.

We start by deciding on what interpolation between the nodes to use. We will use polynomial interpolation here. Then $$q_{j}\phi _{j}\left ( x\right )$$ will become the interpolation function.

In addition, the coeﬃcients $$q_{j}$$ represent the solution at the node $$j$$. These are the unknowns, which we will solve for by solving the weak form integral equation as many times as there are unknowns to solve for.

By solving for the nodal values, we can then use the interpolating function again to ﬁnd the solution at any point between the nodes.

This diagram illustrates the above, using the ﬁrst example given above to solve a diﬀerential equation $$\frac {du}{dx}-u\left ( x\right ) =0$$ with the given boundary condition of $$u\left ( 0\right ) =1$$ and deﬁned over $$0\leq x\leq 1$$

Using linear interpolation, the solution $$u\left ( x\right )$$ when $$x$$ is located inside ﬁrst element, is found from $u\left ( x\right ) =q_{1}+\text {slope}\ \left ( x-x_{0}\right )$ But the slope of the linear interpolating line over the ﬁrst element is $$\frac {q_{2}-q_{1}}{x_{2}-x_{1}}$$, hence the above becomes\begin {align*} u\left ( x\right ) & =q_{1}+\frac {q_{2}-q_{1}}{x_{2}-x_{1}}\left ( x-x_{1}\right ) \\ & =q_{1}\frac {\left ( x_{2}-x\right ) }{\left ( x_{2}-x_{1}\right ) }+q_{2}\frac {\left ( x-x_{1}\right ) }{\left ( x_{2}-x_{1}\right ) } \end {align*}

The above is the linear interpolating polynomial. We could also have used the formula of Lagrangian interpolation to arrive at the same result.

The above is the approximate solution which is valid over the ﬁrst element only. Using superscript to indicate the element number, and assuming we have equal division between nodes of length say $$h$$ (i.e. element length is $$h$$) then we write$u^{1}\left ( x\right ) =q_{1}\frac {\left ( x_{2}-x\right ) }{h}+q_{2}\frac {\left ( x-x_{1}\right ) }{h}$ Again, the above is valid for $$x_{1}\leq x\leq x_{2}$$. We now do the same for the second element

$u^{2}\left ( x\right ) =q_{2}\frac {\left ( x_{3}-x\right ) }{h}+q_{3}\frac {\left ( x-x_{2}\right ) }{h}$ The above is valid for $$x_{2}\leq x\leq x_{3}$$. And ﬁnally for the 3rd element$u^{3}\left ( x\right ) =q_{3}\frac {\left ( x_{4}-x\right ) }{h}+q_{4}\frac {\left ( x-x_{3}\right ) }{h}$ The above is valid for $$x_{2}\leq x\leq x_{3}$$. This is now illustrated in the following diagram

Since our goal is to express the global approximate solution $$u\left ( x\right )$$ as a series sum of basis functions each multiplied by $$q_{j}$$, we now rewrite each of the $$u^{j}$$ to allow this, as follows

\begin {align*} u^{1}\left ( x\right ) & =q_{1}\overset {N_{1}^{1}\left ( x\right ) }{\overbrace {\frac {\left ( x_{2}-x\right ) }{h}}}+q_{2}\overset {N_{2}^{1}(x)}{\overbrace {\frac {\left ( x-x_{1}\right ) }{h}}}\\ & =q_{1}N_{1}^{1}\left ( x\right ) +q_{2}N_{2}^{1}(x) \end {align*}

The above is valid for $$x_{1}\leq x\leq x_{2}$$. Notice the use of the following notation: Since each element will have deﬁned on it two shape functions, $$N_{1}\left ( x\right )$$ and $$N_{2}\left ( x\right )$$, one per node, then we use a superscript to indicate the element number. Hence for element $$1$$, we will write its two shapes functions as $$N_{1}^{1}\left ( x\right )$$ and $$N_{2}^{1}(x)$$.

We now do the same for the second element\begin {align*} u^{2} & =q_{2}\overset {N_{1}^{2}\left ( x\right ) }{\overbrace {\frac {\left ( x_{3}-x\right ) }{h}}}+q_{3}\overset {N_{2}^{2}\left ( x\right ) }{\overbrace {\frac {\left ( x-x_{2}\right ) }{h}}}\\ & =q_{2}N_{1}^{2}\left ( x\right ) +q_{3}N_{2}^{2}\left ( x\right ) \end {align*}

The above is valid for $$x_{2}\leq x\leq x_{3}$$. And ﬁnally for the 3rd element\begin {align*} u^{3} & =q_{3}\overset {N_{1}^{3}\left ( x\right ) }{\overbrace {\frac {\left ( x_{4}-x\right ) }{h}}}+q_{4}\overset {N_{2}^{3}\left ( x\right ) }{\overbrace {\frac {\left ( x-x_{3}\right ) }{h}}}\\ & =q_{3}N_{1}^{3}+q_{4}N_{2}^{3} \end {align*}

The above is valid for $$x_{2}\leq x\leq x_{3}$$. The global trial function is\begin {align*} u(x) & =u^{1}+u^{2}+u^{3}\\ & =\left ( q_{1}N_{1}^{1}+q_{2}N_{2}^{1}\right ) +\left ( q_{2}N_{1}^{2}+q_{3}N_{2}^{2}\right ) +\left ( q_{3}N_{1}^{3}+q_{4}N_{2}^{3}\right ) \\ & =q_{1}N_{1}^{1}+q_{2}\left ( N_{2}^{1}+N_{1}^{2}\right ) +q_{3}\left ( N_{2}^{2}+N_{1}^{3}\right ) +q_{4}\left ( N_{2}^{3}\right ) \end {align*}

The shape function for node 1 is \begin {align*} \phi _{1} & =N_{1}^{1}\\ & =\frac {x_{2}-x}{h} \end {align*}

And the shape function for node 2 is \begin {align*} \phi _{2} & =N_{2}^{1}+N_{1}^{2}\\ & =\frac {\left ( x-x_{1}\right ) }{h}+\frac {\left ( x_{3}-x\right ) }{h} \end {align*}

The shape function for node 3 is \begin {align*} \phi _{3} & =N_{2}^{2}+N_{1}^{3}\\ & =\frac {\left ( x-x_{2}\right ) }{h}+\frac {\left ( x_{4}-x\right ) }{h} \end {align*}

The shape function for the last node is \begin {align*} \phi _{4} & =N_{2}^{3}\\ & =\frac {x-x_{3}}{h} \end {align*}

We see that the shape function for any internal node is $\phi _{j}=\frac {x-x_{j-1}}{h}+\frac {x_{j+1}-x}{h}$ The approximate solution is therefore$u\left ( x\right ) ={\sum \limits _{i=1}^{Number\ Nodes}}q_{i}\phi _{i}$ This completes the ﬁrst part, which was to express the global approximate solution as sum of basis functions, each multiplied by an unknowns $$q$$ coeﬃcients.

The diagram below illustrates the above.

The diagram below illustrates the numbering used.

Next the residual $$R\left ( x\right )$$ is found by substituting the above global solution into the original diﬀerential equation as we did before.

Let The diﬀerential equation to solve be$a\frac {du}{dx}+bu(x)=f(x)$ Deﬁned over $$0\leq x\leq 1$$ with the initial condition $$u(0)=u_{0}$$.

$$N$$ is the number of nodes, therefore the residual is \begin {align*} R\left ( x\right ) & =a\frac {d}{dx}{\sum \limits _{i=1}^{N}}q_{i}\phi _{i}+b{\sum \limits _{i=1}^{N}}q_{i}\phi _{i}-f\left ( x\right ) \\ R\left ( x\right ) & ={\sum \limits _{i=1}^{N}}q_{i}\left ( a\phi _{i}^{\prime }+b\phi _{i}\right ) -f\left ( x\right ) \end {align*}

The test functions are $v_{j}=\phi _{j}\qquad j=1\cdots N$ The weak form of the diﬀerential equation is\begin {align*} {\int \limits _{x=x_{1}}^{x=x_{N}}}\phi _{j}\left ( x\right ) R\left ( x\right ) dx & =0\qquad j=1\cdots N\\ {\int \limits _{x=x_{1}}^{x=x_{N}}}\phi _{j}\left ( x\right ) \left ( \left [ {\sum \limits _{i=1}^{N}}q_{i}\left ( a\phi _{i}^{\prime }+b\phi _{i}\right ) \right ] -f\left ( t\right ) \right ) dx & =0\qquad j=1\cdots N \end {align*}

$$N$$ equations are obtained from the above, which are solved for the $$q_{i}$$.

The derivatives of the shape functions are found ﬁrst. Assuming in this example that the domain is divided into 3 elements results in

 $$i$$ $$\phi _{i}$$ $$\phi _{i}^{\prime }$$ $$1$$ $$\frac {\left ( x_{2}-x\right ) }{h}$$ $$\frac {-1}{h}$$ $$2$$ $$\left \{ \frac {\left ( x-x_{1}\right ) }{h},\frac {\left ( x_{3}-x\right ) }{h}\right \}$$ $$\left \{ \frac {1}{h},\frac {-1}{h}\right \}$$ $$3$$ $$\frac {\left ( x-x_{2}\right ) }{h}$$ $$\frac {1}{h}$$

In general, if there are $$N$$ nodes, then for the ﬁrst node \begin {align*} \phi _{1} & =\frac {\left ( x_{2}-x\right ) }{h}\\ \phi _{1}^{\prime } & =\frac {-1}{h} \end {align*}

And for the last node \begin {align*} \phi _{N} & =\frac {\left ( x-x_{N-1}\right ) }{h}\\ \phi _{N}^{\prime } & =\frac {1}{h} \end {align*}

And for any node in the middle \begin {align*} \phi _{i} & =\frac {\left ( x-x_{i-1}\right ) }{h}\\ \phi _{i}^{\prime } & =\frac {1}{h} \end {align*}

for $$x_{i-1}\leq x\leq x_{i}$$ and \begin {align*} \phi _{i} & =\frac {\left ( x_{i+1}-x\right ) }{h}\\ \phi _{i}^{\prime } & =\frac {-1}{h} \end {align*}

for $$x_{i}\leq x\leq x_{i+1}$$.

Looking back at the weak form integral above, it is evaluated as follows${\int \limits _{x=x_{1}}^{x=x_{N}}}\phi _{j}\left ( \left [ {\sum \limits _{i=1}^{N}}q_{i}\left ( a\phi _{i}^{\prime }+b\phi _{i}\right ) \right ] -f\left ( x\right ) \right ) dx=0\qquad j=1\cdots N$ For the ﬁrst node only, $$j=1$$, the following results$\int \limits _{x=x_{1}}^{x=x_{N}}\phi _{1}\left ( \left [ \sum \limits _{i=1}^{N}q_{i}\left ( a\phi _{i}^{\prime }+b\phi _{i}\right ) \right ] -f\left ( x\right ) \right ) dx=0$ Since $$\phi _{1}$$ is not zero only over $$x_{1}\leq x\leq x_{2}$$, and given that $$\phi _{1}=\frac {x_{2}-x}{h},\phi _{1}^{\prime }=\frac {-1}{h},\phi _{2}=\frac {\left ( x-x_{1}\right ) }{h}$$ due to the range of integration limits, and that $$\phi _{2}^{\prime }=\frac {1}{h}$$, then the above can simpliﬁes to\begin {align*} {\int \limits _{x=x_{1}}^{x=x_{2}}}\phi _{1}\left ( \left [ q_{1}\left ( a\phi _{1}^{\prime }+b\phi _{1}\right ) +q_{2}\left ( a\phi _{2}^{\prime }+b\phi _{2}\right ) \right ] -f\left ( x\right ) \right ) dx & =0\ \ \ \ \ \\ {\int \limits _{x=x_{1}}^{x=x_{2}}}\frac {\left ( x_{2}-x\right ) }{h}\left ( \left [ q_{1}\left ( \frac {-a}{h}+b\frac {\left ( x_{2}-x\right ) }{h}\right ) +q_{2}\left ( a\frac {1}{h}+b\frac {\left ( x-x_{1}\right ) }{h}\right ) \right ] -f\left ( x\right ) \right ) dx & =0 \end {align*}

The above simpliﬁes further to \begin {align*} -q_{1}\left ( \frac {\left ( 3a+2b\left ( x_{1}-x_{2}\right ) \right ) \left ( x_{1}-x_{2}\right ) ^{2}}{6h^{2}}\right ) -q_{2}\frac {\left ( -3a+b\left ( x_{1}-x_{2}\right ) \right ) \left ( x_{1}-x_{2}\right ) ^{2}}{6h^{2}}-{\int \limits _{x=x_{1}}^{x=x_{2}}}\phi _{1}f\left ( x\right ) & =0\\ -q_{1}\left ( \frac {\left ( 3a+2b\left ( x_{1}-x_{2}\right ) \right ) \left ( x_{1}-x_{2}\right ) ^{2}}{6h^{2}}\right ) -q_{2}\frac {\left ( -3a+b\left ( x_{1}-x_{2}\right ) \right ) \left ( x_{1}-x_{2}\right ) ^{2}}{6h^{2}}-\frac {1}{h}{\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x_{2}-x\right ) f\left ( x\right ) & =0 \end {align*}

Since $$x_{2}-x_{1}=h$$ and $$x_{1}-x_{2}=-h$$ the above reduces to\begin {align*} -q_{1}\left ( \frac {\left ( 3a-2bh\right ) h^{2}}{6h^{2}}\right ) -q_{2}\frac {\left ( -3a-bh\right ) h^{2}}{6h^{2}}-\frac {1}{h}{\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x_{2}-x\right ) f\left ( x\right ) & =0\\ q_{1}\left ( \frac {2bh-3a}{6}\right ) +q_{2}\left ( \frac {bh+3a}{6}\right ) -\frac {1}{h}{\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x_{2}-x\right ) f\left ( x\right ) & =0 \end {align*}

The above equation gives the ﬁrst row in the global stiﬀness matrix for any ﬁrst order linear ODE of the form $$a\frac {du}{dx}+bu(x)=f(x)$$. The above shows that numerical integration is only needed to be performed on the term $${\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x_{2}-x\right ) f\left ( x\right )$$.

Next the last equation is found, which will be the last row of the stiﬀness matrix. For the last node only $$j=N$$ the following results${\int \limits _{x=x_{N-1}}^{x=x_{N}}}\phi _{N}\left ( \left [ {\sum \limits _{i=1}^{N}}q_{i}\left ( a\phi _{i}^{\prime }+b\phi _{i}\right ) \right ] -f\left ( x\right ) \right ) dx=0$ Since $$\phi _{N}$$ domain of inﬂuence is $$x_{N-1}\leq x\leq x_{N}$$, the above simpliﬁes to${\int \limits _{x=x_{N-1}}^{x=x_{N}}}\phi _{N}\left ( \left [ q_{N-1}\left ( a\phi _{N-1}^{\prime }+b\phi _{N-1}\right ) +q_{N}\left ( a\phi _{N}^{\prime }+b\phi _{N}\right ) \right ] -f\left ( x\right ) \right ) dx=0$ Since $$\phi _{N-1}=\frac {\left ( x_{N}-x\right ) }{h},\phi _{N-1}^{\prime }=\frac {-1}{h},\phi _{N}=\frac {x-x_{N-1}}{h},\phi _{N}^{\prime }=\frac {1}{h}$$ the above becomes${\int \limits _{x=x_{N-1}}^{x=x_{N}}}\frac {x-x_{N-1}}{h}\left ( \left [ q_{N-1}\left ( \frac {-a}{h}+b\frac {\left ( x_{N}-x\right ) }{h}\right ) +q_{N}\left ( a\frac {1}{h}+b\frac {x-x_{N-1}}{h}\right ) \right ] -f\left ( x\right ) \right ) dx=0$ Which simpliﬁes to$N-1\frac {\left ( 3a+b\left ( x_{N-1}-x_{N}\right ) \right ) \left ( x_{N-1}-x_{N}\right ) ^{2}}{6h^{2}}-q_{N}\frac {\left ( -3a+2b\left ( x_{N-1}-x_{N}\right ) \right ) \left ( x_{N-1}-x_{N}\right ) ^{2}}{6h^{2}}-\frac {1}{h}\int \limits _{x=x_{N-1}}^{x=x_{N}}\left ( x-x_{N-1}\right ) f(x)\,dx=0$ Letting $$x_{N-1}-x_{N}=-h$$ in the above becomes$q_{N-1}\left ( \frac {bh-3a}{6}\right ) +q_{N}\frac {\left ( 3a+2bh\right ) }{6}-\frac {1}{h}{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x-x_{N-1}\right ) f\left ( x\right ) =0$ Hence the last row of the stiﬀness matrix can be determined directly except for the term under the integral which needs to be evaluated using integration.

Now the equation that represents any internal node is found. This will be any row in the global stiﬀness matrix between the ﬁrst and the last row.

For any $$j$$ other than 1 or $$N$$ the following results

$\int \limits _{x=x_{j-1}}^{x=x_{j}}\phi _{j}\left ( \left [ {\sum \limits _{i=1}^{N}}q_{i}\left ( a\phi _{i}^{\prime }+b\phi _{i}\right ) \right ] -f\left ( x\right ) \right ) dx+\ \int \limits _{x=x_{j}}^{x=x_{j+1}}\phi _{j}\left ( \left [ {\sum \limits _{i=1}^{N}}q_{i}\left ( a\phi _{i}^{\prime }+b\phi _{i}\right ) \right ] -f\left ( x\right ) \right ) dx=0$ Where the integral was broken into two parts to handle the domain of inﬂuence of the shape functions.${\int \limits _{x=x_{j-1}}^{x=x_{j}}}\phi _{j}\left ( \ \ q_{j-1}\left ( a\phi _{j-1}^{\prime }+b\phi _{j-1}\right ) +q_{j}\left ( a\phi _{j}^{\prime }+b\phi _{j}\right ) -f\left ( x\right ) \right ) dx+{\int \limits _{x=x_{j}}^{x=x_{j+1}}}\phi _{j}\left ( q_{j}\left ( a\phi _{j}^{\prime }+b\phi _{j}\right ) +q_{j+1}\left ( a\phi _{j+1}^{\prime }+b\phi _{j+1}\right ) -f\left ( x\right ) \right ) dx=0$ For \begin {align*} x_{j-1} & \leq x\leq x_{j},\phi _{j}=\frac {x-x_{j-1}}{h},\phi _{j}^{\prime }=\frac {1}{h},\phi _{j-1}=\frac {x_{j}-x}{h},\phi _{j-1}^{\prime }=\frac {-1}{h}\\ x_{j} & \leq x\leq x_{j+1},\phi _{j}=\frac {x_{j+1}-x}{h},\phi _{j}^{\prime }=\frac {-1}{h},\phi _{j+1}=\frac {x-x_{j}}{h},\phi _{j+1}^{\prime }=\frac {1}{h} \end {align*}

Hence the weak form integral can be written as \begin {multline*} {\int \limits _{x=x_{j-1}}^{x=x_{j}}}\frac {x-x_{j-1}}{h}\left ( q_{j-1}\left ( \frac {-a}{h}+b\frac {x_{j}-x}{h}\right ) +q_{j}\left ( \frac {a}{h}+b\frac {x-x_{j-1}}{h}\right ) -f\left ( x\right ) \right ) dx\\ +{\int \limits _{x=x_{j}}^{x=x_{j+1}}}\frac {x_{j+1}-x}{h}\left ( q_{j}\left ( \frac {-a}{h}+b\frac {x_{j+1}-x}{h}\right ) +q_{j+1}\left ( \frac {a}{h}+b\frac {x-x_{j}}{h}\right ) -f\left ( x\right ) \right ) dx=0 \end {multline*} Which simpliﬁes to\begin {multline*} q_{j-1}\left ( \frac {\left ( -3a-b\left ( x_{j-1}-x_{j}\right ) \right ) \left ( x_{j-1}-x_{j}\right ) ^{2}}{6h^{2}}\right ) +q_{j}\frac {\left ( \left ( 3a-2b\left ( x_{j-1}-x_{j}\right ) \right ) \left ( x_{j-1}-x_{j}\right ) ^{2}\right ) }{6h^{2}}-\frac {1}{h}\int \limits _{x=x_{j-1}}^{x=x_{j}}\left ( x-x_{j-1}\right ) f\left ( x\right ) dx+\\ q_{j}\left ( \frac {\left ( -3a-2b\left ( x_{j}-x_{j+1}\right ) \right ) \left ( x_{j}-x_{j+1}\right ) ^{2}}{6h^{2}}\right ) +q_{j+1}\left ( \frac {\left ( 3a-b\left ( x_{j}-x_{j+1}\right ) \right ) \left ( x_{j}-x_{j+1}\right ) ^{2}}{6h^{2}}\right ) -\frac {1}{h}\int \limits _{x=x_{j}}^{x=x_{j+1}}\left ( x_{j+1}-x\right ) f\left ( x\right ) dx=0 \end {multline*} For equal distance between elements, $$x_{j}-x_{j+1}=-h,x_{j-1}-x_{j}=-h$$ the above simpliﬁes to\begin {multline*} q_{j-1}\left ( \frac {\left ( -3a+bh\right ) h^{2}}{6h^{2}}\right ) +q_{j}\frac {\left ( \left ( 3a+2bh\right ) h^{2}\right ) }{6h^{2}}-\frac {1}{h}\int \limits _{x=x_{j-1}}^{x=x_{j}}\left ( x-x_{j-1}\right ) f\left ( x\right ) dx+q_{j}\left ( \frac {\left ( -3a+2bh\right ) h^{2}}{6h^{2}}\right ) \\ +q_{j+1}\left ( \frac {\left ( 3a+bh\right ) h^{2}}{6h^{2}}\right ) -\frac {1}{h}\int \limits _{x=x_{j}}^{x=x_{j+1}}\left ( x_{j+1}-x\right ) f\left ( x\right ) dx=0 \end {multline*} Combining gives$q_{j-1}\left ( \frac {-3a+bh}{6}\right ) +q_{j}\left ( \frac {2bh}{3}\right ) +q_{j+1}\left ( \frac {3a+bh}{6}\right ) -\frac {1}{h}\left ( {\int \limits _{x=x_{j-1}}^{x=x_{j}}}\left ( x-x_{j-1}\right ) f\left ( x\right ) dx+\int \limits _{x=x_{j}}^{x=x_{j+1}}\left ( x_{j+1}-x\right ) f\left ( x\right ) dx\right ) =0$ The above gives the expression for any row in the stiﬀness matrix other than the ﬁrst and the last row. Hence the global stiﬀness matrix is$\begin {bmatrix} \frac {2bh-3a}{6} & \frac {bh+3a}{6} & 0 & 0 & 0 & \ldots & 0\\ \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & 0 & \ldots & 0\\ 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & \ldots & 0\\ 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0\\ 0 & 0 & 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6}\\ 0 & 0 & 0 & 0 & 0 & \frac {bh-3a}{6} & \frac {3a+2bh}{6}\end {bmatrix}\begin {bmatrix} q_{1}\\ q_{2}\\ q_{3}\\ \vdots \\ \vdots \\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} =\begin {bmatrix} \frac {1}{h}{\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x_{2}-x\right ) f\left ( x\right ) \\ \frac {1}{h}\left ( {\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x-x_{1}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{2}}^{x=x_{3}}}\left ( x_{3}-x\right ) f\left ( x\right ) dx\right ) \\ \frac {1}{h}\left ( {\int \limits _{x=x_{2}}^{x=x_{3}}}\left ( x-x_{2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{3}}^{x=x_{4}}}\left ( x_{4}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{j-1}}^{x=x_{j}}}\left ( x-x_{j-1}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{j}}^{x=x_{j+1}}}\left ( x_{j+1}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{N-2}}^{x=x_{N-1}}}\left ( x-x_{N-2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x_{N}-x\right ) f\left ( x\right ) dx\right ) \\ \frac {1}{h}{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x-x_{N-1}\right ) f\left ( x\right ) \end {bmatrix}$ The above shows that the stiﬀness matrix can be build quickly without the use of any numerical integration in this example. Integration is needed to evaluate the force (or load) vector. Depending on the forcing function, this can be simple or diﬃcult to perform.

Once the load vector is calculated, the unknowns $$q_{i}$$ are solved for. But before doing that, $$q_{1}$$ is ﬁrst replaced by the initial condition given in the problem$\begin {bmatrix} \frac {2bh-3a}{6} & \frac {bh+3a}{6} & 0 & 0 & 0 & \ldots & 0\\ \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & 0 & \ldots & 0\\ 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & \ldots & 0\\ 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0\\ 0 & 0 & 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6}\\ 0 & 0 & 0 & 0 & 0 & \frac {bh-3a}{6} & \frac {3a+2bh}{6}\end {bmatrix}\begin {bmatrix} u_{0}\\ q_{2}\\ q_{3}\\ \vdots \\ \vdots \\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} =\begin {bmatrix} \frac {1}{h}{\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x_{2}-x\right ) f\left ( x\right ) \\ \frac {1}{h}\left ( {\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x-x_{1}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{2}}^{x=x_{3}}}\left ( x_{3}-x\right ) f\left ( x\right ) dx\right ) \\ \frac {1}{h}\left ( {\int \limits _{x=x_{2}}^{x=x_{3}}}\left ( x-x_{2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{3}}^{x=x_{4}}}\left ( x_{4}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{j-1}}^{x=x_{j}}}\left ( x-x_{j-1}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{j}}^{x=x_{j+1}}}\left ( x_{j+1}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{N-2}}^{x=x_{N-1}}}\left ( x-x_{N-2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x_{N}-x\right ) f\left ( x\right ) dx\right ) \\ \frac {1}{h}{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x-x_{N-1}\right ) f\left ( x\right ) \end {bmatrix}$ Now the ﬁrst row is removed (and remembering to multiply $$u_{0}$$ by the ﬁrst entry in the second row) gives$\begin {bmatrix} \frac {-3a+bh}{6}u_{0} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & 0 & \ldots & 0\\ 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & \ldots & 0\\ 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0\\ 0 & 0 & 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6}\\ 0 & 0 & 0 & 0 & 0 & \frac {bh-3a}{6} & \frac {3a+2bh}{6}\end {bmatrix}\begin {bmatrix} q_{2}\\ q_{3}\\ \vdots \\ \vdots \\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} =\begin {bmatrix} \frac {1}{h}\left ( {\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x-x_{1}\right ) f\left ( x\right ) dx+x_{3}-xf\left ( x\right ) dx\right ) \\ \frac {1}{h}\left ( {\int \limits _{x=x_{2}}^{x=x_{3}}}\left ( x-x_{2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{3}}^{x=x_{4}}}\left ( x_{4}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{j-1}}^{x=x_{j}}}\left ( x-x_{j-1}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{j}}^{x=x_{j+1}}}\left ( x_{j+1}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{N-2}}^{x=x_{N-1}}}\left ( x-x_{N-2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x_{N}-x\right ) f\left ( x\right ) dx\right ) \\ \frac {1}{h}{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x-x_{N-1}\right ) f\left ( x\right ) \end {bmatrix}$ The ﬁrst column is now removed after moving the ﬁrst entry in the ﬁrst row to the RHS to become part of the load vector$\begin {bmatrix} \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & 0 & \ldots & 0\\ \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0 & \ldots & 0\\ 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6} & 0\\ 0 & 0 & 0 & \frac {-3a+bh}{6} & \frac {2bh}{3} & \frac {3a+bh}{6}\\ 0 & 0 & 0 & 0 & \frac {bh-3a}{6} & \frac {3a+2bh}{6}\end {bmatrix}\begin {bmatrix} q_{2}\\ q_{3}\\ \vdots \\ \vdots \\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} =\begin {bmatrix} \frac {1}{h}\left ( {\int \limits _{x=x_{1}}^{x=x_{2}}}\left ( x-x_{1}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{2}}^{x=x_{3}}}\left ( x_{3}-x\right ) f\left ( x\right ) dx\right ) -\left ( \frac {-3a+bh}{6}u_{0}\right ) \\ \frac {1}{h}\left ( {\int \limits _{x=x_{2}}^{x=x_{3}}}\left ( x-x_{2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{3}}^{x=x_{4}}}\left ( x_{4}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{j-1}}^{x=x_{j}}}\left ( x-x_{j-1}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{j}}^{x=x_{j+1}}}\left ( x_{j+1}-x\right ) f\left ( x\right ) dx\right ) \\ \vdots \\ \frac {1}{h}\left ( {\int \limits _{x=x_{N-2}}^{x=x_{N-1}}}\left ( x-x_{N-2}\right ) f\left ( x\right ) dx+{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x_{N}-x\right ) f\left ( x\right ) dx\right ) \\ \frac {1}{h}{\int \limits _{x=x_{N-1}}^{x=x_{N}}}\left ( x-x_{N-1}\right ) f\left ( x\right ) \end {bmatrix}$ Now the system above is solved for $$q_{2},q_{3},\cdots ,q_{N}$$. Once the $$q_{i}$$ are found, the solution is$u\left ( x\right ) ={\sum \limits _{i=1}^{Number\ Nodes}}q_{i}\phi _{i}$ The Appendix shows a Mathematica code which solve the above general ﬁrst order ODE. 4 ﬁrst order ODE solved for illustration and the solution is compared to the exact solution. Animation was made to help illustrate this. the RMS error was calculated. The animations can be access by clicking on the links below (shows only in the HTML version of this report).

#### 3.2 Example Two. 2nd order ODE, Boundary value problem. Linear interpolation. Symmetric weak form.

Now the FEM formulation is given for a boundary value, second order ODE with constant coeﬃcients of the form$a\frac {d^{2}u}{dx^{2}}+b\frac {du}{dx}+cu\left ( x\right ) =f\left ( x\right )$ with the initial conditions $$u\left ( x_{0}\right ) =u_{0}$$ (called essential Dirichlet condition), and the Neumann boundary condition $$\frac {du}{dx}=\beta$$ at $$L$$, for $$x_{0}\leq x\leq L$$

As before, the solution is assumed to be of the form $u={\sum \limits _{i=1}^{N}}q_{i}\phi _{i}$ And the unknowns $$q_{i}$$ are found by solving $$N$$ algebraic equations${\int \limits _{x_{0}}^{L}}\phi _{j}R\left ( x\right ) dx=0\qquad j=1\cdots N$ Where $$R\left ( x\right )$$ is the ODE residual obtained by substituting the assumed solution into the original ODE. Hence the above becomes (For simplicity, $$j=1\cdots N$$ is not written each time as it is assumed to be the case)\begin {align*} {\int \limits _{x_{0}}^{L}}\phi _{j}\left ( a\frac {d^{2}u}{dx^{2}}+b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx & =0\ \\ {\int \limits _{x_{0}}^{L}}\phi _{j}a\frac {d^{2}u}{dx^{2}}dx+{\int \limits _{x_{0}}^{L}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx & =0 \end {align*}

Applying integration by parts on $${\int \limits _{x_{0}}^{L}}\phi _{j}a\frac {d^{2}u}{dx^{2}}dx$$. Since $\frac {d}{dx}\left ( a\ \phi _{j}\frac {du}{dx}\right ) =a\phi _{j}\frac {d^{2}u}{dx^{2}}+a\phi _{j}^{\prime }\frac {du}{dx}$ then integrating both sides of the above gives\begin {align*} {\int \limits _{x_{0}}^{L}}\frac {d}{dx}\left ( a\ \phi _{j}\frac {du}{dx}\right ) dx & ={\int \limits _{x_{0}}^{L}}a\phi _{j}\frac {d^{2}u}{dx^{2}}dx+{\int \limits _{x_{0}}^{L}}a\phi _{j}^{\prime }\frac {du}{dx}dx\\ \left [ a\phi _{j}\frac {du}{dx}\right ] _{x_{0}}^{L} & ={\int \limits _{x_{0}}^{L}}a\phi _{j}\frac {d^{2}u}{dx^{2}}dx+{\int \limits _{x_{0}}^{L}}a\phi _{j}^{\prime }\frac {du}{dx}dx \end {align*}

Hence ${\int \limits _{x_{0}}^{L}}a\phi _{j}\frac {d^{2}u}{dx^{2}}dx=\left [ a\phi _{j}\frac {du}{dx}\right ] _{x_{0}}^{L}-{\int \limits _{x_{0}}^{L}}a\phi _{j}^{\prime }\frac {du}{dx}dx$ Substituting the above in equation A1 gives\begin {equation} \left [ a\phi _{j}\frac {du}{dx}\right ] _{x_{0}}^{L}-{\int \limits _{x_{0}}^{L}}a\phi _{j}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{0}}^{L}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx=0\ \tag {A2} \end {equation} Considering the term $\left [ a\phi _{j}\frac {du}{dx}\right ] _{x_{0}}^{L}=\left . a\phi _{j}\frac {du}{dx}\right \vert _{x=L}-\left . a\phi _{j}\frac {du}{dx}\right \vert _{x=x_{0}}$ First considering the term $$\left . a\phi _{j}\frac {du}{dx}\right \vert _{x=L}$$. Since $$\frac {du}{dx}=\beta$$ at $$x=L$$ then this term becomes $$a\beta \left [ \phi _{j}\right ] _{x=L}$$. But $$\left [ \phi _{j}\right ] _{x=L}$$ is non zero only for the $$N^{th}$$ shape function evaluated at $$x=L$$. Since linear interpolation is used, the $$N^{th}$$ shape function is $$\phi _{N}=\frac {x-x_{N-1}}{h}$$ which have the value of $$1$$ at $$x=x_{N}$$. Hence $\left . a\phi _{j}\frac {du}{dx}\right \vert _{x=L}=a\beta \qquad \text {when}j=N,\text {otherwise,}0$ Now considering the term $$\left . a\phi _{j}\frac {du}{dx}\right \vert _{x=x_{0}}$$, since at $$x=x_{0}$$ all shape functions will be zero except for $$\phi _{1}$$ which has the value of $$1$$ at $$x=x_{0}$$. Hence this simpliﬁes to $$\left . a\frac {du}{dx}\right \vert _{x=x_{0}}$$

Recalling that $$\frac {du}{dx}=\frac {d}{dx}{\sum \limits _{i=1}^{N}}q_{i}\phi _{i}={\sum \limits _{i=1}^{N}}q_{i}\phi _{i}^{\prime }$$. But at $$x=x_{0}$$ only $$\phi _{1}$$ is deﬁned and its has the slope of $$\frac {-1}{h}$$, hence $\left . a\frac {du}{dx}\right \vert _{x=x_{0}}=\frac {-aq_{1}}{h}\ \text {when}\ j=1\text {, otherwise }0$ However, $$q_{1}=u_{0}$$ since that is by deﬁnition the initial condition. Finally one obtains $\left [ a\phi _{j}\frac {du}{dx}\right ] _{x_{0}}^{L}=\left \{ \begin {array} [c]{c}\frac {aq_{1}}{h}\qquad j=1\\ a\beta \qquad j=N \end {array} \right \}$ Hence the symmetric weak form equation (A2) can now be simpliﬁed more giving\begin {equation} \left \{ \begin {array} [c]{c}\frac {aq_{1}}{h}\qquad j=1\\ a\beta \qquad j=N \end {array} \right \} -{\int \limits _{x_{0}}^{L}}a\phi _{j}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{0}}^{L}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx=0\qquad j=1\cdots N \tag {A3} \end {equation} The trial function obtain by linear interpolation are used. These are shown in this table

 $$i$$ $$\phi _{i}$$ $$\phi _{i}^{\prime }$$ $$1$$ $$\frac {\left ( x_{2}-x\right ) }{h}$$ $$\frac {-1}{h}$$ $$1 The global stiﬀness matrix is now constructed. For the ﬁrst equation (which corresponds to the ﬁrst row in the global stiﬀness matrix) the result is$\frac {aq_{1}}{h}-{\int \limits _{x_{0}}^{L}}a\phi _{1}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{0}}^{L}}\phi _{1}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx=0$ Since \(u={\sum \limits _{i=1}^{N}}q_{i}\phi _{i}$$, and since $$\phi _{1}$$ domain of inﬂuence is only from $$x_{1}$$ to $$x_{2}$$ then above becomes$\frac {aq_{1}}{h}-{\int \limits _{x_{1}}^{x_{2}}}a\phi _{1}^{\prime }\frac {d}{dx}\left ( q_{1}\phi _{1}+q_{2}\phi _{2}\right ) dx+{\int \limits _{x_{1}}^{x2}}\phi _{1}\left ( b\frac {d}{dx}\left ( q_{1}\phi _{1}+q_{2}\phi _{2}\right ) +c\left ( q_{1}\phi _{1}+q_{2}\phi _{2}\right ) -f\left ( x\right ) \right ) dx=0$ But $$\phi _{1}=\frac {\left ( x_{2}-x\right ) }{h}$$ and $$\phi _{1}^{\prime }=\frac {-1}{h}$$ and over the domain from $$x_{1}$$ to $$x_{2},$$ $$\phi _{2}=\frac {\left ( x-x_{1}\right ) }{h},\phi _{2}^{\prime }=\frac {1}{h}$$, hence the above becomes$\frac {aq_{1}}{h}-{\int \limits _{x_{1}}^{x_{2}}}\frac {-a}{h}\left ( \frac {-q_{1}}{h}+\frac {q_{2}}{h}\right ) dx+{\int \limits _{x_{1}}^{x2}}\frac {\left ( x_{2}-x\right ) }{h}\left ( b\left ( \frac {-q_{1}}{h}+\frac {q_{2}}{h}\right ) +c\left ( q_{1}\frac {\left ( x_{2}-x\right ) }{h}+q_{2}\frac {\left ( x-x_{1}\right ) }{h}\right ) -f\left ( x\right ) \right ) dx=0$ The above simpliﬁes to\begin {align*} \frac {aq_{1}}{h}-\frac {q_{1}}{h^{2}}\left ( ax_{1}-\frac {b(x_{1}-x_{2})^{2}}{2}-\frac {c(x_{1}-x_{2})^{3}}{3}-ax_{2}\right ) +\frac {q_{2}}{h^{2}}\left ( -ax_{1}+\frac {b\left ( x_{1}-x_{2}\right ) ^{2}}{2}-\frac {c\left ( x_{1}-x_{2}\right ) ^{3}}{6}+ax_{2}\right ) -\frac {1}{h}\int \limits _{x_{1}}^{x2}\left ( x_{2}-x\right ) f\left ( x\right ) & =0\\ -\frac {q_{1}}{h^{2}}\left ( -ah+ax_{1}-\frac {b\left ( x_{1}-x_{2}\right ) ^{2}}{2}-\frac {c\left ( x_{1}-x_{2}\right ) ^{3}}{3}-ax_{2}\right ) +\frac {q_{2}}{h^{2}}\left ( -ax_{1}+\frac {b\left ( x_{1}-x_{2}\right ) ^{2}}{2}-\frac {c\left ( x_{1}-x_{2}\right ) ^{3}}{6}+ax_{2}\right ) -\frac {1}{h}\int \limits _{x_{1}}^{x2}\left ( x_{2}-x\right ) f\left ( x\right ) & =0 \end {align*}

The above gives the ﬁrst row in the stiﬀness matrix. Since it is assumed that each element will have the same length, hence $$x_{1}-x_{2}=-h,$$ and the above becomes\begin {align*} -\frac {q_{1}}{h^{2}}\left ( -ah+ax_{1}-\frac {bh^{2}}{2}+\frac {ch^{3}}{3}-ax_{2}\right ) +\frac {q_{2}}{h^{2}}\left ( -ax_{1}+\frac {bh^{2}}{2}+\frac {ch^{3}}{6}+ax_{2}\right ) -\frac {1}{h}{\int \limits _{x_{1}}^{x2}}\left ( x_{2}-x\right ) f\left ( x\right ) & =0\\ -q_{1}\left ( -\frac {a}{h}+\frac {ax_{1}}{h^{2}}-\frac {b}{2}+\frac {ch}{3}-\frac {ax_{2}}{h^{2}}\right ) +q_{2}\left ( -\frac {ax_{1}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{2}}{h^{2}}\right ) -\frac {1}{h}{\int \limits _{x_{1}}^{x2}}\left ( x_{2}-x\right ) f\left ( x\right ) & =0 \end {align*}

For the last equation, which will be the last row in the global stiﬀness matrix$a\beta -{\int \limits _{x_{N-1}}^{x_{N}}}a\phi _{N}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{N-1}}^{x_{N}}}\phi _{N}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx=0$ Since $$u={\sum \limits _{i=1}^{N}}q_{i}\phi _{i}$$, and since $$\phi _{N}$$ domain of inﬂuence is only from $$x_{N-1}$$ to $$x_{N}$$ The above becomes$a\beta -{\int \limits _{x_{N-1}}^{x_{N}}}a\phi _{N}^{\prime }\frac {d}{dx}\left ( q_{N-1}\phi _{N-1}+q_{N}\phi _{N}\right ) dx+{\int \limits _{x_{N-1}}^{x_{N}}}\phi _{N}\left ( b\frac {d}{dx}\left ( q_{N-1}\phi _{N-1}+q_{N}\phi _{N}\right ) +c\left ( q_{N-1}\phi _{N-1}+q_{N}\phi _{N}\right ) -f\left ( x\right ) \right ) dx=0$ But $$\ \phi _{N}=\frac {\left ( x-x_{N-1}\right ) }{h}$$ and $$\phi _{N}^{\prime }=\frac {1}{h}$$ and over the domain from $$x_{N-1}$$ to $$x_{N},$$ $$\phi _{N-1}=\frac {\left ( x_{N}-x\right ) }{h},\phi _{N-1}^{\prime }=\frac {-1}{h}$$ hence the above becomes$a\beta -{\int \limits _{x_{N-1}}^{x_{N}}}\frac {a}{h}\left ( q_{N-1}\left ( \frac {-1}{h}\right ) +q_{N}\frac {1}{h}\right ) dx+{\int \limits _{x_{N-1}}^{x_{N}}}\frac {\left ( x-x_{N-1}\right ) }{h}\left ( b\left ( q_{N-1}\left ( \frac {-1}{h}\right ) +q_{N}\frac {1}{h}\right ) +c\left ( q_{N-1}\frac {\left ( x_{N}-x\right ) }{h}+q_{N}\frac {\left ( x-x_{N-1}\right ) }{h}\right ) -f\left ( x\right ) \right ) dx=0$ The above simpliﬁes to$a\beta -\frac {q_{N-1}}{h^{2}}\left ( -ax_{N-1}-\frac {b\left ( x_{N-1}-x_{N}\right ) ^{2}}{2}-\frac {c\left ( x_{N-1}-x_{N}\right ) ^{3}}{6}+ax_{N}\right ) +\frac {q_{N}}{h^{2}}\left ( +ax_{N-1}+\frac {b\left ( x_{N-1}-x_{N}\right ) ^{2}}{2}-\frac {c\left ( x_{N-1}-x_{N}\right ) ^{3}}{3}-ax_{N}\right ) -\frac {1}{h}\int \limits _{x_{N-1}}^{x_{N}}\left ( x-x_{N-1}\right ) f\left ( x\right ) =0$ The above represents the last row in the global stiﬀness matrix. Since it is assumed that each element will have the same length, hence $$x_{N-1}-x_{N}=-h,$$ and the above becomes\begin {align*} a\beta -\frac {q_{N-1}}{h^{2}}\left ( -ax_{N-1}-\frac {bh^{2}}{2}+\frac {ch^{3}}{6}+ax_{N}\right ) +\frac {q_{N}}{h^{2}}\left ( ax_{N-1}+\frac {bh^{2}}{2}+\frac {ch^{3}}{3}-ax_{N}\right ) -\frac {1}{h}{\int \limits _{x_{N-1}}^{x_{N}}}\left ( x-x_{N-1}\right ) f\left ( x\right ) & =0\\ a\beta -q_{N-1}\left ( -\frac {ax_{N-1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{N}}{h^{2}}\right ) +q_{N}\left ( \frac {ax_{N-1}}{h^{2}}+\frac {b}{2}+\frac {ch}{3}-\frac {ax_{N}}{h^{2}}\right ) -\frac {1}{h}{\int \limits _{x_{N-1}}^{x_{N}}}\left ( x-x_{N-1}\right ) f\left ( x\right ) & =0 \end {align*}

The expression for any row in between the ﬁrst and the last rows is now found. For a general node $$j$$ the result is$-{\int \limits _{x_{j-1}}^{x_{j+1}}}a\phi _{j}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{j-1}}^{x_{j+1}}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx=0$ Breaking the integral into halves to make it easier to write the trial functions over the domain of inﬂuence gives$-\left ( {\int \limits _{x_{j-1}}^{x_{j}}}a\phi _{j}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{j}}^{x_{j+1}}}a\phi _{j}^{\prime }\frac {du}{dx}dx\right ) +{\int \limits _{x_{j-1}}^{x_{j}}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx+{\int \limits _{x_{j}}^{x_{j+1}}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx=0$ Considering the ﬁrst domain $$x_{j-1}\leq x\leq x_{j}$$ gives$-{\int \limits _{x_{j-1}}^{x_{j}}}a\phi _{j}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{j-1}}^{x_{j}}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx$ Over this range, $$u=q_{j-1}\phi _{j-1}+q_{j}\phi _{j}$$ hence $$\frac {du}{dx}=q_{j-1}\phi _{j-1}^{\prime }+q_{j}\phi _{j}^{\prime }$$ where $$\phi _{j-1}=\frac {x_{j}-x}{h},\phi _{j-1}^{\prime }=\frac {-1}{h},\phi _{j}=\frac {x-x_{j-1}}{h},\phi _{j}^{\prime }=\frac {1}{h}$$ hence the above becomes\begin {equation} -{\int \limits _{x_{j-1}}^{x_{j}}}a\frac {1}{h}\left ( q_{j-1}\left ( \frac {-1}{h}\right ) +q_{j}\frac {1}{h}\right ) dx+{\int \limits _{x_{j-1}}^{x_{j}}}\frac {x-x_{j-1}}{h}\left ( b\left ( q_{j-1}\left ( \frac {-1}{h}\right ) +q_{j}\frac {1}{h}\right ) +c\left ( q_{j-1}\frac {x_{j}-x}{h}+q_{j}\frac {x-x_{j-1}}{h}\right ) -f\left ( x\right ) \right ) dx \tag {A4} \end {equation} Now considering the second domain $$x_{j}\leq x\leq x_{j+1}$$ gives$-{\int \limits _{x_{j}}^{x_{j+1}}}a\phi _{j}^{\prime }\frac {du}{dx}dx+{\int \limits _{x_{j}}^{x_{j+1}}}\phi _{j}\left ( b\frac {du}{dx}+cu\left ( x\right ) -f\left ( x\right ) \right ) dx$ Over this range, $$u=q_{j}\phi _{j}+q_{j+1}\phi _{j+1}$$ hence $$\frac {du}{dx}=q_{j}\phi _{j}^{\prime }+q_{j+1}\phi _{j+1}^{\prime }$$ where $$\phi _{j}=\frac {x_{j+1}-x}{h},\phi _{j}^{\prime }=\frac {-1}{h},\phi _{j+1}=\frac {x-x_{j}}{h},\phi _{j+1}^{\prime }=\frac {1}{h}$$ hence the above becomes\begin {equation} -\int \limits _{x_{j}}^{x_{j+1}}\left ( \frac {-a}{h}\right ) \left ( \frac {-q_{j}}{h}+\frac {q_{j+1}}{h}\right ) dx+\int \limits _{x_{j}}^{x_{j+1}}\frac {x_{j+1}-x}{h}\left ( b\left ( \frac {-q_{j}}{h}+q_{j+1}\frac {1}{h}\right ) +c\left ( q_{j}\frac {x_{j+1}-x}{h}+q_{j+1}\frac {x-x_{j}}{h}\right ) -f(x)\right ) dx\nonumber \end {equation} Combine A4 and A5 and simplifying gives\begin {multline*} \frac {q_{j-1}}{h^{2}}\left ( -ax_{j-1}-\frac {b\left ( x_{j-1}-x_{j}\right ) ^{2}}{2}-\frac {c\left ( x_{j-1}-x_{j}\right ) ^{3}}{6}+ax_{j}\right ) +\frac {q_{j}}{h^{2}}\left ( ax_{j-1}+\frac {b\left ( x_{j-1}-x_{j}\right ) ^{2}}{2}-\frac {c\left ( x_{j-1}-x_{j}\right ) ^{3}}{3}-\frac {b\left ( x_{j}-x_{j+1}\right ) ^{2}}{2}-\frac {c\left ( x_{j}-x_{j+1}\right ) ^{3}}{3}-ax_{j+1}\right ) \\ +\frac {q_{j+1}}{h^{2}}\left ( -ax_{j}+\frac {b\left ( x_{j}-x_{j+1}\right ) ^{2}}{2}-\frac {c\left ( x_{j}-x_{j+1}\right ) ^{3}}{6}+ax_{j+1}\right ) -{\displaystyle \int \limits _{x_{j-1}}^{x_{j}}} \frac {x-x_{j-1}}{h}f\left ( x\right ) dx-{\displaystyle \int \limits _{x_{j}}^{x_{j+1}}} \frac {x_{j+1}-x}{h}dx=0 \end {multline*} Since we are assuming each element will have the same length, hence $$x_{j-1}-x_{j}=-h,$$ $$x_{j}-x_{j+1}=-h$$ and the above becomes\begin {multline*} q_{j-1}\left ( -\frac {ax_{j-1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{j}}{h^{2}}\right ) +q_{j}\left ( \frac {ax_{j-1}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{j+1}}{h^{2}}\right ) \\ +q_{j+1}\left ( \frac {-ax_{j}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{j+1}}{h^{2}}\right ) -{\displaystyle \int \limits _{x_{j-1}}^{x_{j}}} \frac {x-x_{j-1}}{h}f\left ( x\right ) dx-{\displaystyle \int \limits _{x_{j}}^{x_{j+1}}} \frac {x_{j+1}-x}{h}f\left ( x\right ) dx=0 \end {multline*} The above gives any row in the stiﬀness matrix other than the ﬁrst and the last row. Now we can write the global stiﬀness matrix as\begin {multline*} \begin {bmatrix} \left ( -\frac {a}{h}+\frac {ax_{1}}{h^{2}}-\frac {b}{2}+\frac {ch}{3}-\frac {ax_{2}}{h^{2}}\right ) & \left ( -\frac {ax_{1}}{h^{2}}+\frac {b}{2}-\frac {ch}{6}+\frac {ax_{2}}{h^{2}}\right ) & 0 & 0\\ \left ( -\frac {ax_{1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{2}}{h^{2}}\right ) & \left ( \frac {ax_{1}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{3}}{h^{2}}\right ) & \left ( \frac {-ax_{2}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{3}}{h^{2}}\right ) & 0\\ 0 & \left ( -\frac {ax_{j-1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{j}}{h^{2}}\right ) & \left ( \frac {ax_{j-1}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{j+1}}{h^{2}}\right ) & \left ( \frac {-ax_{j}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{j+1}}{h^{2}}\right ) \\ 0 & 0 & \ddots & \cdots \\ 0 & 0 & \left ( -\frac {ax_{N-1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{N}}{h^{2}}\right ) & \left ( \frac {ax_{N-1}}{h^{2}}+\frac {b}{2}+\frac {ch}{3}-\frac {ax_{N}}{h^{2}}\right ) \end {bmatrix}\begin {bmatrix} q_{1}\\ q_{2}\\ \vdots \\ \vdots \\ q_{j}\\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} \\ =\begin {bmatrix} \frac {1}{h}{\displaystyle \int \limits _{x_{1}}^{x2}} \left ( x_{2}-x\right ) f\left ( x\right ) \\ \vdots \\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{j-1}}^{x_{j}}} \left ( x-x_{j-1}\right ) f\left ( x\right ) dx+\frac {1}{h}{\displaystyle \int \limits _{x_{j}}^{x_{j+1}}} \left ( x_{j+1}-x\right ) f\left ( x\right ) dx\\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{N-1}}^{x_{N}}} \left ( x-x_{N-1}\right ) f\left ( x\right ) -a\beta \end {bmatrix} \end {multline*}

We must ﬁrst replace $$q_{1}$$ by the initial condition given in the problem, and we must remove the ﬁrst row after that as follows\begin {multline*} \begin {bmatrix} \left ( -\frac {a}{h}+\frac {ax_{1}}{h^{2}}-\frac {b}{2}+\frac {ch}{3}-\frac {ax_{2}}{h^{2}}\right ) & \left ( -\frac {ax_{1}}{h^{2}}+\frac {b}{2}-\frac {ch}{6}+\frac {ax_{2}}{h^{2}}\right ) & 0 & 0\\ \left ( -\frac {ax_{1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{2}}{h^{2}}\right ) & \left ( \frac {ax_{1}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{3}}{h^{2}}\right ) & \left ( \frac {-ax_{2}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{3}}{h^{2}}\right ) & 0\\ 0 & 0 & \ddots & \cdots \\ 0 & 0 & \ddots & \cdots \\ 0 & 0 & \left ( -\frac {ax_{N-1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{N}}{h^{2}}\right ) & \left ( \frac {ax_{N-1}}{h^{2}}+\frac {b}{2}+\frac {ch}{3}-\frac {ax_{N}}{h^{2}}\right ) \end {bmatrix}\begin {bmatrix} u_{0}\\ q_{2}\\ \vdots \\ \vdots \\ q_{j}\\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} \\ =\begin {bmatrix} \frac {1}{h}{\displaystyle \int \limits _{x_{1}}^{x2}} \left ( x_{2}-x\right ) f\left ( x\right ) \\ \vdots \\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{j-1}}^{x_{j}}} \left ( x-x_{j-1}\right ) f\left ( x\right ) dx+\frac {1}{h}{\displaystyle \int \limits _{x_{j}}^{x_{j+1}}} \left ( x_{j+1}-x\right ) f\left ( x\right ) dx\\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{N-1}}^{x_{N}}} \left ( x-x_{N-1}\right ) f\left ( x\right ) -a\beta \end {bmatrix} \end {multline*} Multiply the ﬁrst element in the second row by $$u_{0}$$ and removing the ﬁrst row gives\begin {multline*} \begin {bmatrix} \left ( -\frac {ax_{1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{2}}{h^{2}}\right ) u_{0} & \left ( \frac {ax_{1}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{3}}{h^{2}}\right ) & \left ( \frac {-ax_{2}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{3}}{h^{2}}\right ) & 0\\ 0 & \left ( -\frac {ax_{2}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{3}}{h^{2}}\right ) & \left ( \frac {ax_{2}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{4}}{h^{2}}\right ) & \left ( \frac {-ax_{3}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{4}}{h^{2}}\right ) \\ 0 & 0 & \ddots & \ddots \\ 0 & 0 & \left ( -\frac {ax_{N-1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{N}}{h^{2}}\right ) & \left ( \frac {ax_{N-1}}{h^{2}}+\frac {b}{2}+\frac {ch}{3}-\frac {ax_{N}}{h^{2}}\right ) \end {bmatrix}\begin {bmatrix} q_{2}\\ \vdots \\ \vdots \\ q_{j}\\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} \\ =\begin {bmatrix} \frac {1}{h}{\displaystyle \int \limits _{x_{1}}^{x_{2}}} \left ( x-x_{2}\right ) f\left ( x\right ) dx+\frac {1}{h}{\displaystyle \int \limits _{x_{2}}^{x_{3}}} \left ( x_{3}-x\right ) dx\\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{j-1}}^{x_{j}}} \left ( x-x_{j-1}\right ) f\left ( x\right ) dx+\frac {1}{h}{\displaystyle \int \limits _{x_{j}}^{x_{j+1}}} \left ( x_{j+1}-x\right ) f\left ( x\right ) dx\\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{N-1}}^{x_{N}}} \left ( x-x_{N-1}\right ) f\left ( x\right ) -a\beta \end {bmatrix} \end {multline*} Now moving the ﬁrst element of the ﬁrst row above to the RHS and removing the ﬁrst column gives\begin {multline*} \begin {bmatrix} \left ( \frac {ax_{1}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{3}}{h^{2}}\right ) & \left ( \frac {-ax_{2}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{3}}{h^{2}}\right ) & 0 & 0\\ 0 & \left ( -\frac {ax_{2}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{3}}{h^{2}}\right ) & \left ( \frac {ax_{2}}{h^{2}}+\frac {2ch}{3}-\frac {ax_{4}}{h^{2}}\right ) & \left ( \frac {-ax_{j}}{h^{2}}+\frac {b}{2}+\frac {ch}{6}+\frac {ax_{j+1}}{h^{2}}\right ) \\ \vdots & \ddots & \cdots & \ddots \\ 0 & 0 & \left ( -\frac {ax_{N-1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{N}}{h^{2}}\right ) & \left ( \frac {ax_{N-1}}{h^{2}}+\frac {b}{2}+\frac {ch}{3}-\frac {ax_{N}}{h^{2}}\right ) \end {bmatrix}\begin {bmatrix} q_{2}\\ \vdots \\ \vdots \\ q_{j}\\ \vdots \\ q_{N-1}\\ q_{N}\end {bmatrix} \\ =\begin {bmatrix} \frac {1}{h}{\displaystyle \int \limits _{x_{1}}^{x_{2}}} \left ( x-x_{2}\right ) f\left ( x\right ) dx+\frac {1}{h}{\displaystyle \int \limits _{x_{2}}^{x_{3}}} \left ( x_{3}-x\right ) dx-\left ( -\frac {ax_{1}}{h^{2}}-\frac {b}{2}+\frac {ch}{6}+\frac {ax_{2}}{h^{2}}\right ) u_{0}\\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{j-1}}^{x_{j}}} \left ( x-x_{j-1}\right ) f\left ( x\right ) dx+\frac {1}{h}{\displaystyle \int \limits _{x_{j}}^{x_{j+1}}} \left ( x_{j+1}-x\right ) f\left ( x\right ) dx\\ \vdots \\ \vdots \\ \frac {1}{h}{\displaystyle \int \limits _{x_{N-1}}^{x_{N}}} \left ( x-x_{N-1}\right ) f\left ( x\right ) -a\beta \end {bmatrix} \end {multline*} Now we solve for the vector $$\left [ q_{2},q_{3},\cdots ,q_{N}\right ] ^{T}$$ and this completes our solution.

A Matlab implementation is below which solves any second order ODE. This code was used to generate an animation. This animation can be accessed by clicking on the link below.

### 4 References

1. Methods of computer modeling in engineering and the sciences. Volume 1. By Professor Satya N. Atluri. Tech Science Press.
2. Class lecture notes. MAE 207. Computational methods. UCI. Spring 2006. Instructor: Professor SN Atluri.
3. Computational techniques for ﬂuid dynamics, Volume I. C.A.J.Fletcher. Springer-Verlag

### 5 Animation

First order ODE (generated by Mathematica)

Mathematica notebook which solves ﬁrst order by FEM. notebook

second order ODE animation (generated by Matlab) animation

Matlab driver ﬁle This ﬁle sets up the call to do FEM and the animation

function nma_fem_driver
%script to do a simulation of FEM solution of a second order ODE
%by Nasser Abbasi
%This script calls the m file called nma_fem.m
%
%this solves the ODE
%
%   a u''(t) + b u'(t) + c u(t) = f(t)
%
% with t over the range t0 to len
% and with initial conditions u(0)=u0
% and with u'(len)=beta
%
%To use this to solve an ODE, edit the values below for the variables
% a,b,c,f,len,u0,x0

close all; clear all;

%   equation 1. u''+u'+u=sin(x)*cos(x). u(0)=1, u'(len)=beta, len=10
%   maxy=4 miny=-1
%   a=1, b=1, c=1, x0=0, u0=1,len=10

a=1;
b=1;
c=1;
x0=0;   % starting domain point
u0=-10;
len=10;
beta=1;

ymax=4; ymin=-2;  % for plotting only. To make the plot looks better

%This below solves the equation exactly to compare the FEM solution
%against.

sol= dsolve('a*D2y+b*Dy+c*y=cos(t)*sin(t)','y(0)=-10','Dy(10)=1');
sol=subs(sol);

directory='matlab_ANIMATION/';
extension='.png';

MAX_NODES = 100;
for n = 3:MAX_NODES
ezplot(sol,[0:len,-ymin:ymax]);
hold on;

r = nma_fem(a,b,c,@f,x0,u0,beta,len,sol,n);  %SOLVES by FEM

title(sprintf('Solving y''''+y''+y(x)=sin(x)*cos(x).  y(0)=1, y''(10)=1. by FEM. \nN=%d\nRMSERROR=%f',n,r));
drawnow;

image = getframe(gcf);
p=frame2im(image);
file_name=[directory num2str(n) extension] ;
imwrite(p,file_name,'png');

hold off;
% pause(.5)
end
end

%%%%%%%%%%%%%%%%%%%%%%%%
%
% edit this below to change the forcing function.
%%%%%%%%%%%%%%%%%%%%%%%%
function v=f(x)
v=sin(x).*cos(x);
end


Matlab function to solve ODE using FEM This ﬁle does FEM solution on general second order ODE

function rmserror=nma_fem(a,b,c,f,x0,u0,beta,len,exact_sol,nNodes)
% Solves numerically a second order ODE with constant coeff.
% using the symmetric Garlekin Finite Elements Methods.
%
% see the file nma_fem_driver.m on how to call this function.
%
% Solve a u''(t) + b u'(t) + c u(t) = f(t)
%
% with t over the range t0 to len
% and with initial conditions u(0)=u0
% and with u'(len)=beta

%by Nasser Abbasi. Sept 26,2006.

xc=linspace(x0,len,nNodes);

% This plots the shape functions.
%    x=linspace(x0,len,1000);
%    y=linspace(x0,len,1000);
%    for i=1:nNodes
%        for j=1:length(y)
%            [v,d]=phi(i,x(j),nNodes,x0,len,xc);
%            y(j)=v;
%        end
%        plot(x,y);
%        hold on;
%    end

A    = build_stiffness_matrix(nNodes,xc,a,b,c);

% Now remove the first row and column from the stiffness matrix

A(2,1)  = A(2,1)*u0;
A       = A(2:end,2:end);

% SOLVE for unknowns

q = [u0;q];

% Plot the solution
y = zeros(length(xc),1);
for i=1:length(xc)
y(i)=trial(xc(i),x0,len,xc,nNodes,q);
end

plot(xc,y,'ro');
hold on;
line(xc,y);

% Calculate RMSerror. Use 50 points. Should be more than enough.
NPOINTS = 50;
x = linspace(x0,len,NPOINTS);
rmserror = 0;
for i = 1:length(x)
y = trial(x(i),x0,len,xc,nNodes,q);
t=x(i);  % USED for subs below. do not remove.
yexact   = real(double(subs(exact_sol)));
rmserror = rmserror+(y-yexact)^2;
end

rmserror = rmserror/NPOINTS;
rmserror = sqrt(rmserror);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = trial(x,x0,len,xc,nNodes,q)
if x<x0||x>len
error('in Trial. x outside range');
end
v = 0;
for i = 1:nNodes
[s,d] = phi(i,x,nNodes,x0,len,xc); %notice ignore d here.
v     = v+ q(i)*s;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [v,d]=phi(i,x,nNodes,x0,len,xc)
if(i<1 || i>nNodes)
error('node number outside range');
end
if(x<x0 || x>len)
error('x outside range');
end

h = xc(2)-xc(1);

if i==1
if(x>xc(2))
v = 0;
d = 0;
else
v = (xc(2)-x)/h;
d = -1/h;
end
return;
end

if i==nNodes
if(x<xc(nNodes-1))
v = 0;
d = 0;
else
v = (x-xc(nNodes-1))/h;
d = 1/h;
end
return;
end

if(x>xc(i+1) || x<xc(i-1) )
v = 0;
d = 0;
else
if(x<=xc(i))
v = (x-xc(i-1))/h;
d = 1/h;
else
v = (xc(i+1)-x)/h;
d = -1/h;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEE my report for description of this function.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=build_stiffness_matrix(nNodes,xc,a,b,c)

A = zeros(nNodes,nNodes);
h = xc(2)-xc(1);
for i = 1:nNodes
if i==1
A(1,1) = -a/h+ a*xc(1)/h^2 - b/h +c*h/3 - a*xc(2)/h^2;
A(1,2) = -a*xc(1)/h^2 +b/2 -c*h/6 + a*xc(2)/h^2;
else
if i==nNodes
A(nNodes,nNodes-1) = -a*xc(nNodes-1)/h^2 -b/2 +c*h/6 +a*xc(nNodes)/h^2;
A(nNodes,nNodes)   = a*xc(nNodes-1)/h^2 +b/2 +c*h/3 -a*xc(nNodes)/h^2;
else
A(i,i-1) = -a*xc(i-1)/h^2 -b/2 + c*h/6 + a*xc(i)/h^2;
A(i,i)   = a*xc(i-1)/h^2 + 2*c*h/3 -a*xc(i+1)/h^2;
A(i,i+1) = -a*xc(i)/h^2 + b/2 + c*h/6 + a*xc(i+1)/h^2;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Called to integrate the 'force' function
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v=integrand(x,f,xj)
v = (x-xj).*f(x);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SEE my report for more description of this
% This build the load vector.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%