Numerical Solution of Bending of the Beam with Given Friction

Abstract

We are interested in a contact problem for a thin fixed beam with an internal point obstacle with possible rotation and shift depending on a given swivel and sliding friction. This problem belongs to the most basic practical problems in, for instance, the contact mechanics in the sustainable building construction design. The analysis and the practical solution plays a crucial role in the process and cannot be ignored. In this paper, we consider the classical Euler–Bernoulli beam model, which we formulate, analyze, and numerically solve. The objective function of the corresponding optimization problem for finding the coefficients in the finite element basis combines a quadratic function and an additional non-differentiable part with absolute values representing the influence of considered friction. We present two basic algorithms for the solution: the regularized primal solution, where the non-differentiable part is approximated, and the dual formulation. We discuss the disadvantages of the methods on the solution of the academic benchmarks.

1. Introduction

Mathematical modeling plays a crucial role in the simulations and the development of sustainable civil engineering. This paper aims to present the whole story of the solution process of a given problem, instead of focusing only on one specific step of the solution process. Additionally, we would like to show that in the case of practical computation, in the end, the numerical solver always matters. Indeed, even though the solver has to be chosen in a convenient way referring to the properties of the given problem, the efficiency of the solution crucially depends on the choice of the right numerical approach; actually, there will always be pros and cons, which have to be taken into the consideration.

In the paper, we examine this rule on algorithms for the numerical solution of the beam bending problem of a thin fixed beam with an internal point obstacle with possible rotation and shift depending on a given swivel and sliding friction. The beam is a structural element that primarily resists loads applied laterally to the axis. In the case of a thin beam, we suppose that the diameter of the cross section is much smaller than the length of the beam. The thin fixed beam considered in this paper is the thin beam fixed on both sides.

We suppose that the length of the beam is much bigger than its height. The considered load function applies to small deflections of a beam without considering the effects of shear deformations. Such an example belongs to the most popular benchmarks in engineering practice because of its linearity, and it is well known as the Euler–Bernoulli beam model (see, e.g., in [1], where authors present the derivation of the model). For the nonlinear beam model, which supposes large deformations, we can choose, for example, the Gao beam model, which has been originally proposed in [2]. The model is commonly used both in practice and theory, for example, the extension for contact problems with obstacle [3], friction [4], the extension incorporating nonlocality and surface energy effect [5], or the application to optimal control [6].

From the mathematical point of view, the classical formulation of the problem is represented as a linear differential equation of the fourth-order with corresponding additional conditions. In this case, the weak formulation is defined as an elliptical variational inequality of the second type with non-differentiable functional. For more details see, e.g., the comprehensive books in [7,8,9]. In this paper, we shortly review the formulation in Section 2.1, where we also shortly discuss the existence and uniqueness of the solution in the case of the continuous formulations and the approximated problem.

Convergence analysis and an algebraic formulation is presented in Section 2.4. To solve the problem numerically, we adopt the Finite Element Method (FEM) with Hermite spline functions (see, e.g., the introductory book in [10] or the comprehensive numerical analysis in [11]), see Section 2.6.

Two methods are used to solve the problem in Section 2.7. The first one is the method of regularization, which is based on the approximation of the non-differentiable functional by the differentiable one. The second approach is the dualization of the problem and its reformulation into the so-called Quadratic programming (QP) problem [12]. The corresponding dual problem in our case is a minimization of a strictly convex quadratic function on a feasible set defined by box constraints.

As an example of the presented theory and methods of solution, we consider three benchmarks with increasing difficulty, see Section 2.8. These benchmarks show the basic properties of adopted numerical approaches. The results are discussed in Section 3.

Finally, Section 4 concludes the paper. Please see the section titled “Appendix A” for the list of notations used in the paper.

2. Materials and Methods

2.1. Problem Formulation

We consider a thin fixed beam of length $l>0$ with special internal point obstacle $\widehat{x}$ with possible shift and rotation depending on a given sliding $g_1\ge 0$ and swivel $g_2\ge 0$ friction, see Figure 1. The geometry of the beam is defined by the moment of inertia of the cross section $J\left(x\right)$, $x\in [0,l]$, and the material of the body is defined by Young’s modulus of elasticity $E\left(x\right)$. The problem is to find the deflection of the beam u, which is caused by the load function f and limited by the considered type of obstacle.

2.2. Classical Formulation of the Problem

To introduce the continuous classical formulation, we assume a beam of length $l\in R^{+}$, load function $f\in C\left(\right(0,l\left)\right)$, $E,J\in C^2\left((0,l)\right)$, and internal point obstacle $\widehat{x}\in (0,l)$, given sliding friction $g_1\ge 0$ and given swivel friction $g_2\ge 0$. We consider the following problem: find a function u, such that

$$\begin{array}{c}\hfill \left(P_{TM}\right)u\in C^4\left((0,l)\right):D^2\left(E\left(x\right)J\left(x\right)D^{2}u\left(x\right)\right)=f\left(x\right)\forall x\in (0,l),\end{array}$$

(1)

$$\begin{array}{c}\hfill u\left(0\right)=Du\left(0\right)=u\left(l\right)=Du\left(l\right)=0,\end{array}$$

(2)

$$\begin{array}{c}\hfill |(T^{-}-T^{+})\left(u\left(\widehat{x}\right)\right)|\le g_1,\end{array}$$

(3)

$$\begin{array}{c}\hfill g_1\left|u\left(\widehat{x}\right)\right|+(T^{-}-T^{+})\left(u\left(\widehat{x}\right)\right)u\left(\widehat{x}\right)=0,\end{array}$$

(4)

$$\begin{array}{c}\hfill |(M^{-}-M^{+})\left(u\left(\widehat{x}\right)\right)|\le g_2,\end{array}$$

(5)

$$\begin{array}{c}\hfill g_2\left|Du\left(\widehat{x}\right)\right|+(M^{-}-M^{+})\left(u\left(\widehat{x}\right)\right)Du\left(\widehat{x}\right)=0.\end{array}$$

(6)

Function u, corresponding to the problem $\left(P_{TM}\right)$, will be called a classical solution of the considered problem.

The differential equation (1) is a mathematical model of the bending of the Euler–Bernoulli beam. The boundary conditions (2) define the clamping of the beam at both ends, relations (3) and (4) represent the conditions of given sliding friction, and relations (5) and (6) are the conditions of the given swivel friction.

2.3. Weak and Variational Formulation

As the classical formulation of the problem is too strong and therefore not suitable for direct numerical solution, we formulate the considered problem as a variational one.

Let $f\in L_2\left((0,l)\right)$, $g_i\ge 0,i=1,2$, and $E,J\in L_{\infty }\left((0,l)\right)$. We define the space V of test functions, corresponding to the boundary conditions, as a Sobolev space (see, e.g., in [11,13]) as follows:

$$V=H_0^2\left((0,l)\right).$$

(7)

Let us denote

