up
PDF

Examples solving an ODE using finite elements method. Mathematica and Matlab implementation and animation

Nasser M. Abbasi

Sept 28,2006 page compiled on August 19, 2015 at 10:35am

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

pict

defined over 0 ≤ x ≤ 1   with the boundary condition y (0) = y0    . In the above, only when y (x)   is the exact solution, call it ye (x)   , do we have the above identity to be true.

In other words, only when y = y
     e  we can write that adye+ b y  (x ) − f (x) = 0
  dx      e   . Such a differential equations can be represented as an operator L

              dy
L :=  (y,x ) → a--- +by (x) − f (x )
              dx
(1)

If we know the exact solution, call it y (x)
 e   then we write

pict

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 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 y (x) ,  can be written as

      ∑N
y (x) =   qj   j (x)
       j=1

Where qj  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

pict

Where R (x )   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 qj  which will make R (x)   the minimum over the domain of the solution.

The optimal case is for R (x )   to be zero over the domain. One method to be able to achieve this is by forcing R (x )   to meet the following requirement

∫
  R (x) vi (x ) dx = 0

for all possible sets of function vi (x)   which are also defined over the same domain. The functions vi (x )   are linearly independent from each others. If we can make R (x )   satisfy the above for each one of these functions, then this implies that R (x)   is zero. And the solution y (x )   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 vi (x)   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 (x)   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 qi  , then we obtain a set of N algebraic equations, which we can solve for qi  .

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 y (x)   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 descritized, 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 qi  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 displacement at the nodes, 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 a0 + a1x contain 2 unknowns, the a1    and a0    . If in addition to the 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 the minimum interpolating polynomial needed will be a cubic polynomial a  + a x + a x + a x
 0    1     2    3 . In the examples below, we assume that we are only solving for 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.

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 effect of changing N on the result.

2.1 First example. First order ODE

Given the following ODE

dy
---− y (x) = 0
dx

defined over 0 ≤ x ≤ 1   with the boundary condition y (0) = 1   , we wish to solve this numerically using the WRM. This ODE has an exact solution of      x
y = e  .

The solution using WRM will always follow these steps.

STEP 1

Assume a solution that is valid over the domain 0 ≤ x ≤ 1   to be a series solution of trial (basis) functions. We start by selecting a trial functions. The assumed solution takes the form of

           ∑N
y (x ) = y0 +    qj  j (x )
           j=1

Where j (x)   is the trial function which we have to choose, and qj  are the N unknown coefficients to be determined subject to a condition which will be shown below. y0    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 y0 = 1   will satisfy this boundary condition. Hence our trial solution is

y = 1 + N∑ q     (x)
           j j
       j=1

STEP 2

Now decide on what trial function (x )   to use. For this example, we can select the trial functions to be polynomials in x or trigonometric functions. Let use choose a polynomial j (x) = xj  , hence our assumed solution becomes

       N
      ∑      j
y = 1 +   qj x
       j=1
(1)

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

Substitute the above assumed solution into the original ODE, we obtain the residual R (x)

d-y-−  y (x ) = R (x)
dx

R (x)   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

pict

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

∫

  vi (x ) R (x)dx = 0         i = 1 ⋅⋅⋅N
(3)

The above is a set of N equations. The integration is carried over the whole domain, and vi (x)   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

vi (x ) = xi−1
(4)

STEP 4

Decide on a value for N and solve the set of equations generated from (3). Let us pick N = 3,  hence R (x)   becomes

pict

Substitute the above in (3) we obtain

pict

The above generates N equations to solve, there are

pict

Now carry the integration above, we obtain the following 3 equations

pict

Which can be written in matrix form as

⌊| 1  2   3 ⌋|    ⌊|q1⌋|   ⌊|1⌋|
|| 21  35   411||    ||  ||   ||1||
|| 6  12  20||    ||q2|| = ||2||
||-1  3-  13||    ||q ||   ||1||
⌈12  10  30⌉    ⌈ 3⌉   ⌈3⌉

The solution is

     72       30      20
q1 = ---,q2 = --,q3 = ---
     71       71      71

Hence our assumed series solution is now complete, using the above coefficients, and from equation (1) we write

pict

Hence

|-----------------------------|
 |y (x) = 1 + 72x + 30x 2 + 20-x3
------------71----71-----71----

Let use compare the above solution to the exact solution y = ex  by comparing the values of the solution over a number of points. This is done using the following small Mathematica code

pict

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 ⋅⋅⋅5

pict

pict

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.

pict

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

d4y
---4 + y (x) = 1
dx

defined over 0 ≤ x ≤ 1   with the boundary conditions                    ′      ′
y (0) = 0,y (1) = 0,y (0) ,y (1) = 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

y = ∑N q  (sin (2i − 1)  x )
       i
   i=1

Which already satisfies the boundary conditions. For the test function, the same function as above is used hence the test (weight) function is

vj (x ) = sin (2j − 1) x

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

Select the trial solution.

           ∑N
y (x ) = y0 +    qj  j (x )
           j=1

y0 = 0   as this will satisfy the boundary conditions. Hence the trial solution is

      ∑N
y (x) =   qj   j (x)
       j=1

step 2

Select trial basis function j (x )   . As mentioned above, we select j (x) = sin (2j − 1) x , hence the trial solution is

   N∑

y =    qj (sin (2j − 1)   x)
   j=1

step 3

Substitute the above assumed solution into the original ODE, we obtain the differential equation residual R (x)

d4 y
 ----+ y (x) − 1 = R (x)
dx 4

  4  N               N
d---⟨∑          ⟩  ⟨∑          ⟩
dx4 ∕   qj  j (x )\+ ∕   qj   j (x)\− 1 = R (x)
    ∕j=1        \  ∕j=1        \

Notice the requirement above that the trial basis functions must be 4 times differentiable, which is the case here. From above we obtain

         ∑N (                 ) (               )
