Padé-Hermite Approximants#
Introduction#
A given rational function \(F(z)\in k(z)\) over a field \(k\) can be written as a power series \(f(z)\) in the completion \(k((z))\) of \(k(z)\) at the place \((z)\). It thus is a good approximation of \(f(z)\) in the sense that \(F(z)\) is equal to \(f(z)\) up to infinite order.
Padé-Hermite approximants deal with the converse. Given a (formal) power series \(f(z)\in k[[z]]\) and some non- negative integers \(n_P\) and \(n_Q\), find polynomials \(P(z), Q(z)\in k[z]\) of degrees at most \(n_P\) and \(n_Q\), respectively, such that \(P(z)*f(z)-Q(z) = O(z^{n_P+n_Q+1})\) holds. In other words, \(P/Q(z)\) is an approximation for \(f(z)\) for polynomials \(P\) and \(Q\) of limited degree in \(z\). In this notation the pair \([P,Q]\) is known as the Padé-Hermite approximant for the power series tuple \((f,-1)\).
The idea of approximating one power series by two polynomials can be extended to approximating several power series at the same time as follows. Consider the vector \(\underline{f}^T:= (f_1, f_2, \ldots, f_m)^T\) in the \(m\)-dimensional vector space \(k((z))^m\) over a power series ring. Let \(\underline{n}=(n_1, n_2,\ldots, n_m)\) be an m-tuple of non-negative integers. A Padé-Hermite approximant of \(\underline{f}^T(z)\) of type \(\underline{n}\) is a non-zero vector of polynomials \(\underline{P}=(P_1,P_2,\ldots,P_m)\) in \(k[z]^m\) such that \(\underline{P}\cdot\underline{f}^T = O(z^{N})\) , \(N=n_1+n_2+\cdots+n_m+m-1\), is satisfied. A non-trivial Padé-Hermite approximant always exists. The approximant is contained in the sub-space \(V_{\underline{f},N}=\{ \underline{Q}\in k[z]^m : \underline{Q}\cdot\underline{f}^T = O(z^{N})\}\).
The implementation of Padé-Hermite approximants is based on [Derksen, 1994] and [Beckermann and Labahn, 1994]. They are implemented in Magma as sequences rather than vectors.
Ordering of Sequences#
A Padé-Hermite approximant to some sequence of power series does not have to be unique. They can be ordered according to their maximum degree and their type of sequence. A sequence \(\underline{P}=[P_1, P_2,..., P_m]\) is of degree \(i\) if the maximum of the weak degrees of \(P_1, P_2, \ldots, P_m\) is \(i\). An extension of the definition maximum degree is when a distortion of non-negative integers \(\underline{d}=[d_1, d_2, \ldots, d_m]\) on the degrees is allowed. In this case the maximum degree is defined as the maximum of \(\deg{P_i}-d_i\). The type of \(\underline{P}\) identifies the last \(P_i\) in \(\underline{P}\) whose weak degree equals the maximum degree of \(\underline{P}\). Again a distortion on the degrees is allowed.
- MaximumDegree(f): SeqEnum -> RngIntElt#
Distortion: SeqEnum Default: []
MaximumDegree returns the degree of a sequence of polynomials or power series, defined as the maximum of the degrees of \(f[i]-d[i]\), where \(d\) is the distortion. The value -infinity is returned in the case that \(f\) is weakly equal to the zero-sequence.
- Example: Degree Of Sequence (ex-ea74e8)#
> S<u> := PowerSeriesRing(Rationals()); > f := [u+u^2, 2+u^2+u^3,0]; > MaximumDegree(f); 3 > MaximumDegree(f:Distortion:=[]); 3 > MaximumDegree(f:Distortion:=[0,2,1]); 2 > MaximumDegree([S|0,0]); -Infinity > MaximumDegree([O(u)]); -Infinity
- TypeOfSequence(f): SeqEnum -> RngIntElt, RngIntElt#
Distortion: SeqEnum Default: []
Returns the highest index \(i\) for those \(f[i]\) whose (distorted) degree is weakly equal to the maximum of the degrees of all entries. The second integer returned is the maximum degree of the sequence.
- Example: Type Of Sequence (ex-8d75ca)#
> S<u> := PowerSeriesRing(Rationals()); > f := [u+u^2, 2+u^2+u^3,0]; > TypeOfSequence(f); 2 3 > TypeOfSequence(f:Distortion:=[]); 2 3 > TypeOfSequence(f:Distortion:=[0,2,1]); 1 2 > TypeOfSequence([S|0,0]); 2 -Infinity
The Padé-Hermite approximant of \(\underline{f}\) can be seen as an element of the free module \(k[z]^m\) of rank \(m\), or as an element of the sub-space \(V_{\underline{f},N}\).
A free sub-model \(V\subset k[z]^m\) is generated by \(m\) polynomial vectors \(\underline{Q}_i\), \(1=1,2,\ldots,m\), such that \(\underline{Q}_1(z), \underline{Q}_2(z), \ldots, \underline{Q}_m(z)\) form a minimal vector sequence of \(V\). Such a sequence is defined as a sequence \(S\) of \(m\) vectors in \(V\), such that \(S[i]\) is a non-trivial polynomial vector in \(V\) of minimal degree of type \(i\), for \(i=1,2,\ldots,m\). A minimal vector sequence is not unique.
Two variations on a minimal vector sequence are implemented. The first allows a distortion as an attribute. The sequence is the based on the distorted maximal degree. The second attribute sets the positive power \(p\) of \(z\) in each of the \(\underline{Q}_i\) as follows. Instead of considering \(\underline{Q}_i(z)\) for each \(i\), one considers \(\underline{Q}_i(z^p)\).
- MinimalVectorSequence(f, n): SeqEnum, RngIntElt -> SeqEnum#
Distortion: SeqEnum Default: [] Power : RngIntElt Default: 1
A minimal sequence of vectors \(\underline{Q}_1, \underline{Q}_2, \ldots, \underline{Q}_m\) with respect to the sequence \(f\) of length \(m\) whose entries are polynomials or power series. The order of \(\underline{Q}_i\cdot f\) is at least \(m\).
- Example: Minimal Vector Sequence (ex-67b928)#
> S<u> := PowerSeriesRing(Rationals()); > f := [u+u^2, 2+u^2+u^3]; > seq := MinimalVectorSequence(f, 2); > seq; [ (u 0), ( -1 1/2*u) ] > sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq]; > sums; [ u^2 + u^3, -u^2 + 1/2*u^3 + 1/2*u^4 ] > seq := MinimalVectorSequence(f,3); > #seq eq 2, seq[1], seq[2]; true (u^2 0) (-1 + u 1/2*u)
- Example: The Next Example (ex-d1c2b6)#
> L<x> := PolynomialRing(Rationals()); > f := [1+x, 3-x^2, 5+x+x^3-x^5]; > seq := MinimalVectorSequence(f, 10); > seq[1]; (-2*x^2 + 6 -2*x - 2 0) > sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq]; > sums; [ 0, 0, x^10 ]
- Example: Another Example (ex-328015)#
> S<u> := PowerSeriesRing(Rationals()); > f := [2*u^4,2+u^3+u^6]; > seq := MinimalVectorSequence(f, 10); > seq; [ (1/2*u^6 0), (-1/2 - 1/4*u^3 1/2*u^4) ] > sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq]; > sums; [ u^10, 1/2*u^10 ]
- Example: One More (ex-1d53c0)#
> S<u> := PowerSeriesRing(Rationals()); > f := [1+u-7*u^2, 6-3*u+1/2*u^2-u^3, 5-u+u^2]; > seq := MinimalVectorSequence(f, 5); > [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq]; [ 0, u^5, 0 ] > seq := MinimalVectorSequence(f, 5:Distortion:=[2,0,1]); > [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq]; [ -7/2*u^5, u^5, 0 ] > p:=2; > seq := MinimalVectorSequence(f, 5:Distortion:=[2,0,1], Power:=p); > sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq]; > sums; > mp:= map<S->S| x :-> (IsWeaklyZero(x) select 0 > else &+([Coefficient(x,i)*(S.1)^(p*i) : i in Exponents(x)])) > + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v)) > else S!0 where v := AbsolutePrecision(x))>; > sums := [&+([mp(Q[i])*f[i]: i in [1..#f]]) : Q in seq]; > sums; [ u^6 + u^7 - 7*u^8, u^5 - 23/3*u^6, -5/3*u^5 + 35/3*u^6 ]
Approximants#
The Padé-Hermite approximant of type \(\underline{d}=[d_1,d_2,\ldots,d_m]\) with respect to the tuple \(\underline{f}^T\in k((z))^m\) is an element of the space \(V_{\underline{f},N}=\{ \underline{Q}\in k[z]^m : \underline{Q}\cdot\underline{f}^T = O(z^{N})\}\) for \(N\) equal to \(d_1+d_2+\cdots+d_m+m-1\). This space is generated by the vectors in the minimal vector sequence with respect to \(\underline{f}\) with distortion \(\underline{d}\). The routine PadeHermiteApproximant returns one that is smallest with respect to the degree on sequences. The input sequence \(\underline{f}\) must be a sequence of polynomial ring elements, or be a power series sequence. While the Padé Hermite approximants theoretically are polynomials, Magma returns them as elements of the same ring the entries of \(\underline{f}\) are contained in.
- PadeHermiteApproximant(f, d): SeqEnum, SeqEnum -> ModTupRngElt, SeqEnum, RngIntElt#
Power: RngIntElt Default: 1
Returns a Padé-Hermite form \(\underline{P}\) of \(f\) with distortion \(d\), smallest with respect to the degree on sequences, and the corresponding minimal vector sequence. The third argument returned is the order in the order term of \(\underline{P}\cdot f\).
- Example: Pade Hermite Approximants (ex-05a6f0)#
This example can be found on page 813 in [Beckermann and Labahn, 1994].
> S<u> := PowerSeriesRing(Rationals()); > f := [1,u,u/(1-u^4)+u^10+O(u^16),u/(1+u^4)+u^12+O(u^16)]; > pade, padebasis, ord := PadeHermiteApproximant(f,[2,2,2,2]); > pade, ord; ( u -1 0 0) 11 > BaseRing(Parent(pade)) eq S; true > MinimalVectorSequence(f,10); [ ( u -1 0 0), ( 0 u^4 -1/2 1/2), ( 0 1 - u^4 -1/2 + u^4 -1/2), ( 0 -u 1/2*u 1/2*u) ] > > p := 2; > seq := MinimalVectorSequence(f,10: Distortion :=[2,2,2,2],Power := p); > seq; [ (u^5 0 0 0), ( 0 u^2 -1/2 1/2), ( 0 1 - u^2 -1/2 + u^2 -1/2), ( 0 -u 1/2*u 1/2*u) ] > mp:= map<S->S| x :-> (IsWeaklyZero(x) select 0 > else &+([Coefficient(x,i)*(S.1)^(p*i) : i in Exponents(x)])) > + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v)) > else S!0 where v := AbsolutePrecision(x))>; > sums := [&+([mp(Q[i])*f[i]: i in [1..#f]]) : Q in seq]; > [Valuation(v) : v in sums]; [ 10, 10, 10, 11 ] > sums; [ u^10, -1/2*u^10 + 1/2*u^12 - u^13 + O(u^16), -1/2*u^10 - 1/2*u^12 + u^13 + u^14 + O(u^16), u^11 + 1/2*u^12 + 1/2*u^14 + O(u^18) ]
- Example: Last Example (ex-e98d49)#
This example covers the example on page 815 in [Beckermann and Labahn, 1994].
> S<u> := PowerSeriesRing(Rationals()); > f := [1,u, -1-u^4-2*u^8+u^10+u^11-u^12+O(u^16),-u-u^5-u^9-u^14-u^15+O(u^16)]; > dist:=[2,2,3,3]; > seq := MinimalVectorSequence(f,13:Distortion:=dist); > pade, padebasis, ord := PadeHermiteApproximant(f,dist); > pade, ord; (-u 1 0 0) 13 > padebasis; [ (-u 1 0 0), (1/2*u -1/2*u^4 1/2*u + 1/2*u^3 + 1/2*u^4 -1/2*u^2 - 1/2*u^3 - u^4), (-1/2 1/2*u^3 -1/2 - 1/2*u^2 - 1/2*u^3 - u^4 1/2*u + 1/2*u^2 + 2*u^3), ( -u 0 0 -1 + u^4) ] > padebasis eq seq; true > [[Valuation(w[i]): i in [1..Degree(w)]] : w in seq]; [ [ 1, 0, Infinity, Infinity ], [ 1, 4, 1, 2 ], [ 0, 3, 0, 1 ], [ 1, Infinity, Infinity, 0 ] ] > [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq]; [ [ -1, -2, -Infinity, -Infinity ], [ -1, 2, 1, 1 ], [ -2, 1, 1, 0 ], [ -1, -Infinity, -Infinity, 1 ] ] > p:=2; > seq := MinimalVectorSequence(f,12:Distortion:=dist,Power:=p); > seq; [ ( u - u^3 0 u - 2*u^3 0), (-1 - u + u^2 -u^3 -1 - u + 2*u^2 + u^3 -u^3), (u + u^2 - u^3 0 u + u^2 - 2*u^3 - u^4 0), ( 0 1 0 1 - u^2) ] %> [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq]; %[ % [ 1, -Infinity, 0, -Infinity ], % [ 0, 1, 0, 0 ], % [ 1, -Infinity, 1, -Infinity ], % [ -Infinity, -2, -Infinity, -1 ] %] > seq[1]-seq[3]; ( -u^2 0 -u^2 + u^4 0) > mp:= map<S->S| x :-> (IsWeaklyZero(x) select 0 > else &+([Coefficient(x,i)*(S.1)^(p*i) : i in Exponents(x)])) > + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v)) > else S!0 where v := AbsolutePrecision(x))>; > [Valuation(&+([mp(Q[i])*f[i]: i in [1..#f]])) : Q in seq]; [ 12, 12, 13, 13 ]
- Example (ex-ea09f6)#
This example considers the example on page 816 in [Beckermann and Labahn, 1994].
> S<u> := PowerSeriesRing(Rationals()); > f := [1,u,-1-u^4-2*u^8+u^10+O(u^12),-u-u^5-u^9+u^10+O(u^12)]; > dist := [2,2,3,3]; > seq := MinimalVectorSequence(f,12: Distortion := dist); > [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq]; [ [ -1, -2, -Infinity, -Infinity ], [ -2, 1, 0, 0 ], [ -Infinity, -Infinity, 1, 0 ], [ -1, -Infinity, 0, 1 ] ] > [Valuation(&+([(Q[i])*f[i]: i in [1..#f]]) ) : Q in seq]; [ Infinity, 12, 12, 12 ] > [MaximumDegree([ &+([(Q[i])*f[i]: i in [1..#f]]) ]) : Q in seq]; [ -Infinity, -Infinity, 14, -Infinity ] > PadeHermiteApproximant(f,[2,2,3,3]); % error % % redefine with more precision %> g := [1,u,-1-u^4-2*u^8+u^10+O(u^13),-u-u^5-u^9+u^10+O(u^13)]; %> seq := MinimalVectorSequence(g,12:Distortion:=[2,2,3,3]); %> PadeHermiteApproximant(g,[2,2,3,3]); %(-u 1 0 0) %[ % (-u 1 0 0), % (1/4*u -1/4*u^4 1/4*u + 1/4*u^3 - 1/4*u^4 -1/4*u^2 + 1/4*u^3 - 1/2*u^4), % (-1/4 1/4*u^3 -1/4 - 1/4*u^2 + 1/4*u^3 - u^4 1/4*u - 1/4*u^2 + 3/2*u^3), % (1/2 - u -1/2*u^3 1/2 - 1/2*u^2 + 1/2*u^3 -1 + 1/2*u - 1/2*u^2 - u^3 + u^4) %] %13 > p := 2; > seq := MinimalVectorSequence(f,12:Distortion:=[2,2,3,3],Power:=p); > seq; [ ( 1 - u^2 -1 1 - 2*u^2 -1 + u^2), ( 0 -u^4 0 -u^4), ( -u -1 -u + u^3 -1 + u^2), ( 0 u 0 u - u^3) ] > [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq]; [ [ 0, -2, -1, -1 ], [ -Infinity, 2, -Infinity, 1 ], [ -1, -2, 0, -1 ], [ -Infinity, -1, -Infinity, 0 ] ] > mp:= map<S->S| x :-> (IsWeaklyZero(x) select 0 > else &+([Coefficient(x,i)*(S.1)^(p*i) : i in Exponents(x)])) > + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v)) > else S!0 where v := AbsolutePrecision(x))>; > > [Valuation(&+([mp(Q[i])*f[i]: i in [1..#f]])) : Q in seq]; [ 12, 13, 12, 12 ]
A variant of the Padé-Hermite approximant is when the exponent in the order term is set rather that the type of the sequence. It is also possible to let \(\underline{f}\) be a sequence such that its entries themselves are vectors of polynomials or power series.
- PadeHermiteApproximant(f, m): SeqEnum, RngIntElt -> ModTupRngElt, SeqEnum#
Power: RngIntElt Default: 1
Returns a Padé-Hermite form of minimal degree in the corresponding minimal vector sequence, such that its inproduct with \(f\) has order at least \(m\). The second argument returned is the corresponding minimal vector sequence.
- Example: Pade Hermite Approximants Vectors (ex-b5d830)#
This example can be found on page 813 in [Beckermann and Labahn, 1994].
> S<u> := PowerSeriesRing(Rationals()); > f := [Vector([1]), Vector([u])]; > pade, seq := PadeHermiteApproximant(f,3); Calculating the Pade'-Hermite approximant for the sequence [ 1, u ] with order term 3 and power 1 . > pade; ( u -1) > seq; [ ( u -1), ( 0 u^2) ] > mat := Matrix([Eltseq(v): v in f]); > pade*mat; (0) > PadeHermiteApproximant([1,u],5); ( u -1) [ ( u -1), ( 0 u^4) ] > PadeHermiteApproximant(f,3:Power:=2); > g:= [Vector([1,0,0]), Vector([0,1,0]), Vector([1+u,2+u^2,u^3])]; > pade := PadeHermiteApproximant(g,5); Calculating the Pade'-Hermite approximant for the sequence [ 1, u, 1 + 2*u + u^3 + u^7 + u^11 ] with order term 15 and power 3 . > pade; (-u^3 - u^4 -2*u^3 u^3) > pade*Matrix([Eltseq(v): v in g]); ( 0 u^5 u^6)
- Example (ex-15b388)#
This example considers Padé-Hermite approximants for some series that have non-trivial power series expansions.
> S<u> := PowerSeriesRing(Rationals()); > f := [Sin(u), Cos(u), Exp(u)]; > [Valuation(f[i]) : i in [1..#f]], [Degree(f[i]) : i in [1..#f]]; [ 1, 0, 0 ] [ 19, 20, 20 ] > [AbsolutePrecision(f[i]) : i in [1..#f]]; [ 21, 22, 21 ] > dist := [3,2,5]; > pade, seq, ord := PadeHermiteApproximant(f,dist); > 1/420*pade; (-1275 - 255*u + 45*u^2 + 5*u^3 120 + 495*u + 75*u^2 -120 + 900*u - 600*u^2 + 160*u^3 - 20*u^4 + u^5) > ord eq &+(dist)+#f-1, ord; true 12 > [Degree(pade[i]) : i in [1..Degree(pade)]]; [ 3, 2, 5 ] % The next one is in the Maple package nummapprox[hermite-pade]. > g:= [Cos(2*u)*(u+1)+3,Cos(u)^2+u*Cos(u)+1,Cos(2*u)+1,Cos(u)]; > pade, basis := PadeHermiteApproximant(g,20); > 131/75880*pade; ( 2 -4 + 2*u -3*u 4*u - 2*u^2) % This is one from Maple's OrderBasis. % It returns a different minimal vector sequence. > h := [ 1+u^2-u^7+u^12, Sin(u), Exp(u) ]; > dist:=[3,1,2]; > seq := MinimalVectorSequence(h,8:Distortion := dist); > sums := [&+([Q[i]*h[i]: i in [1..#f]]) : Q in seq]; > [Valuation(s) : s in sums]; [ 8, 8, 8 ] > [[MaximumDegree([w[i]]): i in [1..Degree(w)] ] : w in seq]; [ [ 4, 1, 2 ], [ 4, 2, 2 ], [ 3, 1, 2 ] ] > [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq]; [ [ 1, 0, 0 ], [ 1, 1, 0 ], [ 0, 0, 0 ] ]