# The Tensor Pade´-Type Approximant with Application in Computing Tensor Exponential Function

Authors: Chuanqing Gu; Yong Liu
Journal: Journal of Function Spaces
 (2018)

---

## Abstract

Tensor exponential function is an important function that is widely used. In this paper, tensor Pade´-type approximant (TPTA) is defined by introducing a generalized linear functional for the first time. The expression of TPTA is provided with the generating function form. Moreover, by means of formal orthogonal polynomials, we propose an efficient algorithm for computing TPTA. As an application, the TPTA for computing the tensor exponential function is presented. Numerical examples are given to demonstrate the efficiency of the proposed algorithm.


---

## Body

## 1. Introduction

Tensor exponential function is an important function that is widely used, owing to its key role in the solution of tensor differential equations [1–4]. For instance, Markovian master equation can be written as tensor differential equations ∂/∂tP(t)=A·P(t), where the probabilities tensor P(t)∈Rn1×n2×⋯×nd [5]. Consider the initial value problem defined by the tensor ordinary differential equation [6, 7](1)Y˙t=AYt,Yt0=Y0,where the superimposed dot denotes differentiation with respect to t and A and Y0 are given constant tensors. The solution to system (1) is Y(t)=exp[(t-t0)A]Y0, and exp(·) is the tensor exponential function.To solve (1), we need to calculate an exponential function about tensor A. Recently, tensor computation, especially eigenvalues of tensors, has attracted attention of many scholars; some important results have been found in the current literature [8–15]. For instance, computation of exp(At) is required in applications such as finite strain hyperelastic-based multiplicative plasticity models [7, 16–19]. Explicitly, for a generic tensor A, the tensor exponential can be expressed by means of its series representation [7] (see p.749).(2)expA=∑n=0∞1n!An,and the preceding series is absolutely convergent for any argument A and, as its scalar counterpart, can be used to evaluate the tensor exponential function to any prescribed degree of accuracy [16]. The computation of (2) is carried out by simply truncating the infinite series with nmax terms: (3)expA=∑n=0nmax1n!An, with nmax being such that 1/nmaxAnmax≤ϵtol.However, the accuracy and effectiveness of the preceding algorithm is limited by round-off and choice of termination criterion [16]. Pade´ approximant has become by far the most widely used one in calculation of exponential function or formal power series due to the following reasons: first, the series may converge too slowly to be of any use and the approximation can accelerate its convergence; second, only few coefficients of the series may be known and a good approximation to the series is needed to obtain properties of the function that it represents [20]. For instance, matrix Pade´-type approximant (MPTA) [21] can be used to simplify the high degree multivariable system by approximating transfer function matrix G(s) that can be expanded into a power series with matrix coefficients, i.e., G(s)=∑i=0∞cisi, where ci∈Cs×t. The key to construct TPTA is to maintain the same order of tensor A with different powers. For this issue, we introduce t-product [22, 23] of two tensors to solve it. In addition, in order to give the definition of TPTA, we introduce a generalized linear functional in the tensor space for the first time.This paper is organized as follows. In Section2, we provide some preliminaries. First, we introduce the t-product of two tensors; then, we show the definitions of tensor exponential function and the Frobenius norm of a tensor. In Section 3.1, we define the tensor Pade´-type approximant by using generalized linear functional; the expression of TPTA is of the form of tensor numerator and scalar denominator; and then we introduce the definition of orthogonal polynomial with respect to generalized linear functional and sketch an algorithm to compute the TPTA. Numerical examples are given and analyzed in Section 4. Finally, we finish the paper with concluding remarks in Section 5.


## 2. Preliminaries

