Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Tables of special functions

• To: mathgroup <mathgroup at yoda.physics.unc.edu>
• Subject: Tables of special functions
• From: S.Arridge at cs.ucl.ac.uk
• Date: Mon, 26 Oct 92 18:46:48 +0000

I've got a problem with the dereaded Mathematica precision racheting. Using
version 2.0 on Sparc 2

I need to produce large (~2000 terms) tables of special functions, specifically 
Bessels and Legendres, at quite high precision. Easy enough :
Version 1 )
tab = Table[N[Bessel[i - 1/2,x],P], {i, M+1}]
tab = Table[N[LegendreP[i,x],P], {i, M+1}]

However, it takes a long time. When I write my own functions in C, at normal
precision, I use the reccurrence relations for these functions, so that the
whole array is produced in the same time as the highest order member. I
would like to do the same in Mathematica, to take advantage of high precision,
since presumably, inside the function generation the recurrence relations are
being used.

Naively this should work:
Version 2 )
SetLegendreTables/: SetLegendreTables[ tab_, M_, x_] :=
       {
       	tab = Table[0, {i, M+1}];
       	tab[[1]]  = 1;
	tab[[2]]  = x;
	Do[
	   {tab[[k+1]] =  ( (2 k -1) x tab[[k-1]] - (k - 1)tab[[k-2]] )/(k) } ,
	   {k,2,M}
	]
       }  
If x is symbolic, (e.g. Pi/10) the calculation is symbolic and takes forever.
If x is say N[Pi/10, 40] then the classic precision ratcheting occurs.
Reading the Mathematica journal I tried this :
Version 3)
SetLegendreTables/: SetLegendreTables[ tab_, M_, x_, P_] :=
       {
       	tab = Table[0, {i, M+1}];
       	tab[[1]]  = 1;
	tab[[2]]  = N[x, P];
	Do[
	   {tab[[k+1]] = ((2 k -1) x tab[[k]] - (k -1)tab[[k-1]])/k ;
	    tab[[k+1]] = SetPrecision[tab[[k+1]],P] },
	   {k,2,M}
	]
       } 

I.e. keep x symbolic, reset precision to P on each iteration. This seems to
simply zero pad the entries in tab[[k+1]] to the requested precision, so the
accuracy is still disintegrating. The only one that works is this :
Version 4)
SetLegendreTables/: SetLegendreTables[ tab_, M_, x_, P_] :=
       {
       	tab = Table[0, {i, M+1}];
       	tab[[1]]  = 1;
	tab[[2]]  = N[x, P+M];
	Do[
	   {tab[[k+1]] = N[((2 k -1) x tab[[k]] - (k -1)tab[[k-1]]),P+M-k]/k ;
	    tab[[k-1]] = N[tab[[k-1]],P] },
	   {k,2,M}
	]
       }  

I.e. start with such high precision that the ratcheting only brings you down
to the requested precision at the highest order member. For Bessels the
equivalent is (using downward recursion for stability) :
SetBesselTables/: SetBesselTables[ tab_, M_, x_, P_] :=
       {
       	tab = Table[0, {i, M+1}];
       	tab[[M+1]] = BesselI[M + 1/2, N[x, P+M] ];
	tab[[M ]]  = BesselI[M - 1/2, N[x, P+M] ];
	Do[
	   {k = M-j; 
	    tab[[k]] = N[(tab[[k+2]] + (2 k + 1)tab[[k+1]]/x), P+k];
	    tab[[k+2]] = N[tab[[k+2]],P]} ,
	   {j,1,M-1}
	]
       }

The trouble with these solutions is that carrying the massive precision around
makes the solution longer than version 1).

So my question is - assuming that N[LegendreP[2000,x],P] or 
N[BesselI[2000 + 1/2, x],P ] has already computed the terms 1-1999 at at least 
precision P, how can I get hold of them ?

Simon Arridge