EduPanda Logo

Numerical methods

Finite Difference Method (FDM)

  1. Introduction to the finite difference method - basic information
  2. Central difference formulas for the one-dimensional problem
  3. Solving differential equations using the finite difference method
    • Example 1 - boundary problem
  4. Solving bending beams using the finite difference method
  5. Boundary conditions for beams
    • Fixed support
    • Hinged support
    • Sliding support
    • Free end
  6. Example 2 - solving a beam
    • Bonus - calculating shear force

Introduction to the finite difference method

A method that approximates the derivative of a function using appropriate difference formulas. In the following study, we will focus on two simple applications of the finite difference method:

  1. Solving differential equations
  2. Solving bending beams

We replace individual derivatives with difference formulas using the formulas below.

Central difference formulas for the one-dimensional problem

Metoda różnic skończonych - Centralne wzory różnicowe dla zagadnienia jednowymiarowego

\[ \begin{aligned} &f^I=\frac{1}{2 h}\left(-v_{i-1}+v_{i+1}\right) \\ &f^{II}=\frac{1}{h^2}\left(v_{i-1}-2 v_i+v_{i+1}\right) \\ &f^{III}=\frac{1}{2 h^3}\left(-v_{i-2}+2 v_{i-1}-2 v_{i+1}+v_{i+2}\right) \\ &f^{IV}=\frac{1}{h^4}\left(v_{i-2}-4 v_{i-1}+6 v_i-4 v_{i+1}+v_{i+2}\right) \end{aligned} \]

where \( h \) - length of the element.

Solving differential equations using the finite difference method

The general procedure in this case:

  1. Dividing the investigated interval into n elements of length h
  2. Replacing derivatives in the equation with appropriate difference formulas
  3. Writing difference equations for individual points
  4. Considering boundary conditions
  5. Solving the system of equations in matrix form

Example 1

Content

Solve the boundary problem with the data as below

