Fractional Mathematical Operators and Their Computational Approximation

Usual applied mathematics employs three fundamental arithmetical operators: addition, multiplication, and exponentiation. However, for example, transcendental numbers are said not to be attainable via algebraic combination with these fundamental operators. At the same time, simulation and modelling frequently have to rely on expensive numerical approximations of the exact solution. The main purpose of this article is to analyze new fractional arithmetical operators, explore some of their properties, and devise ways of computing them. These new operators may bring new possibilities, for example, in approximation theory and in obtaining closed forms of those approximations and solutions. We show some simple demonstrative examples.


Introduction
Science mostly hinges on mathematics which has been built on three fundamental operations: the addition, the multiplication, and the exponentiation.We have become extremely used to them.However, the rigidity of these preset tools may at times make it difficult to find ways and solutions for nature problems that are not easily expressed in terms of traditional operators.For example, transcendental numbers are said not to be attainable via algebraic combination with these fundamental operators, which are a long standing human invention.At the same time, simulation and modelling have to deal with some situations in which the only solution is to compute an expensive numerical approximation of the exact solution, which frequently cannot be written in closed form in terms of traditional operators.
The purpose of this work is to present and apply some operations beyond the rigidity of the known fundamental operations, such as an intermediate operation between addition and multiplication.The ideas have some parallelism with fractional derivatives [1] which, for example, allow the accurate characterization of viscous materials with just two parameters [2,3] rendering linear equations.We will also see possible applications of the presented operators in constitutive modelling.
Consider the classification described in Table 1.
As observed in Table 1, in the last column, exponentiation is expressed in terms of a higher operation.This operation is called tetration and is commonly considered as the first hyperoperation because it has a higher rank ( = 4) than exponentiation ( = 3).Preliminary and fundamental works in this line have been performed by Peano [4] with his Peano Arithmetic axioms, Ackermann [5] with the nonprimitive recursive functions, and Goodstein [6] with the hyperoperational sequence, among others.
However, intermediate-rank operations are also possible and may be useful.The problem we address is how to compute a hyperoperation of rank  = 3/2 (from now on,  3/2 ).To this respect, this was already an open question formulated by Rutsov and Romerio [7] or Williams [8].We will first set the proper environment, namely, the arithmetical continuous spectrum.Some definitions are useful to understand some properties that are going to be necessary to compute  3/2 , as well as to establish a comfortable notation for our purposes.

The Arithmetical Continuous Spectrum
The arithmetical continuous spectrum focuses on the idea of having a continuous spectrum of arithmetical operations as it occurs with real numbers.Here we note that there are many possibilities to define hyperoperations of intermediate ranks.
The reader is referred to the explanation and discussion in [9].The fractional operators that are herein defined have the purpose of seeking improved and more compact interpolations in engineering fields in general and in constitutive modelling in particular.
In order to connect operations of different ranks, the environment set of hyperoperations H must be previously defined.For this purpose, a set of definitions and lemmas will facilitate the task of establishing fractional hyperoperators and their computational approximations.Definition 1.Let H be the set that contains every binary hyperoperation of rank ,  ∈ R. H  ⊂ H is the space of integer rank hyperoperations, and H  ⊂ H is the space of fractional rank hyperoperations; that is, H  ⊂ H  ⊂ H.
The suggested notation described in Table 2 will be employed in order to extend a universal notation for hyperoperations of rank  ∈ R. Consider , ,  ∈ R, which resembles the notion of base, exponent, and result, respectively.
The names of R-and L-inverses remind those of root-kind inverse and log-kind inverse, but note that the notation is also such that the bar is slanted right for the R-inverse operator and left for the L-inverse operator.Some trivial particular cases are those for  = 1 and  = 2