R (x) = ⟨∕    1 +  ((2j − 1)  )4   qj sin (2j − 1) x  ⟩\− 1
        ∕j=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

∫
  v  (x ) R (x)dx = 0         i = 1 ⋅⋅⋅N
    i
(3)

The above is a set of N equations. The integration is carried over the whole domain, and vi (x)   is a weight (test) function, which we have to also select. As mentioned above, in this problem we select the test function to be

vi (x ) = (sin (2i − 1)  x)
(4)

step 4

Decide on a value for N and solve the set of equations generated from (3). Let us pick N = 3,  hence R (x)   becomes

pict

Hence (3) becomes

                               ∫

                                 vi (x ) R (x)dx = 0         i = 1 ⋅⋅⋅N

                         ∫ 1

                            (sin(2i − 1)  x ) R (x) dx = 0         i = 1⋅⋅⋅N
                          0
∫1
                   { ((      4)             (         4)               (         4)           )    }
   (sin (2i − 1)  x)     1 +       (q1 sin   x ) + 1 +  (3   )  (q2sin3   x) +  1 +  (5   )  (q3 sin 5  x ) − 1  dx = 0         i = 1 ⋅⋅⋅N
0

The above generates N equations to solve for the coefficients qi

pict

Carry the integration above and simplify and solve for qi  we obtain 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.

pict

3 Finite element method

3.1 Example one. First order ODE, linear interpolation

Let us first summarize what we have discussed so far. Given a differential equation defined over domain    , we assume its solution to be of the form        N
      ∑
u (x ) =   qj   j (x)
       j=1   .

The function j (x)   is called the jth  basis function. qj   j (x)   is called a trial function.

The function j (x)   is made up of functions called the shape functions Nk  as they are normally called in structural mechanics books.

qj   are the unknown coefficients which are determined by solving N set of equations generated by setting N integrals of the form ∫
   R (x ) v (x) dx
          i
 to zero. Where R (x )   is the differential equation residual and vi (x )   is the ith  weight function where i = 1 ⋅⋅⋅N .  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:

We divide the domain itself into a number of elements. Next, the j (x)   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 qi  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. Then qj  j (x)   will become the interpolation function.

In addition, the coefficients qj  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 du − u (x ) = 0
dx   with the given boundary condition of u (0) = 1   and defined over 0 ≤ x ≤ 1

pict

Using linear interpolation, then the solution u (x) ,  when x is located in first element, is found from

u (x) = q1 +slope (x − x0)

But the slope of the linear interpolating line over the first element is qx2−−qx1
 2  1     , hence the above becomes

pict

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 , then we write

u1 = q1 (x2-−-x)-+ q2(x-−-x1)
          h           h

Again, the above is valid for x1 ≤ x ≤ x2

We now do the same for the second element

 2     (x3-−-x)-    (x-−-x2)
u = q2    h    + q3   h

The above is valid for x ≤ x ≤ x
 2        3

and finally for the 3rd element

 3     (x4 − x)     (x − x3)
u = q3 --------+ q4--------
          h           h

The above is valid for x2 ≤ x ≤ x3

This is now illustrated in the following diagram

pict

Since our goal is to express the global approximate solution u as a series sum of basis functions each multiplied by q
 j  , we now rewrite each of the uj  to allow this, as follows

pict

Again, the above is valid for x1 ≤ x ≤ x2    . Notice the use of the following notation: Since each element will have defined on it 2 shape functions, N1 (x )   and N2 (x)   , then we use a superscript to indicate the element number. Hence for element 1   , we will write   1 ( )
N 1 x   and   1
N 2(x)   .

We now do the same for the second element

pict

The above is valid for x2 ≤ x ≤ x3

and finally for the 3rd element

pict

The above is valid for x2 ≤ x ≤ x3

The global trial function is

pict

The shape function for node 1 is

pict

And the shape function for node 2 is

pict

The shape function for node 3 is

pict

The shape function for the last node is

pict

Hence for the first node, the shape function is 1 = (x2−x),
     h  for the last node n = (x−xn−1)
      h  and the shape function for any internal node is

   x-− xj−1  xj+1 −-x
j =    h    +    h

The approximate solution is

        Numbe∑r Nodes
u (x) =           qi  i
            i=1

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.

pict

The diagram below illustrates the numbering used.

pict

Next the residual R (x)   is found by substituting the above global solution into the original differential equation.

The differential equation to solved is

adu-+ bu (x ) = f (x )
 dx

Defined over 0 ≤ x ≤ 1   with the initial condition u (0) = u0    .

N is the number of nodes, therefore the residual is

pict

The test functions are

vj =   j    j = 1 ⋅⋅⋅N

The weak form of the differential equation is

pict

N equations are obtained from the above, which are solved for the qi

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




i i   ′
i



1   (x2−x)
  h  −1
 h



2   {            }
 -(x−hx1), (x3−hx)   {    }
 1h,−h1



3   (x−x2)
  h  1
h



In general, if there are N nodes, then for the first node

pict

And for the last node

pict

And for any node in the middle

pict

for xi−1 ≤ x ≤ xi  and

pict

for x ≤  x ≤ x
 i        i+1    .

Looking back at the weak form integral above, it is evaluated as follows

x=∫xN   ⌊ N              ⌋
      ⟨||∑     (  ′     )||        ⟩
      j ||   qi a   i + b i || − f (x ) dx = 0     j = 1 ⋅⋅⋅N
x=x1   ∕⌈ i=1             ⌉        \

For the first node only, j = 1   , the following results

 x∫=xN    ⌊|N∑               ⌋|
       1⟨ ||   qi (a ′i +b   i)||−  f (x)⟩ dx = 0
       ∕ |⌈i=1              |⌉       \
x=x1

