Thiele 型張量連分式插值及其在張量指數(shù)計算中的應(yīng)用

蔣祥龍，

(上海大學(xué)理學(xué)院， 上海200444)

張量常微分方程的初值問題[1]可以表述為

式中，A和Y0是給定的常張量.該張量常微分方程的解為

式中， exp((t-t0)A)就是著名的張量指數(shù)函數(shù).對于給定的常張量A， 張量指數(shù)一般表示為如下的級數(shù)表達(dá)式:

在文獻(xiàn)[1]中， 式(3)被用來近似計算或者逼近張量指數(shù)， 即

式中， 截斷的最高項nmax滿足

本研究提出了一種張量廣義逆Thiele 型連分式插值方法， 用來近似計算式(2)中的張量指數(shù)函數(shù).該方法可以看作矩陣廣義逆Thiele 型連分式插值方法[2–5]在張量上的推廣， 且具有如下特點: ①在計算過程中，不必計算張量的乘積; ②該方法用連分式表示， 是一種遞推計算方法.

1 有關(guān)張量的幾個基本概念

張量是一個多維數(shù)組.一階張量是一個向量， 一個二階張量是一個矩陣， 并且是一階張量.3 個或3 個以上的張量稱為高階張量.所謂張量矩陣化， 就是把張量用矩陣的形式展開[6].例如， 考慮一個三階張量A= (ai1，i2，i3)∈R2×3×3.如果指定張量A的第3 個指標(biāo)， 可以得到3 個2×3 矩陣， 由此該張量可以表示為

定義1 和定義2 為3 階張量t乘積(t-product)的定義.該定義自然地以遞歸的方式推廣到高階張量.

定義1[7]設(shè)A ∈Rl×m×n.張量A的塊循環(huán)矩陣為

定義張量的展開(unfold(·))[7]為

再定義張量的折疊(fold(·))[7]， 二者互為相反逆運算， 從而滿足fold(unfold(A))=A.

定義2[7]設(shè)A和B分別是l×p×n和p×m×n張量.它們的3 階張量乘積A*B是如下的l×m×n張量:

定義3 和定義4 為張量的內(nèi)積和張量的范數(shù)， 可以看作是矩陣的直接推廣.

定義3[8]設(shè)A=a(i1，i2，···，ip)∈Cn1×n2×···×np，B=b(i1，i2，···，ip)∈Cn1×n2×···×np.張量A和B的內(nèi)積定義為

定義4[8]設(shè)A=a(i1，i2，···，ip)∈Cn1×n2×···×np， 它的范數(shù)定義為

根據(jù)式(5)和(6)， 不難證明

式中，A*表示張量A的共軛轉(zhuǎn)置張量.

2 Thiele 型廣義逆張量連分式插值

本研究首先給出基于式(7)的一種張量廣義逆的定義， 然后在此基礎(chǔ)上構(gòu)造了Thiele 型張量連分式插值公式.這種張量廣義逆可以看作相應(yīng)的矩陣廣義逆[2-5]在張量上的推廣， 并已被成功應(yīng)用于張量Pad′e 型逼近方法[9]、張量廣義逆Pad′e 逼近的連分式方法[10]和APCILON方法[11].

定義5設(shè)張量A ∈Cn1×n2×···×np且A/=0， 張量范數(shù)‖A ‖2=A·A*由式(7)給出.張量A基于內(nèi)積的廣義逆定義為

或

式中， 廣義逆表示張量A的倒數(shù).

按照張量廣義逆的計算公式(8)或(9)， 定義Thiele 型張量截斷連分式如下: 設(shè)插值點的集合為Φ={xi， i=1，2，··· ，n， xi ∈R}.根據(jù)式(8)或(9)， 遞推地構(gòu)造第n個Thiele 型截斷連分式為

式中，

根據(jù)表達(dá)式(10)和(11)， 給出如下的遞推算法1.

算法1計算截斷插值連分式(10).

步驟1 給定初始值(xi，Ai)，i=0，1，··· ，n.

步驟2 Forl=0，1，··· ，n.

步驟3 計算Bl(x0，x1，··· ，xl).

步驟4 如果Bl(x0，x1，··· ，xl)=0， 轉(zhuǎn)步驟1.結(jié)束.

步驟5 結(jié)束.

步驟6 通過式(10)計算Rn(x).

定理1～定理3 分別為Thiele 型廣義逆張量連分式插值公式(10)的存在性、整除性和特征性定理.考慮到它們的證明方法與相應(yīng)的矩陣方法基本相同， 此處只給出有關(guān)結(jié)論.

