1.        Introduction

The elastica task is formulated as follows: what shape do we obtain when an elastica is bent? Here the term elastica means an ideal infinitesimally thin elastic rod without stretching on a plane. Leonard Euler essentially solved this problem in 1744 but similar problems are still being studied by researchers [1]. Euler obtained differential equations for stationary rod configurations and described their possible qualitative types. These configurations are called Euler elastica [2, 3].

Since then up to now, the problem solution related to Euler elastica energy has had a wide range of applications in computer graphics, engineering geometry, in mechanics, engineering, control theory, approximation theory, molecular biology, nanotechnology and other art and industrial areas. Euler elastica energy functional have also wide application in computer vision and image processing to restore the image distorted by multiplicative noise [4].

Work [5] presents a statistical application of Euler elastica to create a model of DNA and molecules of large polymer set. The authors of work [6] describe a new method for creating double curvature surfaces for architectural design by means of Euler elastica. This method provides a direct design of the workflow using a hot-blade robotic cutting, a new robotic manufacturing method that provides high-speed production of double curvature molds. Additionally to perfect the structure, when arbitrary surfaces are converted to a geometry suitable for cutting with a hot blade, the method allows architects and designers to explore the unique architectural potential of this approach in manufacturing. Thus, the development of methods for constructing Euler elastica is of inescapable interest.

Let us briefly consider the formulation of the problem of Euler elastica form finding. Assume that Cartesian system (x; y) in two-dimensional planeis given [2]. Let an arbitrary curve be parameterized as γ(t) = (x(t); y(t)), t  [0; t₁], and let its endpoints have coordinates a_i = (x_i; y_i), i = 0, 1. Denote by θ the angle between the tangent vector to the curve γ and the positive direction of the axis x. Further, let the tangent vector at the endpoints of γ have coordinates v_i = (cos θ_i; sin θ_i), i = 0, 1, see Fig. 1.

 

Fig. 1. Stationary configuration of elastic rod

 

Then the required curve γ(t) = (x(t); y(t)) is determined by a trajectory of the following control system:

(1)

 

(2)

(3)

For an arc-length parametrized curve, the curvature is, up to sign, equal to the angular velocity:  whence we obtain the functional [3]

(4)

The problem of functional (4) analytical solution regarding the relations (1) - (3) is a difficult mathematical issue. Analytical approaches are often based on the use of various mathematical procedures leading to Jacobi functions, that is, to complex combinations of elliptic functions [1-3] and elliptic integrals of the 1st and 2nd kind [7]. Such representation is extremely inconvenient for visualizing Euler elastica by computer means. On the other hand, since the end of the last century the fast-paced electronic computing facilities have provided a solid basis for applying numerical methods to solve various problems [8]. The efficient numerical solution of differential and equations plays an ever-increasing role in state-of-the-art technology. Both the demand and the computational power available from current computer hardware have stimulated the rapid development of numerical methods for this kind of equations [9], [10]. The discrete numerical approach to solve the Euler elastica problem is proposed in this work.

2.        Finite difference formulation

Obviously, the finite difference scheme can be used to give numerical implementation for the variational problem (1-3, 4) solution. When a finite difference method (FDM) is used to treat numerically a differential equation, the differentiable solution is approximated by some grid function, i.e., by a function that is defined only at a finite number of points of collocation that lie in the studied domain and its boundary. Each derivative that appears in the differential equation has to be replaced by a suitable divided difference of function values at the chosen collocation points [9]. Thus, FDM uses the approximation of derivatives by the set of algebraic forms. This converts the differential equation into an approximated algebraic problem, which can be solved by a finite computational procedure.

Let the required curve have length L, then integral (4) is given in (0; L) interval. To replace integral (4) by sum of n integrals, we subdivide further this interval into n sub-intervals with length l = L / n.

(5)

 

Fig. 2. Discreet model of elastic rod

 

Assume that boundary conditions of two endpoints are ,  and the coordinates x_n, y_n of right hand endpoint are unknown. If we neglect curve linearity of each curve interval (see Fig. 2), we can replace function  by the following finite difference

(6)

This allows us to obtain the presentation of energy integral (4) by the following approximate finite sum

(7)

To find unknown tangent angles we should consistently take the derivatives of function J by θ_i (see form (7)). Finally we can obtain the resolving system.

(8)

The sparse symmetry matrix [A] has a 3-diagonal band structure and is shown in the Table 1 together with the right hand vector {B}.

 

Table 1

The structure of matrix [A] and vector {B}

 

It is easy to notice that there is no need to use special procedures to solve the equation system in Table 1 because each tangent angle is the arithmetic mean of two neighboring angles as can be seen below

(9)

This circumstance provides the possibility to apply a very efficient and fast computational scheme. If endpoint tangent angles θ₀ and θ_n are given, we can calculate θ₁ angle by the formula that as appears from the recurrent application of expression (9) starting from endpoint n in the opposite direction of the elastic rod

(10)

Since θ₁ angle is known then the entire angles can be calculated step by step by the following recurrent formula

(11)

Finally, expression (11) allows obtaining the entire angle vector {θ}. Further, we can calculate nodal co-ordinates of collocation points a_i by the discreet analogue of expression (1)

(12)

3.        Numerical Samples