Since 1    is none zero only over x1 ≤ x ≤ x2    , and given that   x −x              (x−x)
1 =-2h-,  ′1 = −1h , 2 = --h-1  due to the range of integration limits, and that ′2 = 1h  , then the above can simplified to give

pict

The above simplifies further to

pict

Since x2 − x1 = h and x1 − x2 = −h the above reduces

pict

The above equation gives the first row in the global stiffness matrix for any first order linear ODE of the form   du
a dx + bu (x ) = f (x )   . The above shows that numerical integration is only needed to be performed on the term  x∫=x2
     (x2 − x) f (x)

x=x1

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

  x=xN
  ∫        ⌊|∑N   (  ′     )⌋|
         N ⟨||   qi a  i + b  i|| − f (x)⟩dx =  0
x=xN−1    ∕|⌈i=1            |⌉        \

Since N  domain of influence is xN −1 ≤ x ≤ xN  , the above simplifies to

  x∫=xN
         (^     (  ′           )      (  ′       )^        )
        N   qN −1 a  N −1 + b  N−1  +qN  a   N + b  N    − f (x ) dx = 0
x=xN−1

Since N−1 = (xN-−x),  ′   = −1,  N =  x−xN-−1,  ′  = 1
       h     N−1   h          h    N    h  the above becomes

 x=∫xN          ([     (               )      (               )]        )
      x-− xN-−1  qN−1  −a-+ b (xN--− x-) + qN  a1-+ bx-−-xN−1-   − f (x ) dx = 0
         h              h        h             h        h
x=xN−1

Which simplifies to

                                                                       x∫=xN
     (3a-+b-(xN-−1 − xN-)) (xN-−1-− xN-)2 (−3a+-2b-(xN−1 − xN-)) (xN-−1-− xN-)2 1
−qN−1            6h2             − qN             6h2              − h     (x − xN−1) f (x)dx = 0
                                                                      x=xN−1

Letting x    − x   = −h
  N−1   N in the above becomes

     (       )                     x∫=xN
q     bh-−-3a  + q  (3a +-2bh)-−-1      (x − x    ) f (x) = 0
 N−1     6        N     6       h            N −1
                                 x=xN−1

Hence the last line of the stiffness matrix can be determined directly except for the term under the integral which needs to be evaluated using numerical 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

  x=xj                                    x=xj+1
 ∫       ⌊|∑N    (        )⌋|                ∫      ⌊|∑N   (         )⌋|
       j ⟨ || qi  a  ′i + b i || − f (x)⟩dx +         j ⟨||  qi a  i′+b   i || − f (x)⟩ dx = 0
x=x    ∕ |⌈i=1             |⌉        \      x=x    ∕|⌈i=1             |⌉       \
   j−1                                       j

Where the integral was broken into two parts to handle the domain of influence of the shape functions.

   ∫x=xj
          (     (  ′         )     (   ′     )        )
         j   qj−1 a  j−1 + b  j−1  + qj a   j + b j −  f (x)  dx +
  x=xj−1
 x=xj+1
 ∫      (   (         )       (            )        )
       j   qj a   ′j + b j  +qj+1 a   ′j+1 + b  j+1 −  f (x)  dx   = 0
x=x
   j

For xj−1 ≤ x ≤ xj,  j = x−xj−1,  ′= 1,   j−1 = xj−x,  ′  = −1
    h     j  h         h    j− 1   h

For xj ≤ x ≤ xj+1,     xj+1−x-  ′  −1        x−xj  ′    1
j =  h   ,  j = h ,  j+1 =  h  ,  j+1 = h

Hence the weak form integral can be written as

   x∫=xj        (      (            )     (             )        )
       x-−-xj−1-  q     −a-+ bxj-−-x  + q  a-+ bx-−-xj−1  − f (x )   dx+
          h        j−1  h       h        j h       h
 x=xj−1
 x=∫xj+1       (    (              )       (           )        )
     xj+1 −-x      −a-    xj+1-−-x-        a-   x-−-xj
        h      qj   h + b    h     + qj+1 h + b  h     − f (x )   dx  = 0
x=xj

Which simplifies to

     (       (      )) (      )2     ((3a− 2b (x  − x)) (x  − x )2)     x∫=xj
 q  ⟨∕-−3a−-b-xj−1 − xj-xj−1 − xj-⟩\+ q----------j−1---j----j−1---j---− 1     (x − x  ) f (x)dx+
  j−1             6h2                j             6h2               h          j−1
   (∕      (       )) (      )  \     (     (       )) (      )     x=xj−1
    −3a− 2b xj − xj+1 xj − xj+1 2      3a− b xj − xj+1  xj − xj+1 2  1 x∫=xj+1(      )
qj ⟨∕-------------2-------------⟩\+qj+1⟨∕-------------2------------⟩\− -      xj+1 − x f (x )dx = 0
  ∕            6h              \     ∕            6h             \  hx=xj

For equal distance between elements, xj − xj+1 = − h,xj−1 − xj = − h the above simplifies to

      (             )     (           2)      x∫=xj
       (−3a-+-bh)-h2      -(3a +-2bh)-h--  1-     (       )
  qj−1       6h2       + qj     6h2       − h       x − xj−1 f (x )dx+
                                            x=xj−1
                                            x=xj+1
   ((−3a + 2bh)h2 )       ((3a + bh) h2)  1  ∫   (       )
qj  -------2------  +qj+1  ------2----- − --      xj+1 − x f (x )dx = 0
         6h                    6h         h x=x
                                               j

Combining gives

    (         )     (    )       (       )       x∫=xj                     x=∫xj+1
     −-3a +-bh       2bh-         3a +-bh    1-⟨∕     (       )                (        )        ⟩\
qj−1     6      +qj   3   + qj+1     6     − h ∕      x − xj−1  f (x )dx +       xj+1 − x f (x)dx \=  0
                                               ∕x=xj−1                     x=xj                   \

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

