Well, first of all, I think the usual way to calculate bessel functions for
large argument is NOT taylor series at all, but asymptotic forms. I don't
know about negative argument though.
Next, a correct taylor series for bessel_j(1/6,x) should look approximately
like
0.960311*x^(1/6)-0.205781*x^(13/6)+0.011872*x^(25.6) - ... according to
Mathematica.
So if you compute x^(1/6) first, you are talking about 16 multiplies and
adds.
So you are clearly not computing that if you take 33 seconds.
When I try computing a taylor series taylor(bessel_j(1/6,x),x,0,16), I
get a mess of
stuff with bessel_j( n/6, 0) for various n. But these are, I think, all
zeros, and
the taylor series is exactly zero, if you evaluate it numerically.
So I know why you got the time to be zero.
It could be that defining Bessel_J of fractional order by taylor series is
not a good idea.
The usual stuff people use is, I think recurrences of some sort.
> -----Original Message-----
> From: maxima-bounces at math.utexas.edu
> [mailto:maxima-bounces at math.utexas.edu] On Behalf Of Richard Hennessy
> Sent: Thursday, June 12, 2008 8:01 AM
> To: fateman at EECS.Berkeley.EDU
> Cc: maxima at math.utexas.edu
> Subject: Re: [Maxima] ansi gcl in rpm
>
> Sorry there are a few extra steps I left out
>
> (%i9) float(mybess(1/6,-10));
> Evaluation took 0.2600 seconds (0.2600 elapsed)
> (eq 9) -0.22332422870625*(-1)^(1/6)
>
> then I have to use rectform(%) to get the answer in standard
> form for complex numbers.
>
> (%i14) rectform(-0.22332422870625*(-1)^(1/6));
> Evaluation took 0.1500 seconds (0.1500 elapsed)
> (eq 14) -0.11166211435312*%i-0.11166211435312*sqrt(3)
>
> (%i15) float(-0.11166211435312*%i-0.11166211435312*sqrt(3)), numer;
> Evaluation took 0.0700 seconds (0.0700 elapsed)
> (eq 15) -0.11166211435312*%i-0.19340445534017
>
> That is in agreement with the bessel run at 100 and 1600
> terms so you don't need all 160 of those terms except for
> large z values, but 16 terms using taylor gives:
>
> (%i16) taylor(bessel_j(1/6,x),x,0,16);
> horner(%,x);
> optimize(%);
> mybess(x):=''%;
> compile(mybess);
> Evaluation took 33.5100 seconds (33.5100 elapsed)
>
> (%i22) rectform(float(mybess(-10)));
> Evaluation took 0.0000 seconds (0.0000 elapsed)
> (eq 22) 0
>
> Now I am really confused. I don't know which is right now.
> Maybe the series is not the definition of bessel_j or maybe I
> made a typo somewhere. Also I don't know why I got 44
> seconds for the run time before, now it is almost 0 seconds.
>
> One thing though that has bothered me is the fact that there
> are six possible answers to (-1)^(1/6) and that might make
> the series inaccurate and a poor way to define the bessel function.
>
> Rich
>
>
>
> ------------Original Message------------
> From: "Richard Fateman" <fateman at cs.berkeley.edu>
> To: "'Richard Hennessy'" <rvh2007 at comcast.net>
> Cc: maxima at math.utexas.edu
> Date: Thu, Jun-12-2008 9:07 AM
> Subject: RE: [Maxima] ansi gcl in rpm
>
> I get quite different results. I tried the expr to only 16,
> rather than
> 160.
> Then mybess(1/6,-10) gives
> -1082845135081018553168355738*(-1)^(1/6)*10^(1/6)/(11391066481
247301946\
> 59276527*2^(1/6)*gamma(1/6))
>
> which is not a number at all, unless you do something more.
>
> the second, taylor() approach gives an expression with a
> large number of
> constants such as bessel_j(7/6,0) which are, in my maxima, evaluated
> numerically
> to 0, so the 2nd mybess() evaluates to 0.
>
> Now I don't know what you really expect, but if the taylor()
> expression has
> lots of constant bessel_j
> expressions,then evaluating the whole mess numerically before
> you compile,
> should result in a polynomial of degree
> 16 in x, wit h floating-point coefficients, and no summation
> of 160 terms
> should be faster.
>
> As for evaluating accurately, the expression expr to 16th
> term, evaluated at
> 1/6, and then converted to floats, looks like this:
> 2.5254888134545135E-116 z
> 32 30
> (1.174882790201127E+79 z - 1.2156120602614335E+82 z
> 28 26
> + 1.1062069748379136E+85 z - 8.7759086670474326E+87 z
> 24 22
> + 6.0085721340384927E+90 z - 3.5090061262784811E+93 z
> 20 18
> + 1.7240916767114858E+96 z - 7.0113061519600392E+98 z
> 16 14
> + 2.3137310301468166E+101 z - 6.0465504254503521E+103 z
> 12 10
> + 1.2133411187070328E+106 z - 1.7957448556864159E+108 z
> 8 6
> + 1.8556030175426283E+110 z - 1.2370686783617525E+112 z
> 4 2
> + 4.700860977774658E+113 z - 8.1481590281427346E+114 z
> + 3.80247421313328E+115)
>
> So there are terms like 10^79*(-10)^320 and 10^114*(-10)^20
> being added
> together.
> The second term is quite negligible as a float. I do not
> know why it should
> be advantageous to compute it at all, unless you are
> interested in an exact
> rational number of huge size.
>
> What am I missing here??
>
> RJF
>
>
> > -----Original Message-----
> > From: maxima-bounces at math.utexas.edu
> > [mailto:maxima-bounces at math.utexas.edu] On Behalf Of
> Richard Hennessy
> > Sent: Thursday, June 12, 2008 4:58 AM
> > To: fateman at EECS.Berkeley.EDU
> > Cc: maxima at math.utexas.edu
> > Subject: Re: [Maxima] ansi gcl in rpm
> >
> > This is how I am running the function.
> > expr:sum((-1)^k*2^(-v-2*k)*z^(v+2*k)/(k!*gamma(v+k+1)),k,0,160)$
> > mybess(v,z):=''%$
> >
> > mybess(1/6,-10) // takes about .270 seconds
> >
> > Now if I do this
> >
> > taylor(bessel_j(1/6,x),x,0,16); //last parameter of 20 runs
> > out of memory so I lowered it to 16
> > horner(%,x);
> > optimize(%);
> > mybess(x):=''%;
> > compile(mybess);
> >
> > now
> > mybess(-10) takes about 44 seconds.
> > It seems like a factor of about 200 difference and the first
> > way is much more accurate.
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > ------------Original Message------------
> > From: "Richard Fateman" <fateman at cs.berkeley.edu>
> > To: "'Richard Hennessy'" <rvh2007 at comcast.net>
> > Cc: maxima at math.utexas.edu
> > Date: Wed, Jun-11-2008 9:38 PM
> > Subject: RE: [Maxima] ansi gcl in rpm
> >
> > there are other things you should do, like converting to
> > floats, declaring x
> > to be a float, perhaps other issues.
> >
> > what exactly are you doing to generate the series?
> >
> >
> > > -----Original Message-----
> > > From: Richard Hennessy [mailto:rvh2007 at comcast.net]
> > > Sent: Wednesday, June 11, 2008 2:42 PM
> > > To: fateman at EECS.Berkeley.EDU
> > > Cc: 'Maxima List'
> > > Subject: RE: [Maxima] ansi gcl in rpm
> > >
> > > This is worse than just using the series directly. It may be
> > > fast for the sin function but not the bessels.
> > >
> > >
> > >
> > >
> > > ------------Original Message------------
> > > From: "Richard Fateman" <fateman at cs.berkeley.edu>
> > > To: "'Richard Hennessy'" <rvh2007 at comcast.net>
> > > Cc: "'Maxima List'" <maxima at math.utexas.edu>
> > > Date: Tue, Jun-10-2008 6:40 PM
> > > Subject: RE: [Maxima] ansi gcl in rpm
> > >
> > > try taylor(sin(x),x,0,20);
> > > horner(%,x);
> > > optimize(%);
> > > f(x):=''%;
> > > compile(f);
> > >
> > >
> > > > -----Original Message-----
> > > > From: maxima-bounces at math.utexas.edu
> > > > [mailto:maxima-bounces at math.utexas.edu] On Behalf Of
> > > Richard Hennessy
> > > > Sent: Tuesday, June 10, 2008 3:30 PM
> > > > To: fateman at EECS.Berkeley.EDU; 'Michel Talon'
> > > > Cc: 'Maxima List'
> > > > Subject: Re: [Maxima] ansi gcl in rpm
> > > >
> > > > I meant taylor takes too long not the power series.
> > > >
> > > >
> > >
> > >
> > >
> >
> >
> > _______________________________________________
> > Maxima mailing list
> > Maxima at math.utexas.edu
> > http://www.math.utexas.edu/mailman/listinfo/maxima
> >
>
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>