7.41 bug in convert[ratpoly], Maple 6 and Maple 7 (26.3.01)

7.41.1 Florian Dufey

I wanted to calculate the pade approximant to a Taylor series for given value of x.

I specified Digits:=31.

I used for this convert(series, ratpoly,i,i) (or convert(series, ratpoly,i-1,i) ) , pade(series,x,[i,i]) and convert(series,confrac,`subdiagonal`) where series is a series of order 2i+1 (2i). Finally, I substituted \(x=1\).

Independently, I calculated the corresponding approximants using Wynn’s epsilon algorithm.

For \(i <10\) all three methods lead to nearly the same results (up to some 25 digits).

For \(i>10\) , convert(ratpoly) and pade suddenly returned results which where only correct to 4 digits.

I did these calculations under MapleV rel. 2 and under Maple6.

Under Maple 6 however, the convert(confrac) gave an error message "division by zero series" in the case that the series were of order 2*i.

The help page doesn’t mention the algorithm used in the convert(confrac) routine.

7.41.2 Helmut Kahovec (28.3.01)

After I had looked at the source code of `convert/confrac` I tried the undocumented third parameter `superdiagonal` -- and it worked. Here is what I have got in Maple6:

> restart; 
> with(numapprox): 
 
> f:=proc(expr,N,x0,d) 
> local x,a1,a2,a3,a4; 
>   x:=op(indets(expr,name)); 
>   a1:=normal( 
>     pade( 
>       convert(series(expr,x=0,2*(N+1)),polynom), 
>       x, 
>       [N,N] 
>     ) 
>   ); 
>   a2:=normal( 
>     convert(series(expr,x=0,2*(N+1)),ratpoly,N,N)); 
>   a3:=normal( 
>     convert(series(expr,x=0,2*N),confrac)); 
>   a4:=normal( 
>     convert( 
>       series(expr,x=0,2*N), 
>       confrac, 
>       `superdiagonal`     # <===!!!=== 
>     ) 
>   ); 
>   if nops({a1,a2,a3,a4})=1 then 
>     map( 
>       print@(u->evalf(eval(u,x=x0),d)), 
>       [a1,expr] 
>     ) 
>   end if; 
>   NULL 
> end proc: 
 
> f(sin(x),10,1,20); 
 
                        .84147098480789650667 
 
> f(sin(x),20,2,40); 
 
              .9092974268256816953960198659117448426875 
 
> f(sin(x),30,3,60); 
 
    .141120008059867222100744802808110279846933264252265584161912
 

All four approximations a1,a2,a3,a4 were the same if the order parameters were equal to those used in f().

7.41.3 Florian Dufey (28.3.01)

As the error with the convert[ratpoly] routine doesn’t seem to appear for every input, I will give my input. The results of the convert[confrac] routine seem to converge to the correct result (4.931704...) which I also know from calculations which don’t use a series expansion.

The two methods differ first for the [10,10] approximant: confrac: 4.9330598771621... ratpoly: 4.9330573010774... the [9,9] approximants coincide to 26 digits.