Operational nature Notation Operation
Example For instance, in the case of rank  = 3, the relations between the base , the exponent , and the result  are as follows: Definition 2 (hyperoperational mean).All rank  arithmetical operation in H has a corresponding hyperoperational mean, which is defined in the following way: This definition is the natural extension of ranks  = 1 (arithmetical mean) and  = 2 (geometrical mean), The sequence can be run either forwards (direct values) in this way, or backwards (inverse values), in this way Generally, RS is based on the relation that adjacent terms satisfy.This relation is the hyperoperational mean If solved for  +1 the forwards direction of RS is obtained.On the contrary, if solved for  −1 it yields the expression of backwards direction According to Definition 5, an example of rank  = 1 is shown.RS 1 (3, ) satisfy arithmetical mean,  1 , As another example, consider this time rank  = 3. RS 3 (3, ) must satisfy the corresponding hyperoperational mean,  3 , It can be checked that this makes recursive sequences, RS, coincide up to rank  = 3,  ∈ N, with the corresponding hyperoperational sequences, HS.In turn, they also coincide with the Goodstein function if the neutral elements are the same for the corresponding rank of the Goodstein function; that is, (2,4) with , ,  ∈ R,  ̸ = 1.
An example of this property is shown in Figure 1.According to Definition 5, Therefore, Lemma 10 (Tetration Identity for natural ranks).All arithmetical operations in H with natural rank  satisfy 2 ⋅  2 = 4, ∀ ∈ N.
In the case of noninteger ranks Lemma 10 must be formulated as a condition (hypothesis) to be imposed.

Lemma 12 (about arithmetical duality on H). Arithmetical duality is not an inherent property of H.
Proof.From Definition 5, the sequence RS  can be expressed recursively, From Lemma 9, (25) A counter example is enough to prove the lemma.For rank  = 4, (26) Let it be assumed that the lemma is not true and therefore  3 =  ⋅ In other words, Therefore,   =  ⋅

𝑟+1
is not a general rule for a RS of rank .That is, arithmetical duality is not inherent to H.
This package of definitions, conditions, and lemmas enables the proper environment for H and for reaching an approximation of the fractional hyperoperation  3/2 .
Williams [8] suggests a creative technique, using theory of functionals.He is even able to provide discrete numerical values.Unfortunately, William's method does not guarantee a unique solution.Campagnolo et al. [10] also suggest in certain way the possibility of a continuity in hyperoperational hierarchy.
In order to obtain a computable approximation of  3/2 , the Gaussian mean shall be the starting point of this work.The reason behind this is the close relationship between an hyperoperational mean and its direct hyperoperation, as observed in previous lemmas and definitions.

Approximation of 𝐻 3/2
The arithmetic-geometric mean agm(, ) is computed recursively (with very fast convergence) as In this approximation of  3/2 (from now on, simply  3/2 ), the arithmetic-geometric mean, agm, will be the mean of  Therefore, the neutral element  3/2 is readily obtained as The neutral element  3/2 is easy to test as observed in Figure 2.
If setting now  = agm(, ) in (35), an explicit and computable expression of the approximation of  3/2 is obtained, In Figure 3   Equation (38) facilitates the computation of an approximation of  5/2 in the same way multiplication may be computed by addition, In Figure 4 there are examples of hyperoperations of the form  =  ⋅  2; namely,  =  ⋅ 2, and the result is compared to the sum, product, and exponentiation.Some properties of  3/2 are as follows: (i) Its consistency and uniqueness are guaranteed by (38).(ii) It is consistent with its mean definition (iv) It has a neutral element; see (37).
(v) It is not associative (although it produces values of the same magnitude).(vi) It is not distributive with addition (vii) It is not distributive with multiplication

