Finding sequence of algebraic approximants

Mathematica program

(* Finding sequence of algebraic approximants *)

BeginPackage["algappr`"]

algappr0::usage =
"   algappr0[coeff, ndegree] finds sequence of algebraic approximants.\n
   Input:\n
coeff - array of coefficients of the series;\n
ndegree - degree of the approximant.\n
   Returnes {coeff[[1]], C},
where coeff[[1]] - zero-order coefficient of the series,\n
C - coefficients of recurrence relations for finding of algebraic approximants.\n
   Options:\n
ndigits (def. 1024) - extra precision,\n
timeprint (def. 10) - how often to print diagnostic messages."

ndigits::usage =
"ndigits is an option to algappr0 specifying the number of digits during calculations."

timeprint::usage =
"timeprint is an option to algappr0 specifying period of time in sec.\n
   of printing messages indicating progress of calculations."

algappr1::usage =
"   algappr1[a0, cf, x, nm] finds all roots of the polynomial equation for the algebraic approximants.\n
   Input:\n
a0,cf - f(0) and array of coef. of recurrent relat. calculated using algappr1;\n
x - argument;\n
nm - effective number of coef. of the series to be incorporated in alg. appr.\n
   Returnes {y1, y2, ..., y_ndegree}."

algappr::usage =
"   algappr[coeff, ndegree, x, nm] finds algebraic approximant as all roots of the polynomial equation.\n
   Useful for low-order calculations of alg. appr.
   Input:\n
coeff - array of coefficients of the series;\n
ndegree - degree of the approximant.\n
x - argument;\n
nm - effective number of coef. of the series to be incorporated in alg. appr.\n
   Returnes {y1, y2, ..., y_ndegree}."

algapprs::usage =
"   algapprs[coeff, ndegr, ncoeff] finds all singularities of the algebraic approximant.\n
   Input:\n
coeff - array of coefficients of the series;\n
ndegree - degree of the approximant.\n
nm - effective number of coef. of the series to be incorporated in alg. appr.\n
   Returnes {y1, y2, ..., } where |y1| <= |y2| <= ..."

feenb2::usage =
"   feenb2[coeff, lambda1, lambda2] performs generalized Feenberg transformation.\n
   Input:\n
coeff - array of coefficients of the series;\n
lambda1, lambda2 - parameters of transformation.\n
   Returnes coefficients of transformed series."

Begin["`private`"]

Options[algappr0] = {ndigits->1024, timeprint->10};

algappr0[e_, ndegr_, opts___]:=
(
{ndig, timep} = {ndigits, timeprint}/.{opts}/.Options[algappr0];
km = ndegr + 1;
settime := (time = SessionTime[]);
chktime := (SessionTime[] - time >= timep);
settime;
nm = Length[e];
a0=SetPrecision[e[[1]],ndig];
Do[e1[n-1]=SetPrecision[e[[n]],ndig],{n,2,nm}];
ek[1,1]=1;
Do[ek[1,n]=0, {n,2,nm}];
Do[ek[2,n]=e1[n-1], {n,2,nm}];
Do[
    If[chktime, Print["Calculating: Power ", k - 1, " of ", ndegr]; 
      settime];
   Do[ek[k,n]=Sum[e1[n-l]*ek[k-1,l],{l,k-1,n-1}],
      {n,k,nm}],{k,3,km}];
Do[Do[c[k,n]=ek[k,n],{n,k,nm}],{k,1,km}];
Do[
    If[chktime, Print["Calculating: Order ", ka - 1, " of ", nm - 1]; 
      settime];
   nf=ka-km;
   nf2=nf+2;
   Do[n=nm-n5+nf2;
      c[1,n]=c[1,n-1],{n5,nf2,nm}];
   c[1,nf+1]=1;
   Do[n0=nf+k;
      a=c[1,n0]/c[k,n0];
      c[k,nf+1]=-a;
      cf[k,nf+1]=c[k,nf+1];
      Do[c[1,n]=c[1,n]-a*c[k,n], {n,n0,nm}], {k,2,km}];
   Do[a=c[1,n];
      Do[c[k-1,n]=c[k,n], {k,2,km}];
   c[km,n]=a, {n,nf2,nm}], {ka,km,nm}];
{a0, Table[cf[k, n], {k, 2, km}, {n, nm-km+1}]}
)

algappr1[c0_, cfrec_, x_, ncoeff_, opts___]:=
(
{ndig, timep} = {ndigits, timeprint}/.{opts}/.Options[algappr0];
{ndegree, nm1} = Dimensions[cfrec];
nm0 = nm1 + ndegree;
If[ncoeff>nm0, Print["Error in ALGAPPR1: Number of coefficients (", ncoeff, ") should not exceed ", nm0]; Return[]];
km = ndegree + 1;
x1 = SetPrecision[x, ndig];
Do[Do[r[k, n] = (-c0)^(n - k)*Binomial[n - 1, k - 1], {n, km}], {k, km}];
Do[na = n - km + 1;
    Do[a = x1*r[k, 1];
      Do[a = a + cfrec[[m, na]]*r[k, m + 1], {m, ndegree}];
      Do[r[k, m] = r[k, m + 1], {m, ndegree}];
      r[k, km] = Expand[a], {k, km}], {n, km, ncoeff}];
If[x === xnosolve, Null,
pol0 = Sum[r[n +1, km]*t^n, {n, 0, ndegree}];
solv = If[NumberQ[x1] && Precision[x1]=!=Infinity, NSolve[pol0 == 0, t], Solve[pol0 == 0, t]];
t /. # & /@ solv]
)

algappr[coeff_, x_, ndegr_, ncoeff_, opts___]:=
(
nm0 = Length[coeff];
If[ncoeff>nm0, Print["Error in ALGAPPR: Number of coefficients (", ncoeff, ") should not exceed ", nm0]; Return[]];
If[ndegr>ncoeff, Print["Error in ALGAPPR: Degree of approximant (", ndegr, ") should not exceed ", ncoeff]; Return[]];
algappr1[Sequence@@algappr0[Take[coeff,ncoeff], ndegr, opts], x, ncoeff, opts]
)

algapprs[coeff_, ndegr_, ncoeff_, opts___]:=
(
algappr[coeff, xnosolve, ndegr, ncoeff, opts];
pol = Sum[r[n +1, km]*t^n, {n, 0, ndegr}];
a = r[ndegr +1, km];
diskr = If[ndegr == 1, a,
   Cancel[Resultant[pol, D[pol, t], t]/a]
 ]//Expand;
(* Solving polynom. eq. *)
s = NSolve[diskr == 0, xnosolve];
Sort[xnosolve /. # & /@ s, Abs[#1]<Abs[#2]&]
)

feenb2[coeff_, lambda1_, lambda2_]:=
(
nm = Length[coeff];
Prepend[
Table[n = n1 - 1; Sum[k = n - i - j; 
      Binomial[i + j - 1, i - 1]Binomial[i + k - 1, i - 1]*
      If[j==0,1,lambda1^j] If[k==0,1,lambda2^k] (1 - lambda1)^i (1 - lambda2)^i coeff[[i + 1]],
      {i, 0, n}, {j, 0, n - i}], {n1, 2, nm}], coeff[[1]] ]
)

End[]

EndPackage[]

printprogr["algappr"];

Molecule - icon for Allen-dataBlankExamples of MP seriesBlankSource code of Mathematica programBlankMathematica programsBlankWork in UMass DartmouthWork in UMassDBlankWaste iconUnpublished reports

Designed by A. Sergeev.