file perout:

     4.906614316194628967315003238663 
    0.0000000000000000000000000000000E+00 
    0.1002750632479359264874673189523 
   -4.2825129161503777879776561918045E-04 
   -0.9552239179232503295339887788126 
    8.5761622971454992720930111738672E-03 
     17.45154716638185567939630131163 
   -0.2350619740877603889649744854291 
    -398.4547942365813625985522011213 
     7.156242135354573796385215035015 
     10188.89027936986829144576634036 
    -228.7449100349042081957499629455 
    -279146.8038168998944993427064161 
     7520.447864435336221878806554732 
     8011969.521345340499121626178387 
    -251826.2174604548840096464667442 
    -237795243.6058462537008681457654 
     8542006.312669276341086071847936 
     7238714408.232864480377075988925 
    -292532053.5581445007019944814581 
    -224760078091.0711886248616092689 
     10092326226.25219280486220429584 
     7090687290387.637279509929437663 
    -350231863205.1572850777062744521 
    -226638741608644.0867600979233863 
     12212162115859.32689138894896721 
     7323623163746757.774836679568589 
    -427512752999902.9314599007038804 
    -238860835495485678.6412699643812 
     15016037765877188.75356709241297 
     7852792004077195796.414068303233 
    -528932737787606479.5042379263840 
    -259961139015313134301.4401844604 
     18677396449998980331.82972087396 
     8658192149601240597963.646201997 
    -660947566734026317829.4553646339 
    -289917464226794648572825.0770237 
     23433697291895503144416.31715501 
     9754258473319875401003109.610984 
    -832233107565159476893114.9681807 
    -329587816875532286809985509.9060 
     29600674200295888158483002.85808 
     11179473545033574608724817973.04 
    -1054252336266437429891111919.944 
    -380528467764493572539940998621.0 
     37593898674137432962073722880.83 
    1.2993678980473297339005892294466E+31 
    -1342057905681619653173895700494. 
   -4.4497694108761200195896964312433E+32 
    4.7958316783437960933390989409083E+31 
    1.5279150172473972113542239730210E+34 
   -1.7153724800355188679239100162462E+33 
   -5.2592706432122604439688797500260E+35 
    6.1407624004915582720195948692240E+34 
    1.8144093185027082758799589502120E+37 
   -2.2000208838299084562191359760871E+36 
   -6.2727241531979608015589396915461E+38 
    7.8876292709492161887386391025367E+37 
    2.1728198209023460865624145001418E+40 
   -2.8298188245549303287218752167117E+39 
   -7.5401484974481726776559622660949E+41 
    1.0158810895294909921167517066116E+41 
    2.6210333816899045621020637730578E+43 
   -3.6490593379566447634719017835607E+42 
   -9.1254543445223108437205271472162E+44 
    1.3114636758518757265359239280792E+44 
    3.1818737742967605271178961612633E+46 
   -4.7157844399969200588659421059172E+45 
   -1.1110125924381726588637014664441E+48 
    1.6965254498816466525386748271677E+47 
    3.8844281138622831702916066729395E+49 
   -6.1060898693459081479463369491358E+48 
   -1.3597992036879097369167695853004E+51 
    2.1986245691023227400378281099625E+50 
    4.7657610803187585236358092920646E+52 
   -7.9197901648766711781523891855711E+51 
   -1.6721381079756262187612289622793E+54 
    2.8539178911023445611577952235147E+53 
    5.8731240194105115627860349479270E+55 
   -1.0287875901382004577051072130392E+55 
   -2.0649047943159536109647993557915E+57
 

maple program:

> a:=[]; 
> Digits:=31; 
> printlevel:=0; 
> for i from 1 to 82 do 
> S:=readline(`perout`): 
> a:=[op(a),op(1,sscanf(S,`%a`))]; 
> od; 
> #don't care about the error message from sscanf 
> print(a); 
> printlevel:=1; 
> l:='l'; 
> f:='f'; 
> f:=x->sum(op(l+2,a)*x^(l-1),l=1..79); 
> for i from 1 to 38 do 
> print(i-1,i,subs(x=1,(convert(series(f(x),x,2*i),confrac,'subdiagonal'))) 
              +op(1,a)); 
> print(i,i,subs(x=1,(convert(series(f(x),x,2*i+1),confrac,'subdiagonal'))) 
            +op(1,a)); 
> print(i-1,i,subs(x=1,(convert(series(f(x),x,2*i),ratpoly,i-1,i)))+op(1,a)); 
> print(i,i,subs(x=1,(convert(series(f(x),x,2*i+1),ratpoly,i,i)))+op(1,a)); 
> od;
 

7.41.4 Helmut Kahovec (2.4.01)

You certainly mean

ratpoly: \(4.9330730107748\)...

don’t you? Reading in your file perout. I can exactly reproduce your results as follows (note that I am using `superdiagonal` again):

> restart; 
> with(numapprox): 
 
> readText:=proc(fname) 
> local fd,S,line; 
>   fd:=open(fname,READ); 
>   S:=NULL; 
>   line:=readline(fname); 
>   while line<>0 do 
>     S:=S,line; 
>     line:=readline(fname) 
>   end do; 
>   close(fd); 
>   S 
> end proc: 
 