There arise mainly a problem for approximating tensor exponential function. That is how to expandeAt into the power series for order-p (p≥3) tensors. For a Symmetric and second-order tensor A, higher powers of A can be computed by the Cayley-Hamilton theorem [24], but it fails for the order-p (p≥3) tensors. Therefore, we shall utilize the t-product to obtain higher powers of order-p (p≥3) tensors in this section. Firstly, we introduce some notations and basic definitions which will be used in the sequel. Throughout this paper tensors are denoted by calligraphic letters (e.g., A, B), while capital letters represent matrices, and lowercase letters refer to scalars.An order-p tensor, A, can be written as (4)A=ai1i2⋯ip∈Rn1×n2×⋯×np.Thus, a matrix is considered a second-order tensor, and a vector is a first-order tensor [22], for i=1,…,np, denoted by Ai∈Rn1×n2×⋯×np-1, the tensor whose order is p-1 and is created by holding the pth index of A fixed at i. For example, consider a third-order tensor, A=(ai1i2i3)∈R3×3×3. Fixing the 3rd index of A, we can get three 3×3 matrices, namely, 2-order tensor, which are A1, A2, and A3 and with elements (5)A1:a111a121a131a211a221a231a311a321a331,A2:a112a122a132a212a222a232a312a322a332,A3:a113a123a133a213a223a233a313a323a333, respectively.Now, we will define the t-product of two tensors.Definition 1 (see [22, 23]).
LetA∈Rn1×n2×⋯×np. Then the block circulant pattern tensor of A is denoted by (6)circA=A1AnpAnp-1⋯A2A2A1Anp⋯A3⋮⋱⋱⋱⋮AnpAnp-1⋯A2A1, where Ai∈Rn1×n2×⋯×np-1, i=1,2,…,np.Defineunfold(·) to an n1×n2×⋯×np tensor by an n1np×n2×⋯×np-1 block tensor in the following way: (7)unfoldA=A1A2⋮Anp. If A is order-3 tensor, then unfold(A) is a block vector. Similarly, definefold(·) as the inverse operation, which takes an n1np×n2×⋯×np-1 block tensor and returns an n1×n2×⋯×np tensor; then (8)foldunfoldA=A.Definition 2 (see [23]).
LetA be Rn1×n2×⋯×np and B be Rn2×l×n3⋯×np. Then the t-product A*B is the n1×l×n3×⋯×np tensor defined recursively as (9)A*B=foldcircA*unfoldB.Remark 3.
IfA and B are order-2 tensors, then the product “∗” can be replaced by standard matrix multiplication.Remark 4.
Thek times power of A is defined as Ak=A*A*⋯*A (k times); “∗” denotes the t-product.Example 5.
Letting(10)A=10-13102102000-111-12∈R3×3×2,then, from Definition 2, we have (11)A4=A*A*A*A=271113-98237104-11924297-9624296-96-73-4656-94-5243.Remark 6.
One of the characteristic features of t-product is that it ensures that the order of multiplication result of two tensors does not change, whereas other tensor multiplications do not have the feature; that is why we chose the t-product as the multiplication of tensors.The tensor exponential function is a tensor function on tensors analogous to the ordinary exponential function, which can be defined as follows.Definition 7.
LetA be an n1×n2×⋯×np real or complex tensor. The tensor exponential function of t, denoted by eAt or exp(At), is the n1×n2×⋯×np tensor given by the power series (12)eAt=∑k=0∞1k!Atk, where A0 is defined to be the identity tensor I (see Definition 8) with the same orders as A.Definition 8 (see [23]).
Then×n×l1×⋯×lp-2 order-p identity tensor (p>3) I is the tensor such that I1 is the n×n×l1×⋯×lp-3 order-(p-1) identity tensor, and Ij is the order-(p-1) zero tensor, for j=2,3,…,lp-2.By Definition8, we can define tensor inverse, transpose, and orthogonality. However, we do not discuss these works here, as it is beyond the scope of the present work. For the details of these definitions of tensor, we refer to reader to [22, 23, 25] and the references therein.LetA∈Rn1×n2×⋯×np; then the norm of a tensor is the square root of the sum of the squares of all its entries [25]; i.e.,(13)A=∑i1=1n1∑i2=1n2⋯∑ip=1npai1i2⋯ip2.This is analogous to the matrix Frobenius norm. The inner product of two same-sized tensorsA, B∈Rn1×n2×⋯×np is the sum of the products of their elements [25]; i.e.,(14)A,B=∑i1=1n1∑i2=1n2⋯∑ip=1npai1i2⋯ipbi1i2⋯ip.It follows immediately that (A,A)=A2.


## 3. Tensor Pade´-Type Approximant

Letf(x) be a given power series with tensor coefficients; i.e.,(15)fx=A0+A1x+A2x2+⋯+Anxn+⋯,Ai∈Rn1×n2×⋯×np,x∈R.LetP denote the set of scalar polynomials in one real variable whose coefficients belong to the real field R and Pk denote the set of elements of P of degree less than or equal to k.Letϕ:P→Rn1×n2×⋯×np be a linear functional on P. Let it act on t by(16)ϕti=Ai,i=0,1,2,….Then, by the linearity of ϕ, we have (17)ϕ1-tx-1=ϕ1+tx+tx2+⋯+txn+⋯=A0+A1x+A2x2+⋯+Anxn+⋯=fx.


### 3.1. Definition of Tensor Pade´-Type Approximant

