(* 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"]; |

Examples of MP seriesMathematica programsWork in UMassDUnpublished reports |

Designed by A. Sergeev.