⌊2bh−3a  bh+3a                                      ⌋ ⌊    ⌋
||---6--   -6---    0       0      0      ...     0  || || q1 ||
||−3a+bh    2bh-   3a+bh     0      0      ...     0  || || q2 ||
||   6      3       6                                || |||    |||
||  0     −3a6+bh-   2b3h    3a+6bh-    0      ...     0  || || q3 ||
||                −3a+bh-   2bh    3a+bh-              || ||  .. ||
||  0       0       6       3       6     ...     0  || ||  . || =
||   ...      ...       ...       ...       ...     ...     ...  || ||  ... ||
||                                                   || ||  . ||
||  0       0       0     −3a+6bh-   2bh3-   3a+6bh-    0  || ||  .. ||
||  0       0       0       0    −3a+bh   2bh   3a+bh|| ||    ||
|||                                  6      3      6  ||| ||qN−1||
|⌈  0       0       0       0      0     bh−63a  3a+26bh|⌉ || qN ||
                        x=x                           ⌈    ⌉
    ⌊|                   ∫  2                             ⌋|
    ||                 1h    (x2 − x) f (x)                ||
    |||                  x=x                               |||
    ||       x∫=x2          1       x∫=x3                   ||
    ||    1⟨                                          ⟩   ||
    ||    h∕∕    (x − x1) f (x)dx +    (x3 − x ) f (x)dx \\ ||
    ||      x=x1                  x=x2                     ||
    ||     ∕ x∫=x3                  x∫=x4               \   ||
    ||    1⟨∕                                          ⟩\   ||
    ||    h∕    (x − x2) f (x)dx +    (x4 − x ) f (x)dx \ ||
    ||     ∕x=x2                  x=x3                 \   ||
    ||                          ..                         ||
    ||                          .                         ||
    ||    x∫=xj(        )           x∫=xj+1(       )         ||
    ||1 ⟨∕∕      x − xj−1 f (x)dx +       xj+1 − x  f (x) dx⟩\\||
    ||h                                                   ||
    ||  ∕x=xj−1                     x=xj                   \||
    |||                          ...                         |||
    ||   x=xN−1                     x=xN                  ||
    ||  ⟨ ∫                         ∫                    ⟩||
    ||1h ∕∕      (x − xN −2) f (x)dx +     (xN − x) f (x) dx\\||
    ||  x=x                        x=x                     ||
    ||  ∕  N−2          x=∫xN          N−1                 \||
    ||               1                                    ||
    ||               h      (x − xN−1) f (x)              ||
    |⌈                x=xN−1                               |⌉

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 qi  are solved for. But before doing that, q1    is first replaced by the initial condition given in the problem, and the first row is removed as follows

⌊2bh−3a  bh+3a                                      ⌋ ⌊    ⌋
||---6--   -6---    0       0      0      ...     0  || || u0 ||
||−3a+bh    2bh-   3a+bh     0      0      ...     0  || || q2 ||
||   6      3       6                                || |||    |||
||  0     −3a6+bh-   2b3h    3a+6bh-    0      ...     0  || || q3 ||
||                −3a+bh-   2bh    3a+bh-              || ||  .. ||
||  0       0       6       3       6     ...     0  || ||  . || =
||   ...      ...       ...       ...       ...     ...     ...  || ||  ... ||
||                                                   || ||  . ||
||  0       0       0     −3a+6bh-   2bh3-   3a+6bh-    0  || ||  .. ||
||  0       0       0       0    −3a+bh   2bh   3a+bh|| ||    ||
|||                                  6      3      6  ||| ||qN−1||
|⌈  0       0       0       0      0     bh−63a  3a+26bh|⌉ || qN ||
                        x=x                           ⌈    ⌉
    ⌊|                   ∫  2                             ⌋|
    ||                 1h    (x2 − x) f (x)                ||
    |||                  x=x                               |||
    ||       x∫=x2          1       x∫=x3                   ||
    ||    1⟨                                          ⟩   ||
    ||    h∕∕    (x − x1) f (x)dx +    (x3 − x ) f (x)dx \\ ||
    ||      x=x1                  x=x2                     ||
    ||     ∕ x∫=x3                  x∫=x4               \   ||
    ||    1⟨∕                                          ⟩\   ||
    ||    h∕    (x − x2) f (x)dx +    (x4 − x ) f (x)dx \ ||
    ||     ∕x=x2                  x=x3                 \   ||
    ||                          ..                         ||
    ||                          .                         ||
    ||    x∫=xj(        )           x∫=xj+1(       )         ||
    ||1 ⟨∕∕      x − xj−1 f (x)dx +       xj+1 − x  f (x) dx⟩\\||
    ||h                                                   ||
    ||  ∕x=xj−1                     x=xj                   \||
    |||                          ...                         |||
    ||   x=xN−1                     x=xN                  ||
    ||  ⟨ ∫                         ∫                    ⟩||
    ||1h ∕∕      (x − xN −2) f (x)dx +     (xN − x) f (x) dx\\||
    ||  x=x                        x=x                     ||
    ||  ∕  N−2          x=∫xN          N−1                 \||
    ||               1                                    ||
    ||               h      (x − xN−1) f (x)              ||
    |⌈                x=xN−1                               |⌉

Now the first row is removed (and remembering to multiply u0    by the first entry in the second row) gives

