Metoda różnic skończonych (MRS)




Edupanda » Numerical Methods » Finite Difference Method (FDM)

Finite Difference Method (FDM)

A method that involves approximating 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 convert individual derivatives into difference formulas using the following formulas:

Central Difference Formulas for the One-Dimensional Problem

Finite Difference Method - Central Difference Formulas for the One-Dimensional Problem

\begin{aligned} &f^I=\frac{1}{2 h}\left(-v_{i-1}+v_{i+1}\right) \\ &f^{I I}=\frac{1}{h^2}\left(v_{i-1}-2 v_i+v_{i+1}\right) \\ &f^{I I I}=\frac{1}{2 h^3}\left(-v_{i-2}+2 v_{i-1}-2 v_{i+1}+v_{i+2}\right) \\ &f^{I V}=\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 with the Finite Difference Method

The general procedure in this case is as follows:

  1. Divide the interval of interest into n elements of length h
  2. Replace derivatives in the equation with appropriate difference formulas
  3. Write down difference equations for each point
  4. Consider boundary conditions
  5. Solve the system of equations in matrix form
solution-shape

Example 1

Content

Solve the boundary problem with the following data

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

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

Solution

Length of a single element:

\begin{aligned} h=\frac{b-a}{n}=1 \end{aligned}

Position of individual nodes:

\begin{aligned} x_i=x_0+i\cdot h \end{aligned}

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

\begin{aligned} &\frac{1}{h^2}\left(x_{i-1}-2 \cdot x_i+x_{i+1}\right)=2 x \\ &x_{i-1}-2 \cdot x_i+x_{i+1}=h^2 \cdot 2 x \end{aligned}

For each node \(i=1,2,3,4\) we write a difference equation according to the above formula, and we also record the value of x for each node:

\begin{array}{lll} i=1 & x_0-2 \cdot x_1+x_2=h^2 \cdot 2 x_1 & x_1=x_0+h=4 \\ i=2 & x_1-2 \cdot x_2+x_3=h^2 \cdot 2 x_2 & x_2=x_0+2 h=5 \\ i=3 & x_2-2 \cdot x_3+x_4=h^2 \cdot 2 x_3 & x_3=x_0+3 h=6 \\ i=4 & x_3-2 \cdot x_4+x_5=h^2 \cdot 2 x_4 & x_4=x_0+4 h=7 \end{array}

Consider the boundary conditions, also converting them into difference formulas:

\begin{array}{lll} y(a)=4 & x_0=4 & \\ y^{\prime}(b)=2 & \frac{1}{2 h}\left(-x_3+x_5\right)=2 & -x_3+x_5=4 h \end{array}

We write down 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, because the last one is actually outside the interval of interest \(i=5 -> x=8\)

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

\( y_{a n}(x)=\frac{x^3}{3}-47 x+136 \)

The table shows the values along with the errors for each point

\( \begin{array}{|r|r|r|r|r|r|} \hline \mathrm{i} & x & y & y_a & \text { Absolute Error } & \text { Relative Error } \\ \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 be higher, but in such cases it would be necessary to use computational programs like Matlab

Solving Bending Beams with the Finite Difference Method

The procedure is almost identical to solving differential equations because we start with the relationship between the displacement and the load on a straight beam:

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

Which, in difference form, takes the following 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 continuously loaded, but for example with a point force applied to it, trapezoidal load, etc., or when its stiffness varies along the length of the beam. In such cases, we will have to convert the load into a continuous load.
Additionally, after calculating the values of beam deflection at various points, we will use known relationships 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 beams depend on the support method as shown below

finite difference method - boundary conditions for beam - fixed support \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}
finite difference method - boundary conditions for beam - hinged support \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}
finite difference method - boundary conditions for beam - sliding support \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}
finite difference method - boundary conditions for beam - free end \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}
solution-shape

Example 2

Content

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

Finite difference method - example 2 - beam

Solution

We write the difference equation for 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}\)

Then we write the boundary conditions, at x=0, i=0 fixed 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}

Similarly 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] \)

With the deflection values, we can calculate, for example, the bending moment at the fixed support

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

Interestingly, this beam is very typical, so we can easily find an analytical solution and compare it with our numerical solution

\( M_a=\frac{-1}{2} q \cdot L^2=-80kNm \)

As we can see, the solution for the bending moment is identical, but the end deflection is more interesting:

\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 we can see, in this case the relative error is over 6%, similarly as before, increasing n would reduce this error.

Bonus

In the above example, we omitted the value of the shear force at the fixed support. As you can easily notice, it creates a 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 displacement value \(v_{-2}\) which we did not calculate. However, if our task was to calculate the value of the shear force, we would need to additionally write a difference equation for i=0, as 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 remaining displacement values have not changed, but with the calculated \(v_{-2}\) we can verify the value of the shear force at the fixed 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 you can see, the result perfectly matches the analytical value

Łukasz Cichowicz
GRZEGORZ
MAZUR
Korepetytor
200 PLN
60 MIN
Kontakt
WhatsApp grzegorz@edupanda.pl
Zobacz mój profil na
Edupanda Logo