Applications
The existence of fractional hyperoperations enlarges the arithmetical tools of mathematics.Therefore, although these operations are not very intuitive because we are used to traditional ones, it is apparent that fractional hyperoperations may have an impact on applications that rely on every discipline that needs mathematics, particularly interesting in the modelling of physical processes.In fact, as we show below, they may be useful in curve fitting in constitutive modelling of solids [11][12][13], which was one of our interests in this type of operators.In this section we show some applications of the presented framework in this field.In the applications shown, the agmo operator was needed in order to uniquely compute these intermediate hyperoperations.The algorithm we employed is shown as follows.
where  are the true (logarithmic) strains and  are the true (Cauchy) stresses and work-conjugate of the logarithmic strains if the volumetric elastic strains are neglected (they are typically orders of magnitude smaller).
This compact expression in very close agreement with the experimental results has been obtained employing a mixed rank-parameter fitting procedure.When the rank is fixed, the procedure is simpler although the result is not so compact.For example, a possible systematic solution would be the generalization of the Vandermonde matrix to any hyperoperation of rank .Vandermonde matrix is useful for interpolation.In order to explain how the usual Vandermonde matrix works, let us consider the case of a 2D interpolation of  experimental points with coordinates (  ,   ), so that, given any  as input, its interpolated image  is obtained.Then, the Vandermonde matrix is built upon the experimental points.It is the representation of a system of equations, in which the  incognita are monomial powers of , of a polynomic expression of grade , for each point (  ,   ) So, in matrix representation, However, a well-known drawback of this method is that the more the terms in the polynomial are included, the bigger the condition number of the resulting matrix is.The reason behind this is the behaviour of    as the degree  of the polynomial increases.
However    can be reexpressed in hyperoperational notation as , the condition number would become lower as the number of terms, , increases.The more the accuracy in the interpolation fitting is required, the more the terms that the polynomial must include are, and the bigger the differences of using  ⋅ 3/2  and  ⋅ 3  become.Furthermore, the method can be adjusted as desired just by setting the hyperoperational rank  of the hypermonomials,  ⋅  .As an illustrative example, in Figure 6 a modified Vandermonde with hyperpolynomial (instead of the usual polynomial) expansion of rank  = 3/2 is applied in order to fit points extracted from the well-known classical Treloar [15] experiment of simple tension, from a sample cut of a simple sheet of vulcanized natural rubber.In particular, the extension has the form described in where the coefficients   come from the analogous vector in (49), as described in (49) As commented, as the number of terms included in the expansion increases, the differences in the condition number between the usual Vandermonde matrix and the modified one are larger.This makes the use of fractional hyperoperations in the Vandermonde method particularly interesting when dealing with a high number of terms in the Vandermonde expansion.An important advantage of hyperoperations in the Vandermonde method is that the rank  of the hyperoperation does not necessarily have to be  = 3/2, but a free parameter instead, in order to best fit the interpolation behaviour.In Figure 6(c), for 15 terms, the hyperoperational variant ( = 3/2) already offers better fitting than the usual Vandermonde one ( = 3).For higher number of terms, the condition number of the usual Vandermonde matrix becomes so unstable that interpolation application is no longer viable.In Table 3, the evolution of the condition number with the inclusion of new terms has been represented from  = 15 terms up to  = 20 terms.
Consider a virtual experiment with the (presumably unknown) result being of the form of a logarithm (i.e., y = log(x), or y = x \ 3  in hyperoperational notation).The result of the experiment is obtained for a vector of random  values.The aim of this emulated experiment is to contrast the value of the two Vandermonde options of the previous example, with the known exact value for the interpolated points; see Figure 7.
For a particular example, pseudoexperimental scattered points from the analytical function  = log() emulate an experimental event.In this simple numerical experiment we wish to obtain a proper fit for point in the interval  ∈ [1.3 : 0.1 : 7].This numerical experiment has been performed with two different sets of random numbers in the interval [1.3 : 7].A summary of the results is presented in Tables 4 and 5, where we used  = 15 terms in both cases with precision of 10 −30 .It is clearly seen that the precision using  = 3/2 is many orders better than that using the classical Vandermonde interpolation.In order to further explore the applicability of the presented operators to curve fitting, we present another example.This time the number of terms used in the Vandermonde method has been increased up to 1000, and the pseudoexperimental curve employed in this example is () = sin().As it could have been expected, the typical Vandermonde approach fails with such number of terms.In fact, in this test, nonnumerical NaNs have been returned by MATLAB when employing the usual Vandermonde approximation due to the extremely high condition number reached, so the required numerical precision to accomplish the task cannot be attained using this approach.Nevertheless,  3/2 Vandermonde breaks this constraint even with practical machine precision of 1 − 15, as it is shown in Figure 8.Note that the 1000 points are the points employed to build the  3/2 Vandermonde hyperpolynomial, not the interpolated points.Hence, the proposed operators clearly bring new possibilities in the use of polynomial approximations by Vandermonde-type schemes.
The last numerical example we show here is to show some additional properties of this kind of interpolation: smoothing  of noisy experimental data.This is a very important problem in many engineering fields, as for example, in statistics [16] or in the case of curve and surface fitting [17].It is also of practical importance in fitting experimental data from many materials, as polymers and biological tissues.In this case we compare the proposed interpolation method with spline interpolation techniques [17].Spline interpolation is used in the mentioned fields because, among other reasons, the usual Vandermonde approximation loses accuracy when a big amount of points is taken into account.However, splines have their own well-known problems related to large oscillations when noisy experimental data is presented; see, for example, [16][17][18] and references therein.However,  3/2 interpolation is global in nature and uses a householder transformation.The interpolation is naturally not forced to pass through all the points in this case.The smoothness of the approximation can be controlled with the order of the hyperpolynomial, a desirable characteristic when dealing with noisy data.For the purpose of this example, the Treloar experimental data for rubber from the previous example have been distorted with a given amount of noise.We show the altered data presented to the approach in Figure 9.In this figure we also show, for comparison purposes, a piecewise cubic spline interpolation and the proposed  3/2 Vandermonde interpolation.It is seen that the proposed global interpolation naturally filters the experimental noise, whereas the spline interpolation would  need specific smoothing techniques.Furthermore, the global nature of the proposed interpolation may be interesting when addressing some desired properties of stored energy functions of hyperelastic materials as polyconvexity.The inherent local nature of the splines used in the WYPIWYG approaches [11,12] makes in principle the analysis of these properties more difficult.

