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 speciﬁed 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(
>       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 diﬀer ﬁrst for the [10,10] approximant: confrac: 4.9330598771621... ratpoly: 4.9330573010774... the [9,9] approximants coincide to 26 digits.

ﬁle 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 ﬁle perout. I can exactly reproduce your results as follows (note that I am using superdiagonal again):

> restart;
> with(numapprox):

> local fd,S,line;
>   S:=NULL;
>   while line<>0 do
>     S:=S,line;
>   end do;
>   close(fd);
>   S
> end proc:

> Digits:=31:

> 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 ﬁxed, convert/confrac gives much more precise results than convert/ratpoly if the coeﬃcients in the polynomial to be approximated vary as much as in your case.