$$q_1\left(v\right)=g_1|v\left(\widehat{x}\right)|,q_2\left(v\right)=g_2\left|Dv\left(\widehat{x}\right)\right|.$$

(8)

Next, we introduce the functional $q:V\to \overline{\mathbb{R}}$, corresponding to given swivel and sliding frictions, in the form

$$q\left(v\right)=q_1\left(v\right)+q_2\left(v\right).$$

(9)

Using (8), we have

$$q\left(v\right)=g_1|v\left(\widehat{x}\right)|+g_2\left|Dv\left(\widehat{x}\right)\right|.$$

(10)

Let u be a solution of the problem $\left(P_{TM}\right)$ and $v\in V$. By scalar multiplication of Equation (1) by the function $v-u$ and integration (using Green’s theorem—see, e.g., [14]) on the interval $(0,l)$, we obtain

$$\begin{array}{c}\underset{0}{\overset{l}{\int }}{D}^2\left(E\left(x\right)J\left(x\right)D^{2}u\left(x\right)\right)(v-u)\left(x\right)dx=\underset{0}{\overset{l}{\int }}E\left(x\right)J\left(x\right)D^{2}u\left(x\right)D^2(v-u)\left(x\right)dx+\\ \\ -(T^{-}-T^{+})\left(u\left(\widehat{x}\right)\right)(v-u)\left(\widehat{x}\right)+(M^{-}-M^{+})\left(u\left(\widehat{x}\right)\right)D(v-u)\left(\widehat{x}\right)=\end{array}$$

$$=\underset{0}{\overset{l}{\int }}f(v-u)\left(x\right)dx\forall v\in V.$$

(11)

From (4) and (6), we have

$$(T^{-}-T^{+})\left(u\left(\widehat{x}\right)\right)u\left(\widehat{x}\right)=-\left|u\left(\widehat{x}\right)\right|g_1,$$

(12)

$$(M^{-}-M^{+})\left(u\left(\widehat{x}\right)\right)Du\left(\widehat{x}\right)=\left|Du\left(\widehat{x}\right)\right|g_2.$$

(13)

Using (12), (13) and (3), (5), it follows from (11) that

$$\begin{array}{c}0=\underset{0}{\overset{l}{\int }}{D}^2\left(E\left(x\right)J\left(x\right)D^{2}u\left(x\right)\right)(v-u)\left(x\right)dx-\underset{0}{\overset{l}{\int }}f(v-u)\left(x\right)dx=\\ \\ =\underset{0}{\overset{l}{\int }}E\left(x\right)J\left(x\right)D^{2}u\left(x\right)D^2(v-u)\left(x\right)dx-(T^{-}-T^{+})\left(u\left(\widehat{x}\right)\right)v\left(\widehat{x}\right)+\\ \\ +(T^{-}-T^{+})\left(u\left(\widehat{x}\right)\right)u\left(\widehat{x}\right)+(M^{-}-M^{+})\left(u\left(\widehat{x}\right)\right)Dv\left(\widehat{x}\right)+\\ \\ -(M^{-}-M^{+})\left(u\left(\widehat{x}\right)\right)Du\left(\widehat{x}\right)-\underset{0}{\overset{l}{\int }}f(v-u)\left(x\right)dx\le \\ \\ \le \underset{0}{\overset{l}{\int }}E\left(x\right)J\left(x\right)D^{2}u\left(x\right)D^2(v-u)\left(x\right)dx+|(T^{-}-T^{+})\left(u\left(\widehat{x}\right)\right)\left|\right|v\left(\widehat{x}\right)|+\\ \\ -|u\left(\widehat{x}\right)|g_1+|(M^{-}-M^{+})\left(u\left(\widehat{x}\right)\right)\left|\right|Dv\left(\widehat{x}\right)|+\\ \\ -|Du\left(\widehat{x}\right)|g_2-\underset{0}{\overset{l}{\int }}f(v-u)\left(x\right)dx\le \\ \\ \le \underset{0}{\overset{l}{\int }}E\left(x\right)J\left(x\right)D^{2}u\left(x\right)D^2(v-u)\left(x\right)dx+|v\left(\widehat{x}\right)|g_1-\left|u\left(\widehat{x}\right)\right|g_1+\\ \\ +|Dv\left(\widehat{x}\right)|g_2-|Du\left(\widehat{x}\right)|g_2-{\int }_0^{l}f(v-u)\left(x\right)dx\forall v\in V.\end{array}$$

Afterwards, we can define a weak formulation of the problem $\left(P_{TM}\right)$ as a problem of finding a function u such that

$$\begin{array}{c}\hfill \left({\mathcal{P}}_{TM}\right)u\in V:a(u,v-u)+q\left(v\right)-q\left(u\right)\ge F(v-u)\forall v\in V,\end{array}$$

(14)

where

$$a(u,v)=\underset{0}{\overset{l}{\int }}E\left(x\right)J\left(x\right)D^{2}u\left(x\right)D^{2}v\left(x\right)dx,$$

(15)

$$F\left(v\right)=(f,v)={\int }_0^{l}fv\left(x\right)dx.$$

(16)

The sought function u is a weak solution of the problem $\left(P_{TM}\right)$.

Remark 1.

The weak formulation $\left({\mathcal{P}}_{TM}\right)$ of the problem $\left(P_{TM}\right)$ is formulated as an elliptical variational inequality of the second type with non-differentiable functional q (see, e.g., in [8]).

Let form $a:V\times V\to \mathbb{R}$ and functional $F:V\to \mathbb{R}$ be given by relations (15) and (16) and $f\in L_2\left((0,l)\right)$, $E,J\in L_{\infty }\left((0,l)\right)$ such that

$$E\left(x\right)\ge E_0>0almosteverywhereon(0,l),$$

(17)

$$J\left(x\right)\ge J_0>0almosteverywhereon(0,l).$$

(18)

Obviously, the form a is bilinear and symmetric. To prove its continuity we use Hölder’s inequality (see, e.g., in [7])

$$\left|a(u,v)\right|=|{\int }_0^{l}E\left(x\right)J\left(x\right)D^{2}u\left(x\right)D^{2}v\left(x\right)dx|\le c{\int }_0^l|D^{2}u\left(x\right)\left|\right|D^{2}v\left(x\right)|\le$$

$$\le c({\int }_0^l|D^2{u\left(x\right)|}^2{)}^{\frac{1}{2}}({\int }_0^l|D^2{v\left(x\right)|}^2{)}^{\frac{1}{2}}\le \tilde{c}{\parallel u\parallel }_{H^2}{\parallel v\parallel }_{H^2}\forall v,u\in V,$$

where $c={\parallel E\parallel }_{L_{\infty }}{\parallel J\parallel }_{L_{\infty }}$. From the assumptions (17) and (18) of functions E, J and from $v^2\left(0\right)=0$, we have