Figure 4 :
Figure 4: Approximation of the result of some hyperoperations   of the form  =  ⋅  2.

Figure 5 :
Figure 5: Modelling curve based on experimental results from alloy C51100 under uniaxial true stress versus strain after necking.

Table 1 :
Hierarchical classification of fundamental operations.
Definition 5 (recursive sequence, RS).Let a recursive sequence RS  (, ), of rank , base , and neutral element   , be the sequence built by its corresponding mean,   , and the following pair of starting values: Figure 1: Arithmetical reflectivity illustrative example.Definition 6 (arithmetical reflectivity).Let HS  (, ) be a hyperoperational sequence of rank .Then it is said that hyperoperation   owns arithmetical reflectivity if and only if 7 (arithmetical duality or duality RS-HS).It is said that a hyperoperation   has arithmetical duality when its recursive sequence, RS  , and its hyperoperational sequence, HS  , coincide, RS  = HS  .Let RS  (, ) be a recursive sequence of base , with  ∈ R.Then, the precedent term,   , of {,  ⋅ Let RS  (, ) be a recursive sequence of base .Then, the first three terms of RS  (, ) are the term   being the right neutral element of  +1 .Proof.The right neutrality of   with  ≥ 2 is imposed to the value of   = 1 due to multiplicative recursion,  ⋅ +1 2}, is the neutral element of RS  (, ), which works as right identity element of     (  ,  ⋅ +1 2) = .(20)Inthecase of rank  ∈ N, Conjecture 8 becomes a lemma.Lemma 9 (RS starting set).+1 1 = .

Table 4 :
Mean and standard deviation of difference betweenVandermonde interpolation and the exact solution.

Table 5 :
Mean and standard deviation of difference betweenVandermonde interpolation and the exact solution.