Letvnx=b0+b1x+⋯+bnxn(bn≠0) be a scalar polynomial of Pn of exact degree n. In this case, vn(x) is said to be quasi-monic. Define the tensor polynomial Wn-1(x) associated with vn(x) with tensor-valued coefficients, by(18)Wn-1x=ϕvnx-vntx-t.It is easily seen that Wn-1(x) is a tensor polynomial of exact degree n-1 in x. Set(19)v~nx=xnvnx-1,(20)W~n-1x=xn-1Wn-1x-1.Then, the polynomials v~n(x) and W~n-1(x) are obtained from vn(x) and Wn-1(x), respectively, by reversing the numbering of the coefficients. By the procedure given above, the following conclusion is obtained.Theorem 9.
Letv~n(0)≠0; then W~n-1(x)/v~n(x)-f(x)=O(xn).Proof.
Expandingvn(x)-vn(t)/x-t in (18) and applying ϕ yields that (21)W~n-1x=∑l=0n-1∑i=0n-l-1bl+i+1Aixl=∑l=0n-1∑i=0lbn-l+iAixl. Computing v~n(x)f(x), we get (22)v~nxfx=∑j=0nbn-jxj∑i=0∞Aixi=∑l=0∞∑i=0lbn-l+iAixl. Thus, (23)v~nxfx-W~n-1x=∑l=n∞∑i=0lbn-l+iAixl=Oxn.Definition 10.
R n - 1 , n ( x ) = W ~ n - 1 ( x ) / v ~ n ( x ) is called a tensor Pade´-type approximant with order n for the given power series (15) and is denoted by (n-1/n)f(x).Remark 11.
The polynomialvn(x), called the generating polynomial of (n-1/n)f with respect to power series f(x), can be arbitrarily chosen.Remark 12.
The tensor Pade´-type approximant (n-1/n)f possesses the degree constraint, which is caused by its construction process. The constraint implies that the method does not construct tensor Pade´-type approximant of type (m/n) when m is different from n-1.To fill this gap, we define a new tensor Pade´-type approximant by introducing a generalized linear functional.Letϕ(q):P→Rn1×n2×⋯×np be a generalized linear functional on P. Let it act on t by(24)ϕqti=Aq+i,i=0,1,2,….Similarly to what was done for Wn-1(x), we consider the polynomial Wq(x) associated with vn(x), and defined by(25)Wqx=ϕqvnx-vntx-t,q=m-n+1.Set(26)W~qx=xn-1Wqx-1,and define(27)Pmnx=v~nx∑i=0m-nAixi+xm-n+1W~qx,m≥n.Then we have the following conclusion.Theorem 13.
Letv~n(0)≠0; then Pmn(x)/v~n(x)-f(x)=O(xm+1).Proof.
Letfm-n+1 be the formal power series (28)fm-n+1x=∑j=0∞Am-n+1+jxj. Then(29)xm-n+1fm-n+1x=fx-∑i=0m-nAixi.Expanding (25) and using (26) we obtain (30)W~qx=∑q=0n-1∑i=0qbn-q+iAm-n+1+ixq. Computing the product v~n(x)fm-n+1(x), we find that (31)v~nxfm-n+1x=∑q=0∞∑i=0qbn-q+iAm-n+1+ixq. By Theorem 9, one has (32)W~qxv~nx=n-1nfm-n+1x. Then, for m≥n we deduce from (27) and (29) that (33)v~nxfx-Pmnx=v~nx∑i=0m-nAixi+xm-n+1fm-n+1x-v~x∑i=0m-nAixi+xm-n+1W~qx=xm-n+1v~nzfm-n+1-W~qx=xm-n+1∑q=n∞∑i=0qbn-q+iAm-n+1+ixq=Oxm+1.Now, we can achieve(m/n)f(x) by the above procedure, and it will be denoted by Pmn(x)/v~n(x).Definition 14.
R m , n ( x ) = P m n ( x ) / v ~ n ( x ) is called TPTA with order m+1 and is denoted by (m/n)f(x).Algorithm 15 (computeRm,n(x)=Pmn(x)/v~n(x) with vn(x) being arbitrarily chosen).
(1)
Setq=m-n+1 and chose a quasi-monic polynomial vn(x). (2)
Use (19) to compute v~n(x). (3)
ComputeWq(x) and W~q(x) by (25) and (26), respectively. (4)
Substitutev~n(x) and W~q(x) into (27) to compute Pmn(x). (5)
SetRm,n(x)=Pmn(x)/v~n(x).Example 16.
Let(34)fx=10011-102-1221+10211-212-1231x+12110-102-1320x2+⋯=A0+A1x+A2x2+⋯.Now we apply Algorithm15 to compute TPTA of type (2/2) for this example.(1)
Chosev2(x)=x2-2x+4, q=m-n+1=1.(2)
Use (19) to compute v~2(x): (35)v~2x=1-2x+4x2.(3)
By using (25) and (26) we get(36)W1x=x-122x+1x-1x-21-2xx2x-2-x+12x-13x-2x-2,and(37)W1~x=1-x2xx1-x1-2xx-102-2xx-12-x2-2x1-2x.(4)
Substitutingv~2(x) and W~1(x) into (27), we obtain(38)P22x=v~2xA0+xW1~x=3x2-x+12x2-3x2+2x3x2-x+12x2-x+1-x2-1-2x2+x6x2-2x+2-3x2+x-17x2-2x+24x2-x+22x2-x+1.(5)
SetR22(x)=P22(x)/v~2(x). It is easy to verify that (39)P22xv~2x-fx=Ox3.


### 3.2. Algorithm for Computing TPTA