⌊−3a+bh      2bh    3a+bh                              ⌋ ⌊    ⌋
||---6--u0    3--   --6--     0      0      ...     0  || || q2 ||
||   0      −3a+bh   2bh    3a+bh-    0      ...     0  || || q3 ||
||            6       3       6                        || ||  . ||
||   0        0     −3a+6bh-   2b3h    3a+b6h-   ...     0  || ||  .. ||
||    ..       ..       ..       ..       ..     ..      ..  || ||  .. ||
||    .       .       .       .       .       .     .  || ||  . || =
||   0        0       0     −3a+bh-   2bh-   3a+bh-    0  || ||  .. ||
||                            6       3      6         || ||  . ||
|||   0        0       0       0    −3a+6bh   2b3h   3a+6bh||| |||qN−1|||
||   0        0       0       0      0     bh−3a  3a+2bh|| ||    ||
⌈                                           6      6  ⌉ ⌈ qN ⌉
     ⌊       x∫=x2                  x∫=x3                   ⌋
     ||    1⟨∕                                          ⟩\   ||
     ||    h∕    (x − x1) f (x)dx +    (x3 − x ) f (x)dx \ ||
     |||     ∕x=x1                  x=x2                 \   |||
     ||       x∫=x3                  x∫=x4                   ||
     ||    1⟨∕    (x − x ) f (x)dx +    (x − x ) f (x)dx ⟩\  ||
     ||    h∕          2                 4             \   ||
     ||     ∕x=x2                  x=x3                 \   ||
     ||                          ...                         ||
     ||    x=xj                     x=xj+1                  ||
     ||  ⟨ ∫   (        )           ∫    (       )        ⟩||
     ||1h ∕∕      x − xj−1 f (x)dx +       xj+1 − x  f (x) dx\\||
     ||  x=x                       x=x                     ||
     ||  ∕  j−1                   .    j                   \||
     ||                          ..                         ||
     ||   x=∫xN−1                     x∫=xN                  ||
     ||1 ⟨∕                                                ⟩\||
     ||h ∕      (x − xN −2) f (x)dx +     (xN − x) f (x) dx\||
     |||  ∕x=xN−2                     x=xN−1                 \|||
     ||                  x=∫xN                               ||
     ||               1      (x − x   ) f (x)              ||
     ||               h            N−1                     ||
     ⌈                x=xN−1                               ⌉

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

   ⌊  2bh    3a+bh                              ⌋ ⌊    ⌋
   ||  3--   --6--     0       0     ...     0  || || q2 ||
   ||−3a+bh   2bh    3a+bh     0     ...     0  || || q3 ||
   ||  6       3       6                        || ||  . ||
   ||  0     −3a6+bh-   2b3h    3a+6bh-   ...     0  || ||  .. ||
   ||   ..      ..       ..       ..     ..      ..  || ||  .. ||
   ||   .      .       .       .       .     .  || ||  . || =
   ||  0       0     −3a+bh-   2bh-   3a+bh    0  || ||  .. ||
   ||                  6       3      6         || ||  . ||
   |||  0       0       0    −3a+6bh   2b3h   3a+6bh||| |||qN−1|||
   ||  0       0       0       0    bh−3a  3a+2bh|| ||    ||
   ⌈                                 6      6  ⌉ ⌈ qN ⌉
⌊   ∫x=x2                  x∫=x3                             ⌋
|| 1⟨∕                                         ⟩\   (−3a+bh  )||
||h ∕    (x − x1) f (x )dx +    (x3 − x) f (x )dx\ −   6   u0 ||
|||  ∕x=x1                  x=x2                \             |||
||          x∫=x3                  x∫=x4                      ||
||       1⟨∕    (x − x ) f (x)dx +    (x − x ) f (x)dx ⟩\     ||
||       h∕          2                 4             \      ||
||        ∕x=x2                  x=x3                 \      ||
||                             ...                            ||
||       x=xj                     x=xj+1                     ||
||     ⟨ ∫   (        )           ∫    (       )        ⟩   ||
||   1h ∕∕      x − xj−1 f (x)dx +       xj+1 − x  f (x) dx\\  ||
||     x=x                       x=x                        ||
||     ∕  j−1                   .    j                   \   ||
||                             ..                            ||
||      x=∫xN−1                     x∫=xN                     ||
||   1 ⟨∕                                                ⟩\   ||
||   h ∕      (x − xN −2) f (x)dx +     (xN − x) f (x) dx\  ||
|||     ∕x=xN−2                     x=xN−1                 \   |||
||                     x=∫xN                                  ||
||                  1      (x − x   ) f (x)                 ||
||                  h            N−1                        ||
⌈                   x=xN−1                                  ⌉

Now the system above is solved for q2,q3, ⋅⋅⋅,qN  . Once the qi  are found, the solution is

        Numbe∑r Nodes
u (x) =           q
                    ii
            i=1

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

?? ?? ?? ??

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 coefficients of the form

  2
ad-u-+ bdu- +cu (x) = f (x )
 dx 2   dx

with the initial conditions u (x0) = u0    (called essential Dirichlet condition), and the Neumann boundary condition du =
dx at L , for x0 ≤ x ≤ L

As before, the solution is assumed to be of the form

    ∑N
u =    qi  i
    i=1

And the unknown q are found by solving N algebraic equations

∫ L
      R (x )dx = 0     j = 1 ⋅⋅⋅N
    j
x0

Where R (x )   is the ODE residual obtained by substituting the assumed solution into the original ODE. Hence the above becomes (For simplicity, j = 1⋅⋅⋅N is not written each time as it is assumed to be the case)

pict

Applying integration by parts on ∫L   jadd2xu2dx
x0 . Since

    (      )
-d-     du--      d2u-     ′du-
dx  a    jdx   = a   jdx 2 + a  jdx

then integrating both sides of the above gives

pict

Hence

∫ L            [      ]    ∫ L
      d2u-         du- L        ′du-
   a  jdx 2dx =  a  jdx    −    a  jdx dx
x0                     x0  x0

Substituting the above in equation A1 gives

              L             L
[    du] L   ∫     du     ∫    ( du                 )
 a   j---   −   a   ′j--dx +      j b ---+ cu (x) − f (x ) dx = 0
     dx  x0   x0    dx      x0     dx
(A2)

Consider the term

[      ]L          |           |
 a  jdu-   =  a  jdu-||   − a  jdu-||
    dx  x0      dx ||x=L      dx ||x=x0