$$a(v,v)={\int }_0^{l}E\left(x\right)J\left(x\right){\left(D^{2}v\left(x\right)\right)}^{2}dx\ge E_0{J}_0({\int }_0^l{\left(D^{2}v\left(x\right)\right)}^{2}dx+v^2\left(0\right)).$$

Here, we apply Friedrichs’s inequality (see, e.g., [7]) to obtain the V-ellipticity of the form a

$$a(v,v)\ge \frac{E_0{J}_0}{k}{\parallel v\parallel }_{H^2}^2=c{\parallel v\parallel }_{H^2}^2\forall v\in V,$$

where $c=\frac{E_0{J}_0}{k}>0$ and k is the positive constant from Friedrichs’s inequality. It is also clear that the functional F is linear. Using Cauchy–Schwarz inequality (see, e.g., in [7]), we can prove its continuity, i.e., $F\in V^{*}$.

Let $t\in (0,1)$. Then, for functional q given by (10), it holds that

$$q\left((1-t)v+tu\right)=g_1|(1-t)v\left(\widehat{x}\right)+tu\left(\widehat{x}\right)|+g_2\left|D[(1-t)v\left(\widehat{x}\right)+tu\left(\widehat{x}\right)]\right|\le$$

$$(1-t)g_1|v\left(\widehat{x}\right)|+t{g}_1|u\left(\widehat{x}\right)|+(1-t)g_2|Dv\left(\widehat{x}\right)|+t{g}_2\left|Du\left(\widehat{x}\right)\right|=$$

$$=(1-t)[g_1|v\left(\widehat{x}\right)|+g_2|Dv\left(\widehat{x}\right)\left|\right]+t[g_1|u\left(\widehat{x}\right)|+g_2|Du\left(\widehat{x}\right)\left|\right]=$$

$$=(1-t)j\left(v\right)+tj\left(u\right),$$

which means that q is convex. Immediately from the definition of q, we can see that $q\left(v\right)\ge 0>-\infty$ $\forall v\in V$. Moreover, for zero function $v\in V$, it holds $q\left(v\right)=0\ne \infty$, i.e., $q\left(V\right)\not\equiv \infty$ and thus q is proper on V. It is clear that the functional q is continuous. The convexity and continuity of q imply its semi-continuity from below. Therefore, the existence and uniqueness of the solution of problem $\left({\mathcal{P}}_{TM}\right)$ are guaranteed by the following Theorem 1, the assumptions of which are thus fulfilled.

Theorem 1.

Let $a:V\times V\to \mathbb{R}$ be bounded, V-elliptical bilinear form, $F\in V^{*}$ and let functional $j:V\to \overline{\mathbb{R}}$ be convex, semi-continuous from below and proper on V. Then, there exists an unique solution of

$$\begin{array}{c}\hfill a(u,v-u)+j\left(v\right)-j\left(u\right)\ge F(v-u)\forall v\in V.\end{array}$$

(19)

For more detail, see, e.g., in [8,15].

As a is bounded; the V-elliptical is bilinear and of symmetric form; functional q is convex, semi-continuous from below, and proper on V; and the problem $\left({\mathcal{P}}_{TM}\right)$ is equivalent to the variational problem of minimizing quadratic energy functional, i.e.,

$$\left({\widehat{\mathcal{P}}}_{TM}\right)u\in V:u=arg\underset{v\in V}{min}\mathcal{J}\left(v\right),$$

(20)

where

$$\mathcal{J}\left(v\right)=\frac{1}{2}a(v,v)+q\left(v\right)-F\left(v\right).$$

(21)

2.4. Approximation

Let $f\in L_2\left((0,l)\right)$ and $E,J\in L_{\infty }\left((0,l)\right)$ satisfy relations (17) and (18), respectively. Let the form $a:V\times V\to \mathbb{R}$ and the functional $F:V\to \mathbb{R}$ be given by (15) and (16), respectively. We consider a system of partitions $\left\{D_h\right\}$, $h\to 0^{+}$ of the interval $[0,l]$ into subintervals

$$T_{k\left(h\right)}=[x_{k\left(h\right)-1},x_{k\left(h\right)}],k\left(h\right)=1,\cdots ,n\left(h\right)$$

such that points $x_i,i=1,\cdots ,n\left(h\right)-1$ are nodal points, i.e.,

$$0=x_0<x_{1\left(h\right)}<\cdots <x_{n\left(h\right)}=l\forall D_h,$$

where $n\left(h\right)+1$ is the number of nodal points of $D_h$ and h is the maximum length of intervals $T_{k\left(h\right)}$. Moreover, the point obstacle $\widehat{x}\in D_h$ for any h. We define the following finite-dimensional subspace $V_h\subset V$:

$$\begin{array}{c}{V}_h={\{v_h\in C^1\left([0,l]\right):v_h|}_{T_{k\left(h\right)}}\in {\mathcal{P}}_3\forall T_{k\left(h\right)}\in D_h,\\ \\ v_h\left(0\right)=D{v}_h\left(0\right)=v\left(l\right)=Dv\left(l\right)=0\}\forall D_h.\end{array}$$

Therefore, spaces $V_h$ satisfy classical boundary conditions and additionally

$$\dim V_h=m\left(h\right)<\infty ,\mathrm{where}m\left(h\right)=2(n\left(h\right)-1)$$

(22)

and

$$\underset{h\to 0^{+}}{lim}\dim V_h=\underset{h\to 0^{+}}{lim}m\left(h\right)=\infty .$$

(23)

We approximate the functional q for each h by functional $q_h:V_h\to \overline{\mathbb{R}}$ by

$$q_h\left(v_h\right)=g_1|v_h\left(\widehat{x}\right)|+g_2\left|D{v}_h\left(\widehat{x}\right)\right|.$$

(24)

Functional $q_h$ is convex, semi-continuous from below and proper in the space $V_h$. We can approximate problem $\left({\mathcal{P}}_{TM}\right)$ by the sequence of problems of finding $u_h$ such that

$$\begin{array}{c}\hfill {\left({\mathcal{P}}_{TM}\right)}_h{u}_h\in V_h:a(u_h,v_h-u_h)+q_h\left(v_h\right)-q_h\left(u_h\right)\ge F(v_h-u_h)\forall v_h\in V_h.\end{array}$$

(25)

The existence and uniqueness of the solution of ${\left({\mathcal{P}}_{TM}\right)}_h$ is guaranteed by the Theorem 1.

As the form a is symmetric, the problem ${\left({\mathcal{P}}_{TM}\right)}_h$ is equivalent to the problem of finding $u_h$ such that

$${\left({\widehat{\mathcal{P}}}_{TM}\right)}_h{u}_h\in V_h:u_h=arg\underset{v_h\in V_h}{min}{\mathcal{J}}_h\left(v_h\right),$$

(26)

where

$${\mathcal{J}}_h\left(v_h\right)=\frac{1}{2}a(v_h,v_h)-F\left(v_h\right)+q_h\left(v_h\right).$$

(27)