Generally, the precision of TPTA is limited, since the denominator polynomials of TPTA are arbitrarily prescribed. In this subsection, in order to improve the precision of approximation, we propose an algorithm for computing the denominator polynomials and illustrate the efficiency of this algorithm in next section.First, we give the following conclusion.Theorem 17 (error formula).
(40) f x - n - 1 n f x = x n v ~ n x ϕ v n t 1 - t x .Proof.
Note thatϕ is a linear functional on P, only acting on t. From (18) and (20) we deduced that (41)W~n-1x=xn-1Wn-1x-1=xn-1ϕxvnx-1-xvnt1-tx=ϕxnvnx-1-xnvnt1-tx=ϕxnvnx-11-tx-xnϕvnt1-tx=v~nxfx-xnϕvnt)1-tx, and then this error formula holds.In terms of the error formula, it holds that(42)fx-n-1nfx=xnv~nxϕvnt1-tx=xnv~nxϕvnt+vnttx+vntt2x2+⋯=xnv~nxϕvnt+ϕvnttx+ϕvntt2x2+⋯.If we impose that vn(t) satisfies the condition ϕ(vn(t))=0, then the first term of (42) disappears, and the order of approximation becomes n+1. If, in addition, we also impose the condition ϕ(tvn(t))=0, the second term in the expansion of the error also disappears, and the order of approximation becomes n+2, and so on. We indicate that vn(x) depends on n+1 arbitrary constants; however, on the other side, a rational function is defined apart from a multiplying factor in its numerator and its denominator. It implies that (n-1/n)f(x) depends on n arbitrary constants. So let us take vn(t) such that(43)ϕvnttk=0,k=0,1,2,…,n-1.Definition 18.
v n ( t ) in (43) is called an orthogonal polynomial with respect to the linear functional ϕ and (n-1/n)f(x) in (42) is also called a TPTA for the given power series (15) when (43) is satisfied.From (43) we obtain(44)ϕvnttk=∑i=0nbiAi+k=0,k=0,1,2,…,n-1.Let bn=1 in (44); then it follows that(45)∑i=0n-1biAi+k=-Ak+n,k=0,1,2,…,n-1.Forming the scalar product of both sides of (45) with A0,A1,…,An-1, respectively, we get(46)∑i=0n-1biAi+k,Ak=-Ak+n,Ak,k=0,1,2,…,n-1.Denote(47)HnA0=A0,A0A1,A0⋯An-1,A0A1,A1A2,A1⋯An,A1⋮⋮⋱⋮An-1,An-1An,An-1⋯A2n-2,An-1,h→n=-An,A0-An+1,A1⋮-A2n-1,An-1,and call det(Hn(A0)) the Hankel determinant of f(x) with respect to the coefficients A0,A1,…,An-1.Then (46) is converted into(48)HnA0x=h→n,x=b0,b1,…,bn-1T.In the case of TPTA,vn(x) is not arbitrarily chosen any more but is determined by the preceding system. The choice of vn(x) can help to improve the accuracy of approximation, but unfortunately we have not been able to guarantee that the solution of system (48) comes into existence, so far. We only give the following basic theorem about system (48) on the basis of linear algebra.Theorem 19.
The solution of (48) exists if and only if rank(Hn(A0))=rank(Hn(A0):h→n). Moreover, the solution is unique if det(Hn(A0))≠0.Proof.
The proof of the assertion follows from the simple fact that, for a system of linear equations, described byAx=b, where A is matrix and x, b are vectors, the solution of system (48) comes into existence for x if and only if rank(A)=rank(A:b); i.e., the right-hand vector must be in the vector space spanned by the columns of the coefficient matrix A. Moreover, if det(Hn(A0))≠0, according to Cramer’s rule, the solution is unique.Theorem 20 (existence).
Letf(x) be given power series (15); then (n-1/n)f and (m/n)f exist and are unique if and only if det(Hn(A0))≠0.Proof.
“⇐” By Theorem 19, if det(Hn(A0))≠0, it means that nonhomogeneous equation (48) exists as a unique solution b0,b1,…,bn-1 for vn(x). From (19), it also means that v~n(x) and W~n-1(x) exist. Hence, by the construction of (n-1/n)f, existence holds.
“⇒” Let (n-1/n)f=W~n-1(x)/v~n-1(x) exist and be unique; then it implies that (49)rankHnA0=rankHnA0:h→n, and if det(Hn(A0))=0, then rank(Hn(A0):h→n)<n; the fact that equation (46) has n-rank(Hn(A0):h→n) solutions, namely, that we can construct n-rank(Hn(A0):h→n) generating polynomials, which is contradictory to the uniqueness of (n-1/n)f, holds.
The proof of existence and uniqueness of(m/n)f is similar to the preceding process.Theorem 21.
Let det(Hn(A0))≠0; then (n-1/n)f(x)=W~n-1(x)/v~n(x), where the generating polynomial vn(x) is given by(50)vnx=detHnA0γ→δ→Txn,where γ→=An,A0,An+1,A1,…,A2n-1,An-1T, δ→T=(1,x,…,xn-1), and v~n(x) and W~n-1(x) are given by (19) and (20), respectively.Now, we can derive an algorithm to calculate(m/n)f(x) using (26), (27), and (50).Algorithm 22 (compute(m/n)f(x)=Pmn(x)/v~n(x)).
(1)
Use (14) to calculate Hn(A0) and γ→. (2)
Use (50) and (19) to compute vn(x) and v~n(x), respectively. (3)
Setq=m-n+1 and compute Wq(x) and W~q(x) by (51)Wqx=ϕm-n+1vnx-vntx-t,W~qx=xn-1Wqx-1. (4)
Compute the numerator of TPTA by(52)Pmnx=v~nx∑i=0m-nAixi+xm-n+1W~qx. (5)
Obtain(m/n)f(x)=Pmn(x)/v~n(x).