定理1(存在性)[3]設(shè)Bl(x0，x1，··· ，xl)∈Cn1×n2×···×np[x]，0 ≤1 ≤n， 存在且不等于0(except forB0(x0)).同時，

滿足

則式(10)中(x)存在， 并成立插值條件(xi)=Ri，xi ∈Φ.

定理2(整除性)[3]設(shè)Bl(x0，x1，··· ，xl)∈Cn1×n2×···×np[x]，0 ≤l≤n，xi ∈Φ.在式(10)中，從后逐次向前利用張量廣義逆計算(xi)，則一個張量多項式P(x)∈Cn1×n2×···×np[x]和一個實多項式q(x)存在， 并且滿足

式中，q(x)|‖P(x)‖2表示分母實多項式能夠整除分子張量多項式的范數(shù)的平方.

定理3(特征性)[3]根據(jù)定理2， 設(shè)式(10)中截斷連分式為(x) =P(x)/q(x)， 則成立: ①如果n是偶數(shù)，(x)具有階數(shù)為[n/n]型; ②如果n是奇數(shù)，(x)具有階數(shù)為[n/n-1]型， 其中[l/k]型表示分子張量多項式為l階， 分母實多項式為k階.

例1 設(shè)插值點x0=-1，x1= 0，x2= 1， 相應(yīng)的待插值張量值A(chǔ)i ∈C2×2×3，i= 0，1，2，則

根據(jù)連分式插值公式(10)和(11)， 利用張量廣義逆(8)或(9)分別得到

不難驗證插值條件R2(xi)=A(xi)，i=0，1，2.

3 在張量指數(shù)函數(shù)計算中的應(yīng)用

給定張量A和t0= 0， 本研究利用張量連分式插值公式(10)， 計算張量常微分方程的解(2)， 或者在如下解

中， 近似計算張量指數(shù)函數(shù)exp(At)的值.本研究計算exp(At)值的方法分為如下幾個步驟:①利用張量的t乘積(見定義2)來計算張量A的冪; ②導(dǎo)出張量指數(shù)函數(shù)exp(At)的Tailor 展開式; ③選擇合適的插值點， 例如， 考慮在原點附近選擇張量A的特征值; ④計算相應(yīng)的插值函數(shù)的值; ⑤根據(jù)算法1 得到連分式插值函數(shù)R(T)， 并利用其計算張量指數(shù)函數(shù)的值; ⑥考慮計算值是否滿足誤差要求， 具體步驟如算法2 所示.

算法2計算張量指數(shù)函數(shù)的值.

步驟1 利用張量的t乘積計算An， n=0，1，···， 從而得到exp(At)的Tailor 展開式.

步驟2 在t0=0 臨近選擇合適的插值點ti， i=0，1，··· ，n.

步驟3 計算待插值的值exp(Ati)， i=0，1，··· ，n.

步驟4 由算法1 得到連分式插值公式Rn(t).

步驟5 給定誤差限ε， 計算誤差Res(ti)=‖exp(Ati)-Rn(ti)‖.

步驟6 如果Res(t)≥ε， 轉(zhuǎn)步驟2， 否則， 輸出Rn(t).結(jié)束.

例2 在張量指數(shù)函數(shù)exp(At)中， 給定計算張量指數(shù)函數(shù)exp(At)在臨近t0=0 的值.

步驟1 利用張量的t乘積計算An， n=0，1，···， 得到

步驟2 計算張量A的特征值[12]為0，1/3，1/2，5/6， 選擇插值點為t0= 0， t1= 1/3， t2=1/2.

步驟3 計算待插值的張量指數(shù)函數(shù)exp(Ati)， i=0，1，2 的值為

步驟4 根據(jù)步驟2 和步驟3， 由算法1 得到

于是得到連分式插值公式

式中:

步驟5 設(shè)誤差限ε= 1×10-3， 計算誤差Res(ti) =‖exp(Ati)-R2(ti)‖.這里， 取張量指數(shù)函數(shù)Tailor 展開式前15 項的和作為exp(At)的精確值.同時， 取連分式插值函數(shù)R2(t)計算的值作為exp(At)的近似值， 結(jié)果如表1 所示， 其中exp(At)的精確值與近似值分別用Evalue 和A-value 表示.

表1 R2(t)近似計算exp(At)的數(shù)值結(jié)果Table 1 Comparisons of numerical results of R2(t) for exp(At)

從表1 可以看出， 當(dāng)插值點在0～0.6 之間時， 連分式插值函數(shù)R2(t)計算張量指數(shù)函數(shù)的值的時候， 逼近效果比較好， 特別當(dāng)插值點t在0.3～0.5 之間時， 數(shù)值結(jié)果更加精確.