First consider the term a    du-||
   jdx |x=L  . Since du =
dx at L then this term becomes a   ^   ^
     j x=L  . But ^   ^
   jx=L  is non zero only for the Nth  shape function evaluated at x = L . Since linear interpolation is used, the   N th  shape function is N = x−xhN−1  which have the value of 1   at x = xN  . Hence

    du||
a   jdx|||    = a      whenj  = N otherwise0
       x=L

Now consider the term    du|
a  jdx||x=x0     , since at x = x0    all shape functions will be zero except for 1    which has the value of 1   at x = x0    . Hence this simplifies to a du-||
  dx |x=x0

recalling that du-   d-N∑        ∑N    ′
dx = dx   qi  i =    qi  i
       i=1       i=1  .

But at x = x0    only 1    is defined and its has the slope of −h1  , hence

    |
adu-||   =  −aq1-when  j = 1, otherwise 0
 dx ||x=x0    h

However, q  = u
 1    0    since that is by definition the initial condition. Finally one obtains

[     ] L   (|  aq1    j = 1 )|
a   jdu-   = {   h           }
    dx  x0   |( a       j = N  |)

Hence the symmetric weak form equation A2 can now be simplified more giving

(|  aq1-    j = 1 )|   ∫L           ∫L   (                    )
{  h            } −   a   ′jdu-dx +      j bdu- +cu (x) − f (x ) dx = 0     j = 1 ⋅⋅⋅N
|( a       j = N  |)        dx             dx
                    x0            x0
(A3)

The trial function obtain by linear interpolation are used. These are shown again in this table

|-----------|------------------|--------|
|i          |                   |  ′      |
|-----------|-i----------------|-i------|
|1          |(x2−hx)             |−1h      |
|-----------|{-(x−xi−1)-(xi+1−x)}--|{1-−-1}-|
|1-<-i-<-N--|----h---,---h-----|-h,-h---|
|N          |(x−xN−1)           |1       |
----------------h---------------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

       L             L
aq    ∫     du     ∫     ( du                )
--1-−   a   ′1--dx +      1  b---+ cu (x ) − f (x) dx = 0
 h          dx             dx
      x0            x0

Since     ∑N
u =    q
        i i
     i=1  , and since 
1    domain of influence is only from x
 1    to x
 2    then above becomes

     ∫ x2                        ∫x2   (                                        )
aq1-      ′ d--                         -d-
h  −    a  1dx  (q1   1 + q2 2)dx +      1  bdx (q1  1 +q2   2) + c (q1 1 + q2 2) − f (x) dx = 0
     x1                         x1

But    (x2−x)
1 =   h  and ′   −1
 1 = h  and over the domain from x1    to x2,     (x−x1)  ′   1
2 =   h  ,  2 = h  , hence the above becomes

      ∫x2                   ∫x2        (                 (                      )        )
aq1-     −a-(−q1-  q2 )        (x2-− x-)  (−-q1  q2-)       (x2-− x-)    (x-−-x1)-
 h  −    h    h  +  h  dx +       h     b   h  +  h   +c  q1   h    + q2    h     − f (x)  dx = 0
      x1                     x1

The above simplifies to

pict

The above gives the first row in the stiffness matrix. Since it is assumed that each element will have the same length, hence x1 − x2 = −h,  and the above becomes

pict

For the last equation, which will be the last row in the global stiffness matrix

     ∫xN            ∫xN   (                    )
           ′du-              du-
a   −    a   Ndx dx +      N  b dx + cu (x ) − f (x) dx = 0
    xN−1           xN−1

Since     ∑N
u =    qi  i
     i=1  , and since N  domain of influence is only from xN −1    to xN  The above becomes

    ∫xN                              ∫xN    (                                                    )
          ′-d-                                -d-
a  −    a   Ndx (qN −1  N−1 +qN   N )dx +      N  bdx (qN−1   N−1 + qN  N ) + c (qN −1 N−1 + qN N ) − f (x) dx = 0
   xN−1                             xN−1

But    (x−xN−1)
N = ---h---  and ′    1
 N = h  and over the domain from xN−1    to xN ,       (xN−x)  ′     −1
N−1 = --h---,  N−1 = h-  hence the above becomes

                                ∫xN   (     (   )       )
                                   a-       −1-      1-
                          a   −     h  qN−1   h  + qN h  dx +
                               xN−1
 ∫xN            ( (     (   )       )    (                             )       )
    (x-−-xN−1)          −-1       1-          (xN-−-x)-     (x-−-xN-−1)
        h       b qN −1  h   + qN h  + c qN −1   h     + qN     h      −  f (x) dx = 0
xN−1

The above simplifies to

                   (                                              )
              qN −1           b (xN −1 − xN )2 c (xN−1 − xN)3
         a   − --2-- − axN−1 − --------------−  --------------+ axN  +
               h                    2                6
q  (          b (x    − x  )2   c (x   − x )3       )   1∫xN
-N- +axN −1 + ----N−1---N---−  ---N−1----N---− axN  − --    (x − xN −1) f (x) = 0
h2                  2                3                h
                                                       xN−1

The above represents the last row in the global stiffness matrix. Since it is assumed that each element will have the same length, hence xN −1 − xN = −h,  and the above becomes

pict

The expression for any row in between the first and the last rows is now found. For a general node j the result is

  xj+1           xj+1
  ∫    du      ∫    (  du                )
−   a   ′j--dx +      j b ---+ cu (x) − f (x ) dx = 0
 x     dx     x       dx
  j−1           j−1

Breaking the integral into halves to make it easier to write the trial functions over the domain of influence gives

    xj           xj+1            xj                              xj+1
  ⟨∫     du     ∫     du   ⟩   ∫    ( du                 )     ∫    (  du                )