\[ \begin{aligned} &y''(x)=2 x \quad x \in[a, b] \\ &a=3 \quad b=7 \\ &y(a)=4 \\ &y'(b)=2 \end{aligned} \]

Dividing the interval into \( n=4 \) elements

Solution

Length of a single element:

\[ h=\frac{b-a}{n}=1 \]

Position of individual nodes

\[ x_i=x_0+i\cdot h \]

We replace the equation \(y''(x)=2x\) using the difference formulas:

\[ \begin{aligned} &\frac{1}{h^2}\left(i_{-1}-2 \cdot i+i_{p1}\right)=2 x \\ &i_{-1}-2 \cdot i+i_{p1}=h^2 \cdot 2 x \end{aligned} \]

For each node \(i=1,2,3,4\) we write the difference equation according to the above formula, additionally recording the values of x for individual nodes:

\[ \begin{array}{lll} i=1 & i_0-2 \cdot i_1+i_2=h^2 \cdot 2 x_1 & x_1=x_0+h=4 \\ i=2 & i_1-2 \cdot i_2+i_3=h^2 \cdot 2 x_2 & x_2=x_0+2 h=5 \\ i=3 & i_2-2 \cdot i_3+i_4=h^2 \cdot 2 x_3 & x_3=x_0+3 h=6 \\ i=4 & i_3-2 \cdot i_4+i_5=h^2 \cdot 2 x_4 & x_4=x_0+4 h=7 \end{array} \]

We consider the boundary conditions, also replacing them with difference formulas

\[ \begin{array}{lll} y(a)=4 & y_0=4 & \\ y'(b)=2 & \frac{1}{2 h}\left(-i_3+i_5\right)=2 & -i_3+i_5=4 h \end{array} \]

We write the above equations in matrix form \(A\cdot y = B\)

\[ \left[\begin{array}{cccccc} 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 1 \end{array}\right] \cdot\left[\begin{array}{l} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \end{array}\right]=\left[\begin{array}{c} h^2 \cdot 2 \cdot x_1 \\ h^2 \cdot 2 \cdot x_2 \\ h^2 \cdot 2 \cdot x_3 \\ h^2 \cdot 2 \cdot x_4 \\ 4 \\ 4 \cdot h \end{array}\right]=\left[\begin{array}{c} 8 \\ 10 \\ 12 \\ 14 \\ 4 \\ 4 \end{array}\right] \]

And finally the solution (in this case obtained in Mathcad):

\[ y=A^{-1} B=\left[\begin{array}{r} 4 \\ -31 \\ -58 \\ -75 \\ -80 \\ -71 \end{array}\right] \]

In practice, we are only interested in the first 5 results, as the last one is actually outside the investigated interval \(i=5 \to x=8\)

The analytical solution of this equation obtained by classical methods of solving differential equations:

\[ y_{an}(x)=\frac{x^3}{3}-47 x+136 \]

The table summarizes the values along with errors for individual points

\[ \begin{array}{|r|r|r|r|r|r|} \hline \mathrm{i} & x & y & y_a & \text{Błąd bezwzględny} & \text{Błąd względny} \\ \hline 0 & 3 & 4 & 4 & 0 & 0,00 \% \\ \hline 1 & 4 & -31 & -30,667 & 0,333 & 1,09 \% \\ \hline 2 & 5 & -58 & -57,333 & 0,667 & 1,16 \% \\ \hline 3 & 6 & -75 & -74 & 1 & 1,35 \% \\ \hline 4 & 7 & -80 & -78,667 & 1,333 & 1,69 \% \\ \hline \end{array} \]

If we divided the interval into more elements, the accuracy would of course be greater, but in such a case it would be necessary to use computational programs like Matlab

Solving bending beams using the finite difference method

The procedure is almost identical to that of solving differential equations, as we start from the relationship between displacement and load of a simple beam:

\[ \frac{d^4 v(x)}{d x^4}=\frac{q(x)}{E J} \]

Which in difference notation takes the form:

\[ \frac{d^4 v(x)}{d x^4}=\frac{p_y(x)}{E J} \Longrightarrow \frac{v_{i-2}-4 v_{i-1}+6 v_i-4 v_{i+1}+v_{i+2}}{h^4}=\frac{q_{i}}{E J} \]

Or in a simpler notation:

\[ v_{i-2}-4 v_{i-1}+6 v_i-4 v_{i+1}+v_{i+2}=b_i, \quad b_i=h^4 \cdot \frac{q_{i}}{E J} \]

Problems arise when the beam is not loaded continuously, but for example by applying a point force, trapezoidal load, etc., or when its stiffness varies along the length of the beam. In such cases, we will have to reduce the load to a continuous load.


Additionally, after calculating the deflection values of the beam at individual points, we will use known relationships that will allow us to calculate the shear force and bending moment:

\[ \begin{aligned} &M(x)=-E J \frac{d^2 v(x)}{d x^2} \Rightarrow \mathrm{M}=-\mathrm{EI} \cdot \frac{v_{i-1}-2 v_i+v_{i+1}}{h^2} \\ &T(x)=-E J \frac{d^3 v(x)}{d x^3} \Rightarrow T=-\mathrm{EI} \cdot \frac{-v_{i-2}+2 v_{i-1}-2 v_{i+1}+v_{i+2}}{2 h^3} \end{aligned} \]

Boundary conditions for beams

The boundary conditions for the beam will result from the way it is supported as shown below

metoda różnic skończonych - warunki brzegowe dla belki - utwierdzenie pełne \[ \begin{aligned} &v=0 \rightarrow v_{i}=0 \\ &\frac{d v}{d x}=0 \rightarrow \frac{-v_{i-1}+v_{i+1}}{2 h}=0 \end{aligned} \]
metoda różnic skończonych - warunki brzegowe dla belki - podpora przegubowa \[ \begin{aligned} &v=0 \rightarrow v_{i}=0\\ &M=-E J \frac{d^2 v}{d x^2}=0 \rightarrow-E J \frac{v_{i-1}-2 v_i+v_{i+1}}{h^2}=0 \end{aligned} \]
metoda różnic skończonych - warunki brzegowe dla belki - utwierdzenie przesuwne \[ \begin{aligned} &\frac{d v}{d x}=0 \rightarrow \frac{-v_{i-1}+v_{i+1}}{2 h}=0 \\ &T=-E J \frac{d^3 v}{d x^3}=0 \rightarrow-E J \frac{-v_{i-2}+2 v_{i-1}-2 v_{i+1}+v_{i+2}}{2 h^3}=0 \end{aligned} \]
metoda różnic skończonych - warunki brzegowe dla belki - swobodny koniec \[ \begin{aligned} &M=-E J \frac{d^2 v}{d x^2}=0 \rightarrow-E J \frac{v_{i-1}-2 v_i+v_{i+1}}{h^2}=0 \\ &T=-E J \frac{d^3 v}{d x^3}=0 \rightarrow-E J \frac{-v_{i-2}+2 v_{i-1}-2 v_{i+1}+v_{i+2}}{2 h^3}=0 \end{aligned} \]

Example 2

Content

Solve the given beam using the FDM, \(EI = 8000 kNm^2\)

Metoda różnic skończonych - przykład 2 - belka

Solution

We write the difference equation for the points i=1,2,3,4

\[ \begin{array}{ll} i=1 & 1 v_{-1}-4 v_0+6 v_1-4 v_2+1 v_3=b \\ i=2 & 1 v_0-4 v_1+6 v_2-4 v_3+1 v_4=b \\ i=3 & 1 v_1-4 v_2+6 v_3-4 v_4+1 v_5=b \\ i=4 & 1 v_2-4 v_3+6 v_4-4 v_5+1 v_6=b \end{array} \]

where \(b=\frac{h^4 \cdot q}{E I}\)

Next, we write the boundary conditions, at x=0, i=0 full support, at x=4, i=4 free end

\[ \begin{aligned} &v_0=0 \\ &-1 v_{-1}+1 v_1=0 \\ &1 v_3-2 v_4+1 v_5=0 \\ &-1 v_2+2 v_3-2 v_5+1 v_6=0 \end{aligned} \]

Similar to the previous example, we write everything in matrix form and solve using a computational program

\[ \left[\begin{array}{cccccccc} 1 & -4 & 6 & -4 & 1 & 0 & 0 & 0 \\ 0 & 1 & -4 & 6 & -4 & 1 & 0 & 0 \\ 0 & 0 & 1 & -4 & 6 & -4 & 1 & 0 \\ 0 & 0 & 0 & 1 & -4 & 6 & -4 & 1 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & -1 & 2 & 0 & -2 & 1 \end{array}\right] \cdot\left[\begin{array}{c} v_{-1} \\ v_0 \\ v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \\ v_6 \end{array}\right]=\left[\begin{array}{c} b \\ b \\ b \\ b \\ 0 \\ 0 \\ 0 \\ 0 \end{array}\right]=\left[\begin{array}{c} 1.25 \cdot 10^{-3} \\ 1.25 \cdot 10^{-3} \\ 1.25 \cdot 10^{-3} \\ 1.25 \cdot 10^{-3} \\ 0 \\ 0 \\ 0 \\ 0 \end{array}\right] \]
\[ \left[\begin{array}{l} v_{-1} \\ v_0 \\ v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \\ v_6 \end{array}\right]=A^{-1} B=\left[\begin{array}{l} 5 \cdot 10^{-3} \\ 0 \\ 5 \cdot 10^{-3} \\ 1.562 \cdot 10^{-2} \\ 2.875 \cdot 10^{-2} \\ 4.25 \cdot 10^{-2} \\ 5.625 \cdot 10^{-2} \\ 7.062 \cdot 10^{-2} \end{array}\right] \]

Having the deflection values, we can calculate, for example, the value of the bending moment at the support

\[ M_0=\frac{-E I}{h^2}\left(v_{-1}-2 \cdot v_0+v_1\right)=-80 kNm \]

This beam is very typical, so we can easily find the analytical solution and compare it with our numerical one

\[ M_a=\frac{-1}{2} q \cdot L^2=-80kNm \]

As can be seen, the solution for the bending moment is identical; the situation is somewhat more interesting with the deflection at the end of the beam:

\[ \begin{aligned} &v_a=\frac{q \cdot L^4}{8 E I}=4 \cdot 10^{-2} \\ &\frac{v_4-v_a}{v_a}=6.25 \% \end{aligned} \]

As can be seen, in this case, the relative error is over 6%; similarly, a higher n would have made this error smaller.

Bonus

In the above example, we omitted the value of the shear force at the support. As can be easily seen, it creates a certain problem because the formula:

\[ T_0=-\frac{E I}{2 h^3}\left(-v_{-2}+2 \cdot v_{-1}-2 v_1+v_2\right) \]

contains the value of displacement \(v_{-2}\) which we did not calculate. However, if the task were to calculate the value of the shear force, we would additionally need to write the difference formula for i=0, co illustrated below

\[ \left[\begin{array}{ccccccccc} 1 & -4 & 6 & -4 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & -4 & 6 & -4 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & -4 & 6 & -4 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & -4 & 6 & -4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 & -4 & 6 & -4 & 1 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 0 & -1 & 2 & 0 & -2 & 1 \end{array}\right] \cdot\left[\begin{array}{l} v_{-2} \\ v_{-1} \\ v_0 \\ v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \\ v_6 \end{array}\right]=\left[\begin{array}{l} b \\ b \\ b \\ b \\ b \\ 0 \\ 0 \\ 0 \\ 0 \end{array}\right] \]
\[ \left[\begin{array}{l} v_{-2} \\ v_{-1} \\ v_0 \\ v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \\ v_6 \end{array}\right]:=A^{-1} B=\left[\begin{array}{l} 2.562 \cdot 10^{-2} \\ 5 \cdot 10^{-3} \\ 0 \\ 5 \cdot 10^{-3} \\ 1.562 \cdot 10^{-2} \\ 2.875 \cdot 10^{-2} \\ 4.25 \cdot 10^{-2} \\ 5.625 \cdot 10^{-2} \\ 7.062 \cdot 10^{-2} \end{array}\right] \]

Of course, the other displacement values have not changed, but having calculated \(v_{-2}\) we can verify the value of the shear force at the support:

\[ \begin{aligned} &T_a=q \cdot L=40 kN \\ &T_0=-\frac{E I}{2 h^3}\left(-v_{-2}+2 \cdot v_{-1}-2 v_1+v_2\right)=40kN \end{aligned} \]

As can be seen, the result is perfectly consistent with the analytical value