To approximate the function $v\in H^2\left((0,l)\right)$, we use the Hermitian cubic spline function for each $D_h$ such that

$$v_h\left(x_i\right)=v\left(x_i\right),D{v}_h\left(x_i\right)=Dv\left(x_i\right)\forall x_i\in D_h,\forall v\in H^2\left((0,l)\right).$$

(28)

Additionally,

$$\underset{h\to 0^{+}}{lim}{v}_h=v\forall v\in H^2\left((0,l)\right).$$

(29)

Remark 2.

For the given function $v\in H^2\left((0,l)\right)$, the Hermitian cubic spline function is given by

$$\begin{array}{c}{v}_h\left(x\right)=v\left(x_k\right)[2(x_k-x)-(x_{k+1}-x_k)]\frac{{(x-x_{k+1})}^2}{{(x_k-x_{k+1})}^3}+Dv\left(x_k\right)(x-x_k)\frac{{(x-x_{k+1})}^2}{{(x_{k+1}-x_k)}^2}+\\ \\ +v\left(x_{k+1}\right)[2(x_{k+1}-x)-(x_{k+1}-x_k)]\frac{{(x-x_k)}^2}{{(x_k-x_{k+1})}^3}+Dv\left(x_{k+1}\right)(x-x_{k+1})\frac{{(x-x_k)}^2}{{(x_{k+1}-x_k)}^2},\\ \\ x\in [x_k,x_{k+1}],k=0,...,n\left(h\right)-1.\end{array}$$

2.5. Convergence Analysis

To prove the convergence of the approximated solutions $u_h$ to the weak solution u of the considered problem, we use the following theorem.

Theorem 2.

Let $a:V\times V\to R$ be bounded, V-elliptical bilinear form on V. Let there exist operator $r_h:X\to V_h$ such that

$$r_{h}v\to vforh\to 0^{+}\forall v\in X,$$

(30)

where $X\subset V$ such that $\overline{X}=V$. Let the system of functionals $j_h$ meet the following two conditions:

$$v_h\rightharpoonup v,v_h\in V_h\Rightarrow \underset{h\to 0^{+}}{lim\; inf}{j}_h\left(v_h\right)\ge j\left(v\right),$$

(31)

$$\underset{h\to 0^{+}}{lim}{j}_h\left(r_{h}v\right)=j\left(v\right)\forall v\in X.$$

(32)

Then, it holds that

$$\parallel u-u_h{\parallel }_V\to 0 for h\to 0^{+}.$$

Moreover,

$$\underset{h\to 0^{+}}{lim}{j}_h\left(u_h\right)=j\left(u\right).$$

For more details, see, e.g., in [16].

As the functional $q_h$ is convex and semi-continuous from below on $V_h$, it is also weakly semi-continuous from below. In our case, operators $r_h$ are operators of the appropriate Hermitian interpolation. From the (28), it follows that

$$q_h\left(v_h\right)=q\left(v\right)\forall v\in V.$$

Therefore,

$$\underset{h\to 0^{+}}{lim}{q}_h\left(r_{h}v\right)=q\left(v\right)\forall v\in V.$$

This means that assumptions of Theorem 2 are satisfied, i.e., the sequence $u_h$ of solution approximations of the problem ${\left({\mathcal{P}}_{TM}\right)}_h$ converges to the solution u of the problem $\left({\mathcal{P}}_{TM}\right)$ in the V-space norm.

2.6. The Algebraic Formulation

Let the form a and the functional F be defined by (15) and (16), respectively. We choose a finite-dimensional space $V_h$, $dim{V}_h=N$ with $N=2\left(n\right(h)-1)$. Let $\left\{{\phi }_i\right\},i=1,\dots ,N$ be a basis of $V_h$ chosen by FEM, i.e., $\forall i,j=1,...,n\left(h\right)-1$

$$\begin{array}{cc}{\phi }_{2i-1}\left(x_j\right)={\delta }_{ij};\hfill & {\phi }_{2i}\left(x_j\right)=0;\hfill \\ D{\phi }_{2i-1}\left(x_j\right)=0;\hfill & D{\phi }_{2i}\left(x_j\right)={\delta }_{ij}.\hfill \end{array}$$

To obtain a shape of basis functions, we apply the Hermite’s cubic splines. We have