− ∕∕   a  ′j--dx +    a   ′j--dx\\ +      j b ---+ cu (x) − f (x ) dx +      j b ---+ cu (x) − f (x ) dx = 0
  x      dx     x     dx      x       dx                       x       dx
  ∕j−1            j         \   j−1                               j

Considering the first domain xj−1 ≤ x ≤ xj  gives

 ∫ xj           ∫xj  (                    )
−   a   ′du-dx +      j bdu--+cu (x) − f (x ) dx
       jdx             dx
 xj−1          xj−1

Over this range, u = qj−1  j−1 + qj  j  hence du = qj−1   ′j−1 +qj  ′j
dx  where j−1 = xj−x, ′j−1 = −1,  j = x−xj−1, ′j = 1
      h         h        h        h  hence the above becomes

  ∫xj                         ∫xj
      1 (    (−1 )    1 )        x − xj−1(  (    (− 1)     1)    (    xj − x    x − xj−1)       )
−   a --qj−1  --- + qj-- dx +    -------- b  qj−1 ---  + qj-- + c qj−1------+ qj--------  − f (x ) dx
 xj−1  h       h       h       xj−1   h              h       h            h          h
(A4)

Now considering the second domain xj ≤ x ≤ xj+1    gives

 ∫xj+1          x∫j+1 (                    )
−   a   ′du-dx +        bdu--+cu (x) − f (x ) dx
       jdx          j  dx
 xj            xj

Over this range, u = qj  j + qj+1 j+1    hence du = qj  ′j + qj+1  ′j+1
dx    where j = xj+1−x-,  ′j = −1,  j+1 = x−xj,  ′j+1 = 1
    h         h         h         h  hence the above becomes

  x∫j+1                        x∫j+1
     (−a ) (−qj   qj+1)         xj+1 − x ( (−qj       1)    (  xj+1 − x      x − xj)       )
−     ---   ----+ ---- dx +     -------- b  ----+ qj+1-- + c qj--------+ qj+1------  − f (x) dx
  xj    h     h     h         xj    h         h        h           h            h
(A5)

Combine A4 and A5 and simplify we obtain

                                     (        )2    (        )3
                     qj−1⟨         b--xj−-1 −-xj-   c-xj−1 −-xj--     ⟩
                      h2 ∕− axj−1 −      2      −       6       + axj\
                         ∕                                           \
                    (        )2    (        )3     (       ) 2    (        )3
       qj-⟨       b--xj−1-−-xj--  c--xj−1 −-xj--  b-xj-−-xj+1---  c--xj −-xj+1-       ⟩
     + h2 ∕axj− 1 +      2      −       3       −       2      −       3      − axj+1\
          ∕                                                                         \
              b (x − x   )2  c (x − x   )3           ∫xj                  x∫j+1
+qj+1-⟨∕−ax  + ----j---j+1--−  ---j---j+1---+ ax  ⟩\ −    x-−-xj−1-f (x) dx −   xj+1-−-xdx = 0
  h2       j        2              6           j+1          h                    h
      ∕                                          \  xj−1                   xj

Since we are assuming each element will have the same length, hence xj−1 − xj = − h,  xj − xj+1 = − h and the above becomes

                           (  axj−1   b   ch   axj)
                       qj−1 − ---2- − --+ ---+ -2-
                             ( h      2   6    h)
                              axj−1  2ch   axj+1
                         +qj  --2--+ ----− ---2-
                               h      3      h
      (−ax             ax   )   ∫xjx − x             x∫j+1x   − x
+qj+1  ---j-+ b-+ ch-+ --j+1  −    -----j−1-f (x) dx −   -j+1----f (x)dx = 0
        h2    2    6    h2            h                    h
                               xj−1                   xj

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

 (                      )  (                  )
⌊| − a+ ax1-− b+ ch-− ax2    − ax1+ b − ch+ ax2-              0                      0          ⌋|
||  (h   h2   2   3   h2)       h(2   2   6    h)2     (                  )                        ||
||   − axh12-− b2 + ch6-+ ahx22         ahx12-+ 2ch3-− axh32         −hax22-+ b2 + ch6-+ ahx23              0          ||
||                          ( axj−1   b   ch-  axj)    (axj−1   2ch-  axj+1)     (−axj-  b   ch-  axj+1)||
||           0               − h2  − 2 + 6 +  h2        h2  +  3 −  h2        h2  + 2 + 6 +  h2  ||
||           0                        0                      ...                    ⋅⋅⋅         ||
||                                                 (  ax      b  ch   ax )  ( ax      b  ch   ax )||
|⌈           0                        0             − -Nh−21-− 2 + 6-+ -hN2     -Nh−21-+ 2 +-3 − -h2N |⌉
                             ⌊                ∫x2                           ⌋
                             ||              -1                              ||
                             ||              h    (x2 − x) f (x )             ||
                             ||                x1                             ||
                   ⌊| q1 ⌋|    ||                       ..                      ||
                   ||    ||    ||                       .                      ||
                   || q2 ||    ||                       ...                      ||
                   ||  .. ||    ||                                              ||
                   |||  . |||    |||                       ...                      |||
                   ||  ... ||    || ∫xj                     ∫xj+1                 ||
                   ||    || =  ||1   (       )           1    (       )        ||
                   || qj ||    ||h    x − xj−1  f (x)dx + h     xj+1 − x  f (x )dx||
                   ||  .. ||    ||xj−1                      xj                   ||
                   ||  . ||    ||                       ..                      ||
                   ||qN− 1||    ||                       .                      ||
                   ||    ||    ||                       ..                      ||
                   |⌈ qN |⌉    ||             x         .                      ||
                             ||            ∫ N                               ||
                             ||          1h    (x − xN−1) f (x ) − a            ||
                             ||           x                                  ||
                             ⌈            N−1                               ⌉

we must first replace q1    by the initial condition given in the problem, and we must remove the first row after that as follows