There are some numerical samples in this chapter. The set of samples is obtained for elastica of the initial rectilinear shape, restrained at the left endpoint and for the following parameters: rod length L = 500 pixels, the number of subintervals n = 200. All the samples are calculated with different values of endpoint angles θ₀ and θ_n in degrees. Each sample includes the family of 60 curves obtained with steps Δθ₀ and Δθ_n of boundary angles. In order to create each family of curves we initially apply formula (10) to calculate θ₁ angle for each curve. Then we can calculate co-ordinates of co-location points by expressions (11) and (12). A very similar and efficient program based on the numerical scheme described in this paper and HTML5 / Javascript languages has been created in the form of a web-page. The program list is presented below.

<!DOCTYPE html>

<html>

<head>

<meta charset=utf-8 />

<title>Draw Euler elastica</title>

</head>

<body>

<canvas id="DemoCanvas" width="1200" height="800"></canvas>

<script>

let cx = document.querySelector("canvas").getContext("2d");

var x_0 = 400;

var y_0 = 300;

var nn = 200;

var len = 500/nn;

for (var j = 0; j < 60; j++) {

     var alfa_n = (180 - 5 * j) * Math.PI/ 180.0 ;

     var alfa_0 = (0 + 5 * j) * Math.PI/ 180.0 ;

     var alfa_1 = (alfa_n + (nn - 1) * alfa_0) / nn;

     var alfa_i = alfa_0;

     var alfa_i_1 = alfa_1;

     var alfa_i_2 = alfa_0;

     var x_beg = x_0;

     var y_beg = y_0;

     cx.beginPath();

     cx.moveTo(x_0, y_0);

     for (var i = 0; i < nn-1; i++) {

         cx.lineTo(x_beg + len * Math.cos(alfa_i), y_beg

                          + len * Math.sin(alfa_i));

         x_beg = x_beg + len * Math.cos(alfa_i);

         y_beg = y_beg + len * Math.sin(alfa_i);

         alfa_i = 2 * alfa_i_1 - alfa_i_2;

         alfa_i_2 = alfa_i_1;

         alfa_i_1 = alfa_i;

     }

  cx.stroke();}

</script>

</body>

</html>

 

a)
b)
c)
d)
e)

Fig. 3. The sample families of elastica with different boundary conditions

The curve families are graphically presented in Figure 3. It is easy to spot that in all the samples the curves converge into an ideal circle the way boundary tangent vectors converge into co-linear vectors that makes the elastica closed. It is consistent with the similar conclusion in works [11, 12].

4.        Conclusion

The paper describes a fast and efficient method that is equivalent to the minimization of Euler elastica energy functional. It is an extremely difficult problem to find the task solution by minimizing the energy functional due to its non-convexity, nonlinearity and higher order with derivatives. The method is fast to produce the solution and simple to implement. It seems quite encouraging that this kind of algorithm could have reliable application in a number of real industrial problems related to scientific visualization and computer graphics.

5.        Acknowledgments

Russian Foundation for Basic Research has supported this work under grant # 17-07-00543.

Reference

1.     Yu.Sachkov, Maxwell strata in Euler’s elastic problem, Journal of Dynamical and Control Systems, Vol. 14 (2008), No. 2 (April), 169–234. Available at: arXiv:0705.0614 [math.OC], 3 May 2007.

2.     V. Jurdjevic, The Geometry of the Plate-Ball Problem, Arch. Ration. Mech. Anal. 124, 305–328 (1993).

3.     V. Jurdjevic, Geometric Control Theory, Cambridge University Press, 1997.

4.     Khan MA, Chen W, Ullah A, Ji L (2018) Euler's elastica and curvature based model for image restoration. PLoS ONE 13(9): e0202464. https://doi.org/10.1371/journal.pone.0202464

5.     Matsutani, Shigeki. (2010). Euler’s elastica and beyond. Journal of Geometry and Symmetry in Physics. 17.

6.     David Brander et al, Designing for Hot-Blade Cutting, in Advances in Architectural Geometry 2016, vdf Hochschulverlag AG an der ETH Zürich, DOI 10.3218/3778-4, ISBN 978-3-7281-3778-4, 2016, http://vdf.ch/advances-in-architectural-geometry-2016.html

7.     K. Mattiasson, Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals, Internat. J. Numerical Methods Engrg. 17 (1981) 145–153.

8.     D.F. Lawden, Elliptic functions and applications, Springer-Verlag, 1989.

9.     Christian Grossmann & Hans-Görg Roos, Numerical Treatment of Partial Differential Equations, Springer Berlin Heidelberg New York, 2005. ISBN: 9783540715825

10. Donald Greenspan, Discrete Numerical Methods in Physics and Engineering, ACADEMIC PRESS, INC., New York, 1974.

11. Piotr Szablewski, Ryszard Korycki  Shape Determination of Elastica Subjected to Bending by Means of Displacements FIBRES & TEXTILES in Eastern Europe 2016, Vol. 24, 6(120)  DOI: 10.5604/12303666.1221748

12. Yu. L. Sachkov, Closed Euler elasticae, Differential equations and dynamical systems, Collected papers, Tr. Mat. Inst. Steklova, 278, MAIK Nauka/Interperiodica, Moscow, 2012, 227–241