$${\phi }_{2i}\left(x\right)=\left\{\begin{array}{cc}0\hfill & x\notin [x_{i-1},x_{i+1}],\hfill \\ \frac{{(x-x_i)}^3}{{(x_{i-1}-x_i)}^2}-2\frac{{(x-x_i)}^2}{x_{i-1}-x_i}+(x-x_i)\hfill & x\in [x_{i-1},x_i],\hfill \\ \frac{{(x-x_i)}^3}{{(x_{i+1}-x_i)}^2}-2\frac{{(x-x_i)}^2}{x_{i+1}-x_i}+(x-x_i)\hfill & x\in [x_i,x_{i+1}],\hfill \end{array}\right.$$

$${\phi }_{2i-1}\left(x\right)=\left\{\begin{array}{cc}0\hfill & x\notin [x_{i-1},x_{i+1}],\hfill \\ 1+2{\left(\frac{x-x_i}{x_{i-1}-x_i}\right)}^3-3{\left(\frac{x-x_i}{x_{i-1}-x_i}\right)}^2\hfill & x\in [x_{i-1},x_i],\hfill \\ 1+2{\left(\frac{x-x_i}{x_{i+1}-x_i}\right)}^3-3{\left(\frac{x-x_i}{x_{i+1}-x_i}\right)}^2\hfill & x\in [x_i,x_{i+1}].\hfill \end{array}\right.$$

For an example of geometric representation, please see Figure 2.

Any function $v_h\in V_h$ can be uniquely written in the form of linear combination of basis functions, i.e.,

$$v_h=\sum _{i=1}^N{c}_i{\phi }_i,$$

(33)

where $c_i\in \mathbb{R}$ are the corresponding coordinates of $v_h$ in the basis. It is clear that ${\phi }_i$, $c_i$ depend on h. If we substitute $v_h$ from (33) into the functional ${\mathcal{J}}_h$ given by (27), we obtain

$${\mathcal{J}}_h\left(v_h\right)={\mathcal{J}}_h(\sum _k{c}_k{\phi }_k)=\frac{1}{2}a(\sum _k{c}_k{\phi }_k,\sum _j{c}_j{\phi }_j)-F(\sum _k{c}_k{\phi }_k)+$$

$$+q_h(\sum _k{c}_k{\phi }_k)=\frac{1}{2}\sum _k\sum _j{c}_k{c}_{j}a({\phi }_k,{\phi }_j)-\sum _k{c}_{k}F\left({\phi }_k\right)+$$

$$+q_h(\sum _k{c}_k{\phi }_k)=\frac{1}{2}(\mathit{c},\mathit{A}\mathit{c})+q_h(\mathsf{\Phi },\mathit{c})-(\mathit{b},\mathit{c}),$$

where

• $\mathit{A}={\left(a_{kj}\right)}_{k,j=1}^N$ is a stiffness matrix with elements $a_{kj}=a({\phi }_k,{\phi }_j)$,
• $\mathit{b}=(b_1,\dots ,b_N)$ is a vector of right-hand sides with elements $b_i=F\left({\phi }_i\right)$,
• $\mathit{c}=(c_1,\dots ,c_N)$ is a vector of unknown coefficients of linear combination (33),
• $\mathsf{\Phi }=({\phi }_1,...,{\phi }_N)$ is vector of base functions ${\phi }_i,i=1,...,N$.

Functional ${\mathcal{J}}_h$ is expressed as a quadratic function of variables $c_1,\dots ,c_N$. Let us define an isomorphism $\mathcal{T}:V_h\to {\mathbb{R}}^N$ by relation

$$\mathcal{T}{v}_h=\mathit{c}=(c_1,\dots ,c_N)\in {\mathbb{R}}^N\forall v_h\in V_h.$$

(34)

Then, the problem ${\left({\widehat{\mathcal{P}}}_{TM}\right)}_h$ can be rewritten in the equivalent form where it is required to find ${\mathit{c}}^{*}$ such that

$${\left({\widehat{\mathcal{P}}}_{TM}\right)}_N{\mathit{c}}^{*}\in {\mathbb{R}}^N:{\mathit{c}}^{*}=arg\underset{\mathit{c}\in {\mathbb{R}}^N}{min}{\tilde{\mathcal{J}}}_h\left(\mathit{c}\right),$$

(35)

where

$${\tilde{\mathcal{J}}}_h\left(\mathit{c}\right)={\mathcal{J}}_h\left({\mathcal{T}}^{-1}\mathit{c}\right)=\frac{1}{2}(\mathit{c},\mathit{A}c)+q_h(\mathsf{\Phi },\mathit{c})-(\mathit{b},\mathit{c}),$$

(36)

$$q_h\left(\mathsf{\Phi },\mathit{c}\right)=g_1|c_{2\widehat{j}-1}|+g_2\left|c_{2\widehat{j}}\right|,$$

(37)

${\mathcal{T}}^{-1}$ is inverse of $\mathcal{T}$, and $\widehat{j}$ is index of the nodal point $x_{\widehat{j}}=\widehat{x}$.

Problem ${\left({\widehat{\mathcal{P}}}_{TM}\right)}_N$ is the algebraic form of problem ${\left({\widehat{\mathcal{P}}}_{TM}\right)}_h$.

2.7. Numerical Realization

In this section, we examine two different approaches for the practical numerical solution of the problem. For simplicity, we consider an equidistant grid

$$0=x_0<x_1<\cdots <x_n=l,x_j=jh,j=0,\dots ,n,$$

where $n+1$ is a number of nodal points and $h=l/n$ is a size of intervals. Following the Ritz–Galerkin method, we approximate V by finite-dimensional space

$$\begin{array}{c}{V}_h={\{v_h\in C^1\left([0,l]\right):v_h|}_{[x_j,x_{j+1}]}\in {\mathcal{P}}_3\forall j=0,\dots ,n-1\\ \mathrm{and}{v}_h\left(0\right)=D{v}_h\left(0\right)=v_h\left(l\right)=D{v}_h\left(l\right)=0\}.\end{array}$$

Additionally, we assume that there exists $\widehat{j}$ such that $\widehat{x}=x_{\widehat{j}}$, i.e., the point with defined conditions (3)–(6) is one of the used nodal points.

Remark 3.

Note that $\mathit{dim}{V}_h=N=2(n-1)$ because the solution in boundary points is given by conditions (2).

The optimization problem ${\left({\widehat{\mathcal{P}}}_{TM}\right)}_N$ can be written in the form

$${\left({\tilde{\mathcal{P}}}_{TM}\right)}_N{\mathit{c}}^{*}=arg\underset{c\in {\mathbb{R}}^N}{min}{\tilde{\mathcal{J}}}_N\left(\mathit{c}\right),\widetilde{{\mathcal{J}}_N}\left(\mathit{c}\right)=\psi \left(\mathit{c}\right)+q_N\left(\mathit{c}\right)$$

(38)

with

$$\psi \left(\mathit{c}\right)=\frac{1}{2}{\mathit{c}}^T\mathit{A}\mathit{c}-{\mathit{b}}^T\mathit{c},q_N\left(\mathit{c}\right)=g_1|c_{2\widehat{j}-1}|+g_2\left|c_{2\widehat{j}}\right|,$$

(39)

symmetric positive definite (SPD) matrix $\mathit{A}\in {\mathbb{R}}^{N,N}$, and vector $\mathit{b}\in {\mathbb{R}}^N$.

The solution ${\mathit{c}}^{*}\in {\mathbb{R}}^N$ of (38) represents the coordinates (33) of unknown discretized deflection $u_h$ in the considered basis.

Although the problem (38) has a simple structure, we are still dealing with a nonlinear optimization problem with a non-differentiable objective function. In our case, we are dealing with absolute values in $q_N\left(\mathit{c}\right)$ (39). The optimization problems with absolute values in objective functions arise in various applications, and they are already well examined. For example, when one uses L1-regularization (i.e., so-called least absolute shrinkage and selection operator-LASSO [17]) for regularization of linear regression, the problem has a similar form to (38).

Generally, there exist two types of numerical methods to solve the problem-approximation of a non-differentiable term in objective function or dualization of the problem.

2.7.1. Method of Regularization

The idea is to replace (approximate) the non-differentiable term in objective function with a more suitable differentiable one. In our case, we decided to use the piecewise quadratic approximation, i.e., we approximate absolute value $\omega \left(z\right)=\left|z\right|$ by

$${\tilde{\omega }}_{\epsilon }\left(z\right)=\left\{\begin{array}{cc}-z-\frac{\epsilon }{2}\hfill & \mathrm{if}z<-\epsilon ,\hfill \\ \frac{z^2}{2\epsilon }\hfill & \mathrm{if}z\in [-\epsilon ,\epsilon ],\hfill \\ z-\frac{\epsilon }{2}\hfill & \mathrm{if}z>\epsilon ,\hfill \end{array}\right.$$

(40)

where $\epsilon >0$ is a sufficiently small regularization parameter. The example of this approximation is presented in Figure 3. The derivative is given by

$${\tilde{\omega }}_{N,\epsilon }^{\prime }\left(z\right)=\left\{\begin{array}{cc}-1\hfill & \mathrm{if}z<-\epsilon ,\hfill \\ \frac{z}{\epsilon }\hfill & \mathrm{if}z\in [-\epsilon ,\epsilon ],\hfill \\ 1\hfill & \mathrm{if}z>\epsilon .\hfill \end{array}\right.$$

(41)

Using this approximation, the non-differentiable function $q_N\left(\mathsf{\Phi },\mathbf{c}\right)$ (39) can be written in form

$${\tilde{q}}_{N,\epsilon }\left(\mathit{c}\right)=g_1{\tilde{\omega }}_{\epsilon }\left(c_{2\widehat{j}-1}\right)+g_2{\tilde{\omega }}_{\epsilon }\left(c_{2\widehat{j}}\right).$$

(42)

AS we use piecewise quadratic approximation (40), the derivative is piecewise linear function (41). Additionally, the gradient of quadratic function $\psi$ (39) is linear, and consequently the necessary optimality condition for unconstrained optimization problem (38) is given by the system of linear equations

$$\mathit{A}\mathit{c}-\mathit{b}+\nabla {\tilde{q}}_{N,\epsilon }\left(\mathit{c}\right)=0.$$

The components of gradient ${\tilde{q}}_{N,\epsilon }\left(\mathit{c}\right)\in {\mathbb{R}}^N$ are equal to zero except for the components corresponding to partial derivatives of $c_{2\widehat{j}-1}$ and $c_{2\widehat{j}}$. These values depend on conditions from definition of (40) with respect to values $c_{2\widehat{j}-1}$ and $c_{2\widehat{j}}$. The values are presented in

Table 1. The non-zero components of gradient of the approximated non-differentiable term in the objective function.

		${\mathit{c}}_{\mathbf{2}\widehat{\mathit{j}}-\mathbf{1}}$
		$(\infty ,-\mathit{\epsilon })$	$[-\mathit{\epsilon },\mathit{\epsilon }]$	$(\mathit{\epsilon },\infty )$
$c_{2\widehat{j}}$	$(-\infty ,-\epsilon )$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=-g_1$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=\frac{g_1}{\epsilon }{c}_{2\widehat{j}-1}$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=g_1$
		${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=-g_2$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=-g_2$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=-g_2$
	$[-\epsilon ,\epsilon ]$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=-g_1$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=\frac{g_1}{\epsilon }{c}_{2\widehat{j}-1}$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=g_1$
		${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=\frac{g_2}{\epsilon }{c}_{2\widehat{j}}$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=\frac{g_2}{\epsilon }{c}_{2\widehat{j}}$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=\frac{g_2}{\epsilon }{c}_{2\widehat{j}}$
	$(\epsilon ,\infty )$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=-g_1$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=\frac{g_1}{\epsilon }{c}_{2\widehat{j}-1}$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}-1}=g_1$
		${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=g_2$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=g_2$	${[\nabla {\tilde{q}}_{N,\epsilon }\left(\mathbf{c}\right)]}_{2\widehat{j}}=g_2$

.

The final system of linear equations can be written in form

$$\widehat{\mathit{A}}\mathit{c}=\widehat{\mathit{b}},$$

(43)

where $\widehat{\mathit{A}}$ is a matrix $\mathit{A}$ with diagonal elements modified by coefficients of linear terms in $\nabla {\tilde{q}}_{N,\epsilon }$ and $\widehat{\mathit{b}}$ is a vector b with additional constant terms from $\nabla {\tilde{q}}_{N,\epsilon }$.

It is necessary to solve the problem for all possible cases with respect to $c_{2\widehat{j}-1},c_{2\widehat{j}}$. The solution of the problem ${\left({\widehat{\mathcal{P}}}_{TM}\right)}_N$ (38) is the vector $\mathit{c}$, for which the corresponding condition on $c_{2\widehat{j}-1},c_{2\widehat{j}}$ with respect to $\epsilon$ is satisfied.

2.7.2. Dual Problem

The idea of this approach is based on simple observation of the following equivalent form of absolute value:

$$\forall \alpha \in {\mathbb{R}}^{+}\forall x\in \mathbb{R}:\alpha \left|x\right|=\underset{-\alpha \le \lambda \le \alpha }{max}\lambda x.$$

(44)

Using this identity, we can rewrite the function $q_N(\mathsf{\Phi },\mathit{c})$ (39) into

$$q_N(\mathsf{\Phi },\mathit{c})=\psi \left(\mathit{c}\right)+\left(\underset{-g_1\le {\lambda }_1\le g_1}{max}{\lambda }_1{c}_{2\widehat{j}-1}\right)+\left(\underset{-g_2\le {\lambda }_2\le g_2}{max}{\lambda }_2{c}_{2\widehat{j}}\right)=\underset{-\mathit{g}\le \mathit{\lambda }\le g}{max}{\mathit{\lambda }}^T\mathit{B}\mathit{c},$$

where we denoted

$$\mathit{\lambda }=\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\end{array}\right],\mathit{g}=\left[\begin{array}{c}{g}_1\\ g_2\end{array}\right],$$

(45)

$$\mathit{B}\in {\mathbb{R}}^{2,N}:B_{1,i}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}i=2\widehat{j}-1,\hfill \\ 0\hfill & \mathrm{elsewhere},\hfill \end{array}\right.\mathrm{and}{B}_{2,i}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}i=2\widehat{j},\hfill \\ 0\hfill & \mathrm{elsewhere}.\hfill \end{array}\right.$$

(46)

Using the separability of variables and the saddle-point property [12], we can write optimization problem (38) in equivalent form

$$min\left(\psi \left(\mathit{c}\right)+\underset{-\mathit{g}\le \mathit{\lambda }\le \mathit{g}}{max}{\mathit{\lambda }}^T\mathit{B}\mathit{c}\right)=min\underset{-\mathit{g}\le \lambda \le \mathit{g}}{max}\left(\psi \left(\mathit{c}\right)+{\mathit{\lambda }}^T\mathit{B}\mathit{c}\right)=\underset{-\mathit{g}\le \mathit{\lambda }\le \mathit{g}}{max}min\underset{=L(\mathit{c},\mathit{\lambda })}{\underbrace{\left(\psi \left(\mathit{c}\right)+{\mathit{\lambda }}^T\mathit{B}\mathit{c}\right)}}.$$

(47)

Function $L:{\mathbb{R}}^N\times 2\to \mathbb{R}$ can be considered as a Lagrange function. From the first Karush–Kuhn–Tucker optimality condition [18], we can derive (note that $\mathit{A}$ is SPD, i.e., non-singular)

$${\nabla }_{c}L(\mathit{c},\mathit{\lambda })=\nabla \psi \left(\mathit{c}\right)+{\mathit{B}}^T\mathit{\lambda }=\mathit{A}\mathit{c}-\mathit{b}+{\mathit{B}}^T\mathit{\lambda }=0\Rightarrow \mathit{c}={\mathit{A}}^{-1}(\mathit{b}-{\mathit{B}}^T\mathit{\lambda })$$

(48)

and substitute into objective function L to get (for more details see in [19])

$$\frac{1}{2}{\mathit{c}}^T\mathit{A}\mathit{c}-{\mathit{b}}^T\mathit{c}+{\mathit{\lambda }}^T\mathit{B}\mathit{c}=-\frac{1}{2}{\mathit{\lambda }}^T\mathit{B}{\mathit{A}}^{-1}{\mathit{B}}^T\mathit{\lambda }+{\mathit{\lambda }}^T\mathit{B}{\mathit{A}}^{-1}\mathit{b}-\frac{1}{2}{\mathit{b}}^T{\mathit{A}}^{-1}\mathit{b}.$$

As the original problem (47) is a maximization problem, we can change the sign of the objective function to obtain

$${\mathit{\lambda }}^{*}=arg\underset{-\mathit{g}\le \mathit{\lambda }\le \mathit{g}}{min}\frac{1}{2}{\mathit{\lambda }}^T\widehat{\mathit{A}}\mathit{\lambda }-{\widehat{\mathit{b}}}^T\mathit{\lambda },\widehat{\mathit{A}}=\mathit{B}{\mathit{A}}^{-1}{\mathit{B}}^T,\widehat{\mathit{b}}=\mathit{B}{\mathit{A}}^{-1}\mathit{b}.$$

(49)

Please, notice that the dual problem (49) is a minimization of a strictly convex quadratic function on a feasible set defined by box constraints. Such a problem always has a unique solution [12]. The dimension of the problem is 2, which is much more lower number then the dimension of original primal problem (38). Moreover, the problem belongs to the most basic nonlinear optimization problems and it is solvable by several types of methods, for example, Interior Point, Active set, or Projected Gradient methods [18].

2.8. Numerical Experiments

As a demonstration of the presented theory and methods of solution, we consider three numerical benchmarks: with analytical solution (Benchmark 1), with non-trivial load function (Benchmark 2), and with more internal points with given sliding and swivel friction (Benchmark 3).

Benchmark 1

We consider a most simple case: we suppose problem $\left(P_{TM}\right)$ (1)–(6), i.e., optimization problem ${\left({\tilde{\mathcal{P}}}_{TM}\right)}_N$ (38), with constant $E,J,f\in \mathbb{R}$. It can be easily shown that this problem has an analytical solution. In this paper, the solution is used for measuring the absolute error of proposed numerical algorithms. The form of solution is given by

$$u\left(x\right)=\left\{\begin{array}{c}\gamma x^2{(x-\widehat{x})}^2+\frac{{\widehat{u}}_1}{{\widehat{x}}^2}{x}^2(x-\widehat{x})-\frac{{\widehat{u}}_0}{{\widehat{x}}^3}{x}^2(2x-3\widehat{x})\mathrm{for}x\in (0,\widehat{x})\hfill \\ \gamma {(x-l)}^2{(x-\widehat{x})}^2+\frac{{\widehat{u}}_1}{{(\widehat{x}-l)}^2}{(x-l)}^2(x-\widehat{x})+\frac{{\widehat{u}}_0}{{(\widehat{x}-l)}^3}{(x-l)}^2(-2x+3\widehat{x}-l)\hfill \\ \mathrm{for}x\in (\widehat{x},l)\hfill \end{array}\right.$$

(50)

with $\gamma =\frac{f}{24EJ}$ and unknown ${\widehat{u}}_0=u\left(\widehat{x}\right)$ and ${\widehat{u}}_1=Du\left(\widehat{x}\right)$. These unknown constants can be computed from the conditions (3), (5), (4), and (6). We derived and simplified the corresponding derivatives and form the system of nonlinear equation and inequalities:

$$\begin{array}{ccc}\hfill |{\alpha }_3{\widehat{u}}_0+{\alpha }_2{\widehat{u}}_1-{\beta }_1|& \le & g_1,\hfill \\ \hfill g_1\left|{\widehat{u}}_0\right|-{\beta }_1{\widehat{u}}_0+{\alpha }_3{\widehat{u}}_0^2+{\alpha }_2{\widehat{u}}_0{\widehat{u}}_1& =& 0,\hfill \\ \hfill |-{\alpha }_2{\widehat{u}}_0-{\alpha }_1{\widehat{u}}_1-{\beta }_2|& \le & g_2,\hfill \\ \hfill g_2\left|{\widehat{u}}_1\right|-{\beta }_2{\widehat{u}}_1-{\alpha }_1{\widehat{u}}_1^2-{\alpha }_2{\widehat{u}}_0{\widehat{u}}_1& =& 0\hfill \end{array}$$

(51)

with

$${\alpha }_1=4EJ\left(\frac{1}{l-\widehat{x}}+\frac{1}{\widehat{x}}\right),{\alpha }_2=6EJ\left(\frac{1}{{(l-\widehat{x})}^2}-\frac{1}{{\widehat{x}}^2}\right),{\alpha }_3=12EJ\left(\frac{1}{{(l-\widehat{x})}^3}+\frac{1}{{\widehat{x}}^3}\right)$$

and

$${\beta }_1=\frac{fl}{2},{\beta }_2=\frac{f}{12}({\widehat{x}}^2-{(l-\widehat{x})}^2).$$

We solve the system (51) by solving only the equations (which consist of the elimination of absolute values, dividing by the ${\widehat{u}}_0$ and ${\widehat{u}}_1$ to obtain systems of linear equations), and afterwards we choose the solution which satisfies all the equations and inequalities (51).

To use the numerical algorithms, we will need to prepare objects in optimization problem ${\left({\tilde{\mathcal{P}}}_{TM}\right)}_N$ (38). It can be easily shown that in the case of constant $E,J,f\in \mathbb{R}$, matrix $\mathit{A}$ is a block tridiagonal Toeplitz matrix and together with vector $\mathit{b}$, it can be written in a form

$$\mathit{A}=EJ\left[\begin{array}{ccccc}\mathit{H}& {\mathit{G}}^T& & & \\ \mathit{G}& \mathit{H}& {\mathit{G}}^T& & \\ & \ddots & \ddots & \ddots & \\ & & \mathit{G}& \mathit{H}& {\mathit{G}}^T\\ & & & \mathit{G}& \mathit{H}\end{array}\right],\mathit{b}=f\left[\begin{array}{c}h\\ 0\\ ⋮\\ h\\ 0\end{array}\right]$$

(52)

with

$$\mathit{H}=\left[\begin{array}{cc}\frac{24}{h^3}& 0\\ 0& \frac{8}{h}\end{array}\right],\mathit{G}=\left[\begin{array}{cc}-\frac{12}{h^3}& -\frac{6}{h^2}\\ \frac{6}{h^2}& \frac{2}{h}\end{array}\right].$$

To provide a specific set of data for benchmark, we examine algorithms on the example of a steel beam (E = 2.15 × 10${}^{11}$ Nm${}^{-2}$) of the length $l=1$ m. The point obstacle with given friction is given by $\widehat{x}=8/10$ m with friction values $g_1={10}^2$ N and $g_2=5\times {10}^4$ N. We suppose the load function f = 50,000 N. The beam has a rectangular cross section of height $v=0.02$ m and width $s=0.02$ m. The analytical solution of this problem is presented in Figure 4.

We start our examination with number of used elements as $n=20$ and analyze the influence of regularization parameter on absolute error of the solution of regularized problem. We introduce the error measure by

$$\mathrm{err}=\frac{1}{n+1}\sqrt{\sum _i^{n+1}\left(\overline{u}\left(x_i\right)-u_{\epsilon }^{*}\left(x_i\right)\right)},$$

where $\overline{u}$ is analytical solution and $u_{\epsilon }^{*}$ is numerical solution computed by solving (43) with regularization parameter $\epsilon$. This formula measures the relative difference between the solutions in used nodal points. The results are presented in Figure 5.

In the case of the dual problem, the absolute error decreases naturally with the increase of the discretization density, see Figure 6. The main bottleneck of this approach is the computation of the inverse.

Benchmark 2

In this benchmark, we consider the same problem parameters as in the case of Benchmark 1 except for a shape of cross section and a load function. We consider a solid circular cross section of radius $r=0.01$ m and a load function

$$f\left(x\right)=f_0+(f_1-f_0)x$$

with given parameters $f_0,f_1\in \mathbb{R}$. Obviously, f is a linear function with $f\left(0\right)=f_0$ and $f\left(l\right)=f\left(1\right)=f_1$.

Computing the corresponding integrals, it is easy to derive the components of a vector of right–hand sides

$$b_{2i}=\frac{h^3}{15}(f_1-f_0),b_{2i-1}=h(f_0+ih(f_1-f_0)),i=1,\dots ,n.$$

We set $f_0={10}^5$ N and solve the problem for various choice of $f_1$, see Figure 7 (left). To compare the regularization method and the dual problem method, we solve the problem with $f_1=1.5\times {10}^5$ N for increasing the problem dimension n. The computational time is demonstrated in Figure 7 (right). To decrease the small oscillations of the results (especially in the case of short time), we compute all problems 100 times and present the average computational time.

Benchmark 3

In our last benchmark, we increase the number of internal point obstacles with given friction coefficients to 4 equidistantly distributed on $[0,1]$. We consider the coefficients of given sliding and swivel friction the same for all points $g_1={10}^4$ N and $g_2=5\times {10}^4$ N. We consider material constants defined in Benchmark 1 and the solid circular cross section and load function defined in Benchmark 2. The solution provided by the dual problem approach is presented in Figure 8 (left). In the case of 4 obstacle points, the number of corresponding Lagrange multipliers $\mathit{\lambda }$ is 8, $g={[g_1,g_2,g_1,g_2,g_1,g_2,g_1,g_2]}^T\in {\mathbb{R}}^8$ and matrix $\mathit{B}$ (45) has 8 rows. The dimension of dual problem (49) is 8; however, the most time-consuming operation remains the inverse of the Hessian matrix, see Figure 8 (right).

3. Discussion

For the demonstration, we plot the largest condition number from nine systems, which has to be solved, see Figure 5. The exact source of the problem is the term $1/\epsilon$ added to diagonal elements of the Hessian matrix of corresponding regularized problem, see Table 1 for exact equations.

Figure 5 demonstrates the typical problem with regularization: the decrease of regularization parameter decreases the error of solution; however, it consequently increases the condition number of the system matrix $\widehat{\mathit{A}}$. The additional term from the gradient of approximated function (Table 1) is shifting some of the diagonal values of the original SPD matrix $\mathit{A}$. Using the Gershgorin circle theorem [20], we can state that with decreasing $\epsilon$, the corresponding centers of Gershgorin discs (which includes eigenvalues) are converging to infinity, keeping the radius of disks (sum of absolute values of non-diagonal elements in the corresponding row) constant. Consequently, some of the eigenvalues are converging to infinity. As other eigenvalues remain unchanged, the condition number of the modified matrix $\widehat{\mathit{A}}$ is converging to infinity. The system becomes ill-conditioned.

The speed of convergence of iterative numerical algorithms for solving SPD systems (e.g., Conjugate Gradient method [18]) depends crucially on condition number. On the other hand, if we decide to solve the problem using the direct method (e.g., Gauss elimination method), we hit the finite limits of arithmetical operation. To be more specific, at some point in the process, we have to divide by large number on a diagonal. We can conclude that with decreasing $\epsilon$, it is harder to solve the system precisely.

On the other hand, the dual algorithm is not using any regularization, and therefore it is not necessary to tune this parameter to obtain the optimal ratio between the condition number of the corresponding matrix and the error of approximation. However, this approach has a big price: the assembly of the dual problem. In dual problem (49), one has to compute the inverse of the Hessian matrix, which typically scales with $\mathcal{O}\left(n^3\right)$. In our benchmarks, the Hessian matrix is sparse with a suitable structure to obtain scaling $\mathrm{O}\left(n^2\right)$, see results from Benchmark 1 presented in Figure 6 (right), as well as solution of Benchmark 3 presented in Figure 8 (right). In comparison to the regularization method, the computation of inverse is at some point the bottleneck of the dual problem approach, and the method is slower than the regularization approach, see Figure 7 (right). It is necessary to mention that there exist techniques to avoid the computation of the inverse, e.g., the factorization or the matrix-free implementation [21]. Additionally, the reconstruction of primal solution from dual solution (48) can be performed as a solution of a system of linear equations instead of the multiplication with a dense inverse matrix. This is the reason why (using our naive implementation) the reconstruction step scales with $\mathcal{O}\left(n^2\right)$ in Figure 8.

It is necessary to highlight that we did not solve Benchmark 3 with the regularization method. For multiple obstacle points, the number of corresponding equations (43) would be equal to the number of all combinations of Table 1 constructed for each obstacle point. This number scales exponentially with the number of obstacle nodes. For instance, in the case of 4 points, we would have to solve $9^4$ systems of linear equations (43). On the other hand, the dual problem is a much more efficient approach: increasing the number of obstacle points increases the size of the dual problem linearly.

4. Conclusions

In this paper, we presented both the theoretical and practical aspects of solving the bending problem of the beam: starting from variational formulation, discussing the discretization by Finite Element Method, and ending with numerical methods for solving the corresponding optimization problem. Our numerical experiments revealed the advantages and disadvantages of the two most basic techniques for dealing with the non-differentiable friction term in objective function: the regularization method and the dual problem approach.

The performance and efficiency of the regularization method depend on the a priori chosen parameter. If this parameter is large, the approximation of the original function is not sufficient, and consequently the solution of the approximated problem can be wrong. If the parameter is small, the corresponding problem is ill-conditioned and hard to solve. Additionally, in the case of the problem with more obstacle points, this approach becomes practically inapplicable.

The method of dual problem transforms the problem onto the problem of Lagrange multipliers. The size of the problem is equal to the number of defined friction conditions in the original problem. However, the transformation requires the computation with the inverse of the Hessian matrix of the original problem.

Nevertheless, the explicit computation of the inverse matrix can be avoided using smart implementation techniques. This will be the topic of our future research.