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
hotblade robotic cutting, a new robotic manufacturing method that provides
highspeed 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 twodimensional planeis given [2].
Let an arbitrary curve be parameterized as γ(t)
= (x(t); y(t)), t [0;
t_{1}], 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 arclength 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 [13] 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 fastpaced 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
everincreasing role in stateoftheart 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.
Obviously,
the finite difference scheme can be used to give numerical implementation for
the variational problem (13, 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 subintervals 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 3diagonal band structure
and is shown in the Table 1 together with the right hand vector {B}.
Table
1
0

1

2

...

i1

i

i+1

...

n


B

1

2

1



0





θ_{0}

2

1

2



0






...





...






i1




2

1






i

0

0

...

1

2

1

...

0


0

I+1





1

2





...





...






n





0



2


θ_{n}

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 θ_{0} and θ_{n} are
given, we can calculate θ_{1 }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 θ_{1 }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 coordinates of collocation points a_{i} by the
discreet analogue of expression (1)
(12)
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 θ_{0}^{ }and θ_{n }in
degrees. Each sample includes the family of 60 curves obtained with steps
Δθ_{0}^{ }and
Δθ_{n} of
boundary angles. In order to create each family of curves we initially apply
formula (10) to calculate θ_{1 }angle
for each curve. Then we can calculate coordinates of colocation 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 webpage. The program list is presented below.
<!DOCTYPE html>
<html>
<head>
<meta charset=utf8 />
<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 <
nn1; 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>
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 colinear vectors that makes the elastica closed.
It is consistent with the similar conclusion in works [11, 12].
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 nonconvexity,
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.
Russian Foundation for Basic Research has
supported this work under grant # 170700543.
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 PlateBall 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
HotBlade Cutting, in Advances in
Architectural Geometry 2016, vdf Hochschulverlag AG an der ETH Zürich, DOI
10.3218/37784, ISBN 9783728137784, 2016, http://vdf.ch/advancesinarchitecturalgeometry2016.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, SpringerVerlag, 1989.
9.
Christian
Grossmann & HansGö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