> Digits:=31: 
> L:=map(op@sscanf,[readText("perout.")],"%a"): 
> P:=sort(add(L[i+2]*x^(i-1),i=1..nops(L)-2)): 
 
> f:=proc(expr,N,x0) 
> local x,a1,a2; 
> global L; 
>   x:=op(indets(expr,name)); 
>   a1:=normal( 
      convert(series(expr,x=0,2*N+1),confrac,`superdiagonal`) 
    ); 
>   a2:=normal( 
      convert(series(expr,x=0,2*N+1),ratpoly,N,N) 
    ); 
>   evalf(eval([a1,a2],x=x0)) 
> end proc: 
 
> map(print@(u->L[1]+u),f(P,9,1)): 
 
                   4.934381129522362692062112024892 
 
> map(print@(u->L[1]+u),f(P,10,1)): 
 
                   4.933059877162190764768477498212 
                   4.933073010774838552651467870788 
 
> map(print@(u->L[1]+u),f(P,11,1)): 
 
                   4.933059894927514843532036653480 
                   4.933069963374504461817256514371 
 
> map(print@(u->L[1]+u),f(P,38,1)): 
 
                   4.931704232773919949343610941915 
                   4.933073010774838552651467870788
 

In order to get a comparable result with `convert/ratpoly` you have to use many more digits than 31. If you use 1000 digits then you get, for example:

> Digits:=1000: 
> map(print@(u->L[1]+u),f(P,38,1)): 
 
  4.93170423288275656381677632185958486834043124693630018968775046\ 
        1234529963373540963584619846083066401123053680085099607942\ 
        0328715961291010659471756681799884495422071473146251070582\ 
        7990428662066810671752383310173210403557784994307118350767\ 
        8139833827704381599458365997173259708069322507871997670569\ 
        3149071862721174701806771009248923026317495785341880626188\ 
        8522420783352995077805682170190288454479764934769003238723\ 
        2260931755373851066407563887505305457430088002269628193883\ 
        2846081319698802111935819200027072977775887603028767314734\ 
        7958387110298720716939533808609153261457530737498342094061\ 
        6855670522324379444449927633878443776415640615557205860961\ 
        0084489661519479440526372806042620798825665827135990702220\ 
        2245495924898086109889710749728154963997248182252220835264\ 
        8125577102002620506128417819889374325041611282003249139108\ 
        5390244509038217863168134174553674483686289068120239790401\ 
        7132676193116806250402677988700631926076531469912753054621\ 
        4729598064263426328050397154067263390649016770881617446965\ 
        770737846 
 
  4.93170423288275656381677632185958486834043124693630018968775046\ 
        1234529963373540963584619846083066401123053680085099607942\ 
        0328715961291010659471756681799884495422071473146251070582\ 
        7990428662066810671752383310173210403557784994307118350767\ 
        8139833827704381599458365997173259708069322507871997670569\ 
        3149071862721174701806771009248923026317495785341880626188\ 
        8522420783352995077805682170190288454479764934769003238723\ 
        2260931755373851066407563887505305457430088002269628193883\ 
        2846081319698802111935819200027072977775887603028767314734\ 
        7958387110298720716939533808609153261457530737498342094061\ 
        6855670522324379444449927633878443776415640615557205860961\ 
        0084489661519479440526372806042620798825665827135990702220\ 
        2245495924898086109889710749728154963997248182252220835264\ 
        8125577102002620506128417819889374325041611282003249139108\ 
        5390244509038217863168134174553674483686289068120239790401\ 
        7132676193116806250402677988700631926076531469912753054621\ 
        4729598064263426328050397154067263390649016770506896526738\ 
        283094823                                     ^^^^^^^^^^^^ 
        ^^^^^^^^^
 

Thus, keeping the number of digits fixed, `convert/confrac` gives much more precise results than `convert/ratpoly` if the coefficients in the polynomial to be approximated vary as much as in your case.