## 3.1. Definition of Tensor Pade´-Type Approximant

Letvnx=b0+b1x+⋯+bnxn(bn≠0) be a scalar polynomial of Pn of exact degree n. In this case, vn(x) is said to be quasi-monic. Define the tensor polynomial Wn-1(x) associated with vn(x) with tensor-valued coefficients, by(18)Wn-1x=ϕvnx-vntx-t.It is easily seen that Wn-1(x) is a tensor polynomial of exact degree n-1 in x. Set(19)v~nx=xnvnx-1,(20)W~n-1x=xn-1Wn-1x-1.Then, the polynomials v~n(x) and W~n-1(x) are obtained from vn(x) and Wn-1(x), respectively, by reversing the numbering of the coefficients. By the procedure given above, the following conclusion is obtained.Theorem 9.
Letv~n(0)≠0; then W~n-1(x)/v~n(x)-f(x)=O(xn).Proof.
Expandingvn(x)-vn(t)/x-t in (18) and applying ϕ yields that (21)W~n-1x=∑l=0n-1∑i=0n-l-1bl+i+1Aixl=∑l=0n-1∑i=0lbn-l+iAixl. Computing v~n(x)f(x), we get (22)v~nxfx=∑j=0nbn-jxj∑i=0∞Aixi=∑l=0∞∑i=0lbn-l+iAixl. Thus, (23)v~nxfx-W~n-1x=∑l=n∞∑i=0lbn-l+iAixl=Oxn.Definition 10.
R n - 1 , n ( x ) = W ~ n - 1 ( x ) / v ~ n ( x ) is called a tensor Pade´-type approximant with order n for the given power series (15) and is denoted by (n-1/n)f(x).Remark 11.
The polynomialvn(x), called the generating polynomial of (n-1/n)f with respect to power series f(x), can be arbitrarily chosen.Remark 12.
The tensor Pade´-type approximant (n-1/n)f possesses the degree constraint, which is caused by its construction process. The constraint implies that the method does not construct tensor Pade´-type approximant of type (m/n) when m is different from n-1.To fill this gap, we define a new tensor Pade´-type approximant by introducing a generalized linear functional.Letϕ(q):P→Rn1×n2×⋯×np be a generalized linear functional on P. Let it act on t by(24)ϕqti=Aq+i,i=0,1,2,….Similarly to what was done for Wn-1(x), we consider the polynomial Wq(x) associated with vn(x), and defined by(25)Wqx=ϕqvnx-vntx-t,q=m-n+1.Set(26)W~qx=xn-1Wqx-1,and define(27)Pmnx=v~nx∑i=0m-nAixi+xm-n+1W~qx,m≥n.Then we have the following conclusion.Theorem 13.
Letv~n(0)≠0; then Pmn(x)/v~n(x)-f(x)=O(xm+1).Proof.
Letfm-n+1 be the formal power series (28)fm-n+1x=∑j=0∞Am-n+1+jxj. Then(29)xm-n+1fm-n+1x=fx-∑i=0m-nAixi.Expanding (25) and using (26) we obtain (30)W~qx=∑q=0n-1∑i=0qbn-q+iAm-n+1+ixq. Computing the product v~n(x)fm-n+1(x), we find that (31)v~nxfm-n+1x=∑q=0∞∑i=0qbn-q+iAm-n+1+ixq. By Theorem 9, one has (32)W~qxv~nx=n-1nfm-n+1x. Then, for m≥n we deduce from (27) and (29) that (33)v~nxfx-Pmnx=v~nx∑i=0m-nAixi+xm-n+1fm-n+1x-v~x∑i=0m-nAixi+xm-n+1W~qx=xm-n+1v~nzfm-n+1-W~qx=xm-n+1∑q=n∞∑i=0qbn-q+iAm-n+1+ixq=Oxm+1.Now, we can achieve(m/n)f(x) by the above procedure, and it will be denoted by Pmn(x)/v~n(x).Definition 14.
R m , n ( x ) = P m n ( x ) / v ~ n ( x ) is called TPTA with order m+1 and is denoted by (m/n)f(x).Algorithm 15 (computeRm,n(x)=Pmn(x)/v~n(x) with vn(x) being arbitrarily chosen).
(1)
Setq=m-n+1 and chose a quasi-monic polynomial vn(x). (2)
Use (19) to compute v~n(x). (3)
ComputeWq(x) and W~q(x) by (25) and (26), respectively. (4)
Substitutev~n(x) and W~q(x) into (27) to compute Pmn(x). (5)
SetRm,n(x)=Pmn(x)/v~n(x).Example 16.
Let(34)fx=10011-102-1221+10211-212-1231x+12110-102-1320x2+⋯=A0+A1x+A2x2+⋯.Now we apply Algorithm15 to compute TPTA of type (2/2) for this example.(1)
Chosev2(x)=x2-2x+4, q=m-n+1=1.(2)
Use (19) to compute v~2(x): (35)v~2x=1-2x+4x2.(3)
By using (25) and (26) we get(36)W1x=x-122x+1x-1x-21-2xx2x-2-x+12x-13x-2x-2,and(37)W1~x=1-x2xx1-x1-2xx-102-2xx-12-x2-2x1-2x.(4)
Substitutingv~2(x) and W~1(x) into (27), we obtain(38)P22x=v~2xA0+xW1~x=3x2-x+12x2-3x2+2x3x2-x+12x2-x+1-x2-1-2x2+x6x2-2x+2-3x2+x-17x2-2x+24x2-x+22x2-x+1.(5)
SetR22(x)=P22(x)/v~2(x). It is easy to verify that (39)P22xv~2x-fx=Ox3.


