This report is a basic review of weighted residual methods and FEM. Some basic differential equations are used to illustrate the method. Mathematica and Matlab code written to solve numerically a first and second order ODE using FEM.
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 first and second order ODE’s. Currently I show how to use FEM to solve first 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 differential equations (ordinary or partial). It can also be used to solve non-linear differential equations but I have not yet studied how this is done. FEM is a more versatile numerical method than the finite difference methods for solving differential equations as it supports more easily different types of geometry and boundary conditions, in addition the solution of the differential 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 finite difference methods. On the other hand FEM is more mathematically complex method, and its implementation is not as straight forward as with finite difference methods.
Considering only ordinary differential equations with constant coefficients 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*}
defined 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 differential 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 differential 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 first step in solving the differential 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 coefficients (the weights) to be determined. Hence the main computational part of FEM will be focused on determining these coefficients.
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 differential 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 coefficients \(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 defined 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 find 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 differential equation) for \(N\) number of test functions, where \(N\) is the number of the unknown coefficients \(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. Different numerical schemes use different 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 defined over each element resulting from the discretization process.
In addition, in FEM, the unknown coefficients, 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 sufficient 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 sufficient.
At first, we will work with global trial functions to illustrate how to use weighted residual method.
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 effect of changing \(N\) on the result.
Given the following ODE \[ \frac {dy}{dx}-y\left ( x\right ) =0 \] defined 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 coefficients 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 coefficients \(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 find 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 satisfies 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 different 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 coefficients, 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.
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 \] defined 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 satisfies 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 coefficients. 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 differential 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 differentiable, 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 satisfies 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 coefficients \(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 different \(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.
Let us first summarize what we have done so far. Given a differential equation defined 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 coefficients 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 differential 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 differences 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 coefficients \(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 find the solution at any point between the nodes.
This diagram illustrates the above, using the first example given above to solve a differential equation \(\frac {du}{dx}-u\left ( x\right ) =0\) with the given boundary condition of \(u\left ( 0\right ) =1\) and defined over \(0\leq x\leq 1\)
Using linear interpolation, the solution \(u\left ( x\right ) \) when \(x\) is located inside first 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 first 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 first 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 finally 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 defined 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 finally 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 first part, which was to express the global approximate solution as sum of basis functions, each multiplied by an unknowns \(q\) coefficients.
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 differential equation as we did before.
Let The differential equation to solve be\[ a\frac {du}{dx}+bu(x)=f(x) \] Defined 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 differential 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 first. 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 first 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 first 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 simplifies 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 simplifies 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 first row in the global stiffness matrix for any first 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 stiffness 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 influence is \(x_{N-1}\leq x\leq x_{N}\), the above simplifies 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 simplifies 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 stiffness 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 stiffness matrix between the first 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 influence 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 simplifies 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 simplifies 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 stiffness matrix other than the first and the last row. Hence the global stiffness 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 stiffness 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 difficult to perform.
Once the load vector is calculated, the unknowns \(q_{i}\) are solved for. But before doing that, \(q_{1}\) is first 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 first row is removed (and remembering to multiply \(u_{0}\) by the first 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 first column is now removed after moving the first entry in the first 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 first order ODE. 4 first 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).
Now the FEM formulation is given for a boundary value, second order ODE with constant coefficients 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 simplifies 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 defined 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 definition 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 simplified 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<i<N\) | \(\left \{ \frac {\left ( x-x_{i-1}\right ) }{h},\frac {\left ( x_{i+1}-x\right ) }{h}\right \} \) | \(\left \{ \frac {1}{h},\frac {-1}{h}\right \} \) |
| \(N\) | \(\frac {\left ( x-x_{N-1}\right ) }{h}\) | \(\frac {1}{h}\) |
The global stiffness matrix is now constructed. For the first equation (which corresponds to the first row in the global stiffness 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 influence 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 simplifies 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 first row in the stiffness 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 stiffness 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 influence 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 simplifies 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 stiffness 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 first 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 influence 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 first 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 stiffness matrix other than the first and the last row. Now we can write the global stiffness 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 first replace \(q_{1}\) by the initial condition given in the problem, and we must remove the first 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 first element in the second row by \(u_{0}\) and removing the first 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 first element of the first row above to the RHS and removing the first 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.
First order ODE (generated by Mathematica)
Mathematica notebook which solves first order by FEM. notebook
second order ODE animation (generated by Matlab) animation
Matlab driver file This file 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 file 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); load = build_load_vector(a,u0,nNodes,xc,f,beta); % Now remove the first row and column from the stiffness matrix A(2,1) = A(2,1)*u0; load(2) = load(2)-A(2,1); A = A(2:end,2:end); load = load(2:end); % SOLVE for unknowns q = A\load; 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. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function load = build_load_vector(a,u0,nNodes,xc,f,beta) load = zeros(nNodes,1); h = xc(2)-xc(1); load(1) = -(1/h)* quad(@(x)integrand(x,f,xc(2)),xc(1),xc(2)) ; load(nNodes) = (1/h)* quad(@(x)integrand(x,f,xc(nNodes-1)),xc(nNodes-1),xc(nNodes)) - a*beta; for i = 2:nNodes-1 load(i) = quad(@(x)integrand(x,f,xc(i-1)),xc(i-1),xc(i)) - ... quad(@(x)integrand(x,f,xc(i+1)),xc(i),xc(i+1)); load(i) = (1/h)*load(i) ; end end