⌊( a   ax1  b   ch   ax2)  (  ax1  b   ch  ax2)                                                ⌋
|| −h( + h2 − 2 +  3 − h)2    − h(2 + 2 − 6 +  h)2     (        0        )              0          ||
||   −ax12 − b-+ ch + ax22        ax12-+ 2ch-− ax32        −ax22+ b + ch+ ax32               0          ||
|||     h    2   6   h          h     3    h          h    2.  6    h                           |||
||           0                       0                      ..                     ⋅⋅⋅         ||
||                                                         ..                                  ||
||           0                       0           (           .         )  (        ⋅⋅⋅       ) ||
||           0                       0            − axN−1− b + ch+ axN     axN−1+ b + ch-− axN  ||
⌈                                                   h2    2   6    h2      h2    2   3    h2  ⌉
                            ⌊|                ∫x2                           ⌋|
                            ||              -1   (x2 − x) f (x )             ||
                            ||              h                               ||
                  ⌊    ⌋    ||                x1                             ||
                  || u0 ||    ||                       ...                      ||
                  || q2 ||    ||                       .                      ||
                  ||  . ||    ||                       ..                      ||
                  ||  .. ||    ||                       ..                      ||
                  |||  .. |||    |||                       .                      |||
                  ||  . || =  || ∫xj(       )            ∫xj+1(       )        ||
                  || q  ||    ||1    x − xj−1  f (x)dx + 1     xj+1 − x  f (x )dx||
                  ||  j. ||    ||h                       h                     ||
                  ||  .. ||    ||xj−1                      xj                   ||
                  ||    ||    ||                       ...                      ||
                  ||qN− 1||    ||                       .                      ||
                  || qN ||    ||                       ..                      ||
                  ⌈    ⌉    ||            ∫xN                               ||
                            ||          1                                   ||
                            ||          h    (x − xN−1) f (x ) − a            ||
                            |⌈           xN−1                               |⌉

Multiply the first element in the second row by u0    and remove the first row we get

⌊(  ax1  b   ch  ax2)       (ax1  2ch   ax3)      (−ax2   b  ch   ax3)                         ⌋
|| − h2 − 2 + 6 +  h2  u0  (  h2 +  3 −  h2  )      h(2 +  2 + 6 + h)2      (        0        ) ||
||           0              −ax22 − b+ ch + ax32-       ax22 + 2ch-− ax42-        −ax23+ b + ch+ ax42  ||
||                            h    2   6   h          h    3.    h           h    2.  6    h   ||
||           0                      0                      ..                      ..         ||
||                                               ( axN−1  b   ch   axN-)  (axN−1  b   ch   axN-)||
|⌈           0                      0             −  h2 − 2 +  6 + h2       h2  + 2 +  3 − h2  |⌉
                            ⌊      ∫x2                   ∫x3               ⌋
                            |||    -1   (x − x ) f (x )dx + 1  (x  − x )dx    |||
                            ||    h         2            h     3            ||
                  ⌊    ⌋    ||      x1                     x2                ||
                  || q2 ||    ||                       ...                      ||
                  ||  .. ||    ||                       .                      ||
                  ||  . ||    ||                       ..                      ||
                  ||  ... ||    || ∫xj                     ∫xj+1                 ||
                  ||    ||    ||1   (       )           1    (       )        ||
                  || qj || =  ||h    x − xj−1  f (x)dx + h     xj+1 − x  f (x )dx||
                  ||  .. ||    ||xj−1                      xj                   ||
                  ||  . ||    ||                       ..                      ||
                  ||qN− 1||    ||                       .                      ||
                  ||    ||    ||                       ...                      ||
                  |⌈ qN |⌉    ||             xN                               ||
                            ||            ∫                                 ||
                            |||          1h    (x − xN−1) f (x ) − a            |||
                            |⌈           x                                  |⌉
                                         N−1

Now move the first element of the first row above to the RHS and remove the first column we obtain

⌊(ax1  2ch  ax3)   (−ax2-  b   ch-  ax3)                                                ⌋
|| h2 +  3 −  h2   (  h2  + 2 + 6 +  h2 )      (      0      )      (         0        ) ||
||       0          − ax22− b + ch+ ax32        ax22+ 2ch− ax42        −ax2j+ b + ch+ axj2+1  ||
||       .            h    2.  6    h         h     3    h          h    2 . 6    h    ||
||       ..                   ..                    ⋅⋅⋅                      ..         ||
||                                       (  axN−1-  b   ch-  axN)   (axN−1-  b   ch-  axN) ||
|⌈       0                   0            −  h2  − 2 + 6 +  h2      h2  + 2 + 3 −  h2  |⌉
                ⌊ ∫x2                    ∫x3             (                 )   ⌋
                |||1   (x − x ) f (x)dx + 1   (x − x)dx  −  −ax1− b + ch + ax2- u |||
                ||h         2           h      3             h2   2    6   h2   0||
       ⌊    ⌋   || x1                     x2                                     ||
       || q2 ||   ||                               ...                              ||
       || ..  ||   ||                               .                              ||
       || .  ||   ||                               ..                              ||
       || ...  ||   ||         ∫xj                     ∫xj+1                         ||
       ||    ||   ||        1   (       )           1    (       )                ||
       || qj || = ||        h    x − xj−1  f (x)dx + h     xj+1 − x  f (x )dx        ||
       || ..  ||   ||        xj−1                      xj                           ||
       || .  ||   ||                               ..                              ||
       ||qN −1||   ||                               .                              ||
       ||    ||   ||                               ...                              ||
       |⌈qN  |⌉   ||                     xN                                       ||
                ||                    ∫                                         ||
                |||                  1h    (x − xN−1) f (x ) − a                    |||
                |⌈                   x                                          |⌉
                                     N−1

Now we solve for the vector [q2,q3,⋅⋅⋅,qN ]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.

animation

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 fluid dynamics, Volume I. C.A.J.Fletcher. Springer-Verlag