## 3.2. Algorithm for Computing TPTA

Generally, the precision of TPTA is limited, since the denominator polynomials of TPTA are arbitrarily prescribed. In this subsection, in order to improve the precision of approximation, we propose an algorithm for computing the denominator polynomials and illustrate the efficiency of this algorithm in next section.First, we give the following conclusion.Theorem 17 (error formula).
(40) f x - n - 1 n f x = x n v ~ n x ϕ v n t 1 - t x .Proof.
Note thatϕ is a linear functional on P, only acting on t. From (18) and (20) we deduced that (41)W~n-1x=xn-1Wn-1x-1=xn-1ϕxvnx-1-xvnt1-tx=ϕxnvnx-1-xnvnt1-tx=ϕxnvnx-11-tx-xnϕvnt1-tx=v~nxfx-xnϕvnt)1-tx, and then this error formula holds.In terms of the error formula, it holds that(42)fx-n-1nfx=xnv~nxϕvnt1-tx=xnv~nxϕvnt+vnttx+vntt2x2+⋯=xnv~nxϕvnt+ϕvnttx+ϕvntt2x2+⋯.If we impose that vn(t) satisfies the condition ϕ(vn(t))=0, then the first term of (42) disappears, and the order of approximation becomes n+1. If, in addition, we also impose the condition ϕ(tvn(t))=0, the second term in the expansion of the error also disappears, and the order of approximation becomes n+2, and so on. We indicate that vn(x) depends on n+1 arbitrary constants; however, on the other side, a rational function is defined apart from a multiplying factor in its numerator and its denominator. It implies that (n-1/n)f(x) depends on n arbitrary constants. So let us take vn(t) such that(43)ϕvnttk=0,k=0,1,2,…,n-1.Definition 18.
v n ( t ) in (43) is called an orthogonal polynomial with respect to the linear functional ϕ and (n-1/n)f(x) in (42) is also called a TPTA for the given power series (15) when (43) is satisfied.From (43) we obtain(44)ϕvnttk=∑i=0nbiAi+k=0,k=0,1,2,…,n-1.Let bn=1 in (44); then it follows that(45)∑i=0n-1biAi+k=-Ak+n,k=0,1,2,…,n-1.Forming the scalar product of both sides of (45) with A0,A1,…,An-1, respectively, we get(46)∑i=0n-1biAi+k,Ak=-Ak+n,Ak,k=0,1,2,…,n-1.Denote(47)HnA0=A0,A0A1,A0⋯An-1,A0A1,A1A2,A1⋯An,A1⋮⋮⋱⋮An-1,An-1An,An-1⋯A2n-2,An-1,h→n=-An,A0-An+1,A1⋮-A2n-1,An-1,and call det(Hn(A0)) the Hankel determinant of f(x) with respect to the coefficients A0,A1,…,An-1.Then (46) is converted into(48)HnA0x=h→n,x=b0,b1,…,bn-1T.In the case of TPTA,vn(x) is not arbitrarily chosen any more but is determined by the preceding system. The choice of vn(x) can help to improve the accuracy of approximation, but unfortunately we have not been able to guarantee that the solution of system (48) comes into existence, so far. We only give the following basic theorem about system (48) on the basis of linear algebra.Theorem 19.
The solution of (48) exists if and only if rank(Hn(A0))=rank(Hn(A0):h→n). Moreover, the solution is unique if det(Hn(A0))≠0.Proof.
The proof of the assertion follows from the simple fact that, for a system of linear equations, described byAx=b, where A is matrix and x, b are vectors, the solution of system (48) comes into existence for x if and only if rank(A)=rank(A:b); i.e., the right-hand vector must be in the vector space spanned by the columns of the coefficient matrix A. Moreover, if det(Hn(A0))≠0, according to Cramer’s rule, the solution is unique.Theorem 20 (existence).
Letf(x) be given power series (15); then (n-1/n)f and (m/n)f exist and are unique if and only if det(Hn(A0))≠0.Proof.
“⇐” By Theorem 19, if det(Hn(A0))≠0, it means that nonhomogeneous equation (48) exists as a unique solution b0,b1,…,bn-1 for vn(x). From (19), it also means that v~n(x) and W~n-1(x) exist. Hence, by the construction of (n-1/n)f, existence holds.
“⇒” Let (n-1/n)f=W~n-1(x)/v~n-1(x) exist and be unique; then it implies that (49)rankHnA0=rankHnA0:h→n, and if det(Hn(A0))=0, then rank(Hn(A0):h→n)<n; the fact that equation (46) has n-rank(Hn(A0):h→n) solutions, namely, that we can construct n-rank(Hn(A0):h→n) generating polynomials, which is contradictory to the uniqueness of (n-1/n)f, holds.
The proof of existence and uniqueness of(m/n)f is similar to the preceding process.Theorem 21.
Let det(Hn(A0))≠0; then (n-1/n)f(x)=W~n-1(x)/v~n(x), where the generating polynomial vn(x) is given by(50)vnx=detHnA0γ→δ→Txn,where γ→=An,A0,An+1,A1,…,A2n-1,An-1T, δ→T=(1,x,…,xn-1), and v~n(x) and W~n-1(x) are given by (19) and (20), respectively.Now, we can derive an algorithm to calculate(m/n)f(x) using (26), (27), and (50).Algorithm 22 (compute(m/n)f(x)=Pmn(x)/v~n(x)).
(1)
Use (14) to calculate Hn(A0) and γ→. (2)
Use (50) and (19) to compute vn(x) and v~n(x), respectively. (3)
Setq=m-n+1 and compute Wq(x) and W~q(x) by (51)Wqx=ϕm-n+1vnx-vntx-t,W~qx=xn-1Wqx-1. (4)
Compute the numerator of TPTA by(52)Pmnx=v~nx∑i=0m-nAixi+xm-n+1W~qx. (5)
Obtain(m/n)f(x)=Pmn(x)/v~n(x).


## 4. Application for Computing the Tensor Exponential Function

The method of truncated infinite series has abroad applications in finite single crystal plasticity for computing tensor exponential function [16]. However, the accuracy and effectiveness of such algorithm are limited by round-off and choice of termination criterion. In this section, we will utilize the method of TPTA to compute tensor exponential function. We start by briefly reviewing some basic equations that model the behaviour of single crystals in the finite strain range [16].Consider a single crystal model(53)F=FeFp,where Fe and Fp denote elastic part and plastic part, respectively.For a single crystal with a total numbernsyst of slip systems, the evolution of the inelastic deformation gradient, Fp, is defined by means of the following rate form:(54)F˙pFp-1=∑α=1nsystγ˙αs0α⊗m0α,where γ˙α denotes the contribution of slip system α to the total inelastic rate of deformation. The vectors s0α and m0α denote, respectively, the slip direction and normal direction of slip system α.The above tensor differential equation can be discretized in an implicit fashion with use of the tensor exponential function. The implicit exponential approximation to the inelastic flow equation results in the following discrete form:(55)Fn+1p=exp∑α=1nsystΔγαs0α⊗m0αFnp.The above formula is analogous to the exact solution of initial value problem (1) and it is necessary to calculate(56)exp∑α=1nsystΔγαs0α⊗m0α.In [7], the author used Algorithm 23 to calculate (56).Algorithm 23 (truncated infinite series method [7] (p.749)).
(1)
Given tensorX, initialise n=0 and expX≔I. (2)
Increment countern≔n+1. (3)
Computen! and Xn. (4)
Add new term to the seriesexpX≔exp(X)+1/n!Xn. (5)
Check convergence, if(57)Xnn!<ϵtol,thenexit,elsegoto1.Example 24.
Consider a tensor exponential functionexp(Ax); the entries of A are a121=1, a221=-2, a122=2, a222=-1, and zero elsewhere.To find a tensor Pade´-type approximation of type (3/3) for the tensor exponential function, first we should expand exp(Ax) into power series by means of Definition 7. We can obtain(58)expAx=10000100+01020-20-1x+0-20-5205202x2+013601460-1460-136x3+0-40240-41240412404024x4+⋯=A0+A1x+A2x2+A3x3+A4x4+⋯. By Algorithm 22, the following can be done.(1)
Use (14) to compute H3(A0) and γ→: (59)H3A0=2-25210-14413412-61336524,γ→=-73-616-54760.(2)
Use (50) to calculate v3(x): (60)v3x=150411080+394760x+118910x2+149318x3,and
computev3~(x) by (19), so we get (61)v~3x=150411080x3+394760x2+118910x+149318.(3)
Setq=m-n+1=1 and compute W1(x), W~1(x): (62)W1x=0149318x2-422990x+1039135014939x2+5479180x+150415400-14939x2-5479180x-150415400-149318x2+422990x-1039135,and(63)W~1x=0149318-422990x+1039135x2014939+5479180x+15041540x20-14939-5479180x-15041540x20-149318+422990x-1039135x2.(4)
Substitutev~3(x) and W~1(x) into (27) to compute P33(x): (64)P33x=150411080x3+394760x2+118910x+1493181039135x3-422990x2+149318x015041540x3+5479180x2+14939x0-150411080x3+318190x2-422990x+1493180-1039135x3+422990x2-149318x.(5)
Obtain(m/n)f(x)=P33(x)/v~3(x).In Table1 we compare the number of exact figures given by the method of TPTA of type (3/3) with corresponding exact value of exp(Ax) referring to the entries of (1,2,1), (2,2,1), (1,2,2), and (2,2,2). We also compute the norm of absolute residual tensor (denoted by Res). Here, (65)Resxi=exp⁡Axi-33eAxxi2, where the operation · is defined by (13).Table 1
Numerical results of Example24 at different points by using Algorithm 22.
x  (1,2,1) ( 2,2 , 1 ) ( 1,2 , 2 ) (2,2,2) R e s 0.1 e A x 0.08200959 0.82282781 0.17717281 -0.08200959 8.34e-12 (3/3)eAx(x) 0.08200778 0.82282688 0.17717311 -0.08200778  0.2 e A x 0.13495955 0.68377119 0.31622880 -0.13495955 1.62e-9 ( 3 / 3 ) e A x ( x ) 0.13493452 0.68375764 0.31624235 -0.13493452  0.3 e A x 0.16712428 0.57369394 0.42630605 -0.16712428 3.14e-8 ( 3 / 3 ) e A x ( x ) 0.16701602 0.57363058 0.42636941 -0.16701602  0.4 e A x 0.18456291 0.48575712 0.51424287 -0.18456291 2.38e-7 ( 3 / 3 ) e A x ( x ) 0.18427224 0.48557038 0.51442961 -0.18427224From Table1, it is observed that the estimates from TPTA can reach the desired accuracy.Example 25.
Letf(x) be given by Example 24.By Algorithm22 for preceding example again, we calculate (m/m)f(1),m=1,2,3,4,5. The exact value and approximant value associated with the entries of (1,2,1), (2,2,1), (1,2,2), and (2,2,2) are listed in Tables 2 and 3, respectively, where x=1.Table 2
The exact value ofexp(A).
( 1,2 , 1 ) ( 2,2 , 1 ) ( 1,2 , 2 ) ( 2,2 , 2 ) e x p ( A ) 0.15904705 0.20883238 0.79116761 -0.15904705Table 3
Numerical approximations ofexp⁡A using Algorithm 22 for Example 25.
m ( 1,2 , 1 ) ( 2,2 , 1 ) ( 1,2 , 2 ) ( 2,2 , 2 ) R e s 1 0.50000000 0.00000000 1.00000000 - 0.50000000 3.19e-1 2 - 2.12500000 0.75000000 0.25000000 2.12500000 1.10e+1 3 0.15503865 0.20377270 0.79622729 - 0.15503865 8.33e-5 4 0.17454584 0.22682585 0.77317481 - 0.17454584 1.12e-3 5 0.17625313 0.19112365 0.80887636 - 0.17625313 1.21e-3From Table3, we can see that (3/3)eAx has the best approximation for this example. We also compute exp(A) by using Algorithm 23, and the corresponding numerical results are listed in Table 4. By comparison of Table 3 with Table 4, we find that it requires at most 6 coefficients (since m=3) of power series expansion of exp(Ax) to achieve an error of 10-5 in Algorithm 22, while requiring 11 coefficients in Algorithm 23. It is straightforward to understand that Algorithm 23 is more expensive than Algorithm 22 especially for higher order tensor exponential function. In practical applications, only few coefficients of the series may be known, so, we may get the desired results by means of TPTA. Thus the effectiveness of the proposed Algorithm 22 is verified.Table 4
Numerical results of Example25 by using Algorithm 23.
n m a x ( 1,2 , 1 ) ( 2,2 , 1 ) ( 1,2 , 2 ) ( 2,2 , 2 ) R e s 1 1 - 1 2 2 4.33 2 - 1 1.5 - 0.5 1 6.02 3 1.16666666 - 0.83333333 1.83333333 - 1.16666666 4.20 4 - 0.50000000 0.87500000 0.12500000 0.50000000 1.75 5 0.50833333 - 0.14166666 1.14166666 - 0.50833333 4.89e-1 6 0.00277777 0.36527777 0.63472222 - 0.00277777 9.77e-2 7 0.21964285 0.14821428 0.85178571 - 0.21964285 1.47e-2 8 0.13829365 0.22958829 0.77041170 - 0.13829365 1.72e-3 9 0.16541280 0.20246638 0.79753361 - 0.16541280 1.62e-4 10 0.15735119 0.21060267 0.78939732 - 0.15735119 1.22e-5 11 0.15957013 0.20838371 0.79161628 - 0.15957013 8.77e-7


## 5. Conclusion

In this paper, we presented tensor Pade´-type approximant method for computing tensor exponential function; the expression of TPTA is of the form of tensor numerator and scalar denominator. In order to have a tensor Pade´-type approximant with the higher possible precision of approximation, we proposed an algorithm for computing denominator polynomials of TPTA, and its effectiveness has been investigated in one example of tensor exponential function. The key to the TPTA to be applied to the tensor exponential function is that it can be expanded into power series with the same order tensors coefficients by means of t-product. Of course, there are several ways to multiply tensors [26–30], but the order of the resulting tensor may be changed. For example, if A is n1×n2×n3 and B is n1×m2×m3, then thecontracted product [26] of A and B is n2×n3×m2×m3. So, the choice of the multiplication of two tensors is an open question for expanding tensor exponential function, and the corresponding tensor Pade´ approximant theoretic is a subject of further research.



---

*Source: 2835175-2018-06-19.xml*