Edwin Woollett | 1 Jan 09:00 2002
Picon

integrate bessel_k(2,%i*y) surprise

I am just playing around with bessel_k and don't
understand the large imaginary part of the integral
of bessel_k(2, %i*y) over the interval
1 <=  y <= 100, based on the table of discrete 
values displayed below:
(f(x) is just fchop(expand(float(x))) )
--------------------------------
(%i1) load(nint);
(%o1) "c:/work2/nint.mac"

(%i25) for y:1 step 10 thru 100 do 
          print("  ",y,"  ",f(bessel_k(2,%i*y)))$

   1    0.18048997206696*%i-2.592886175491198 
   11    0.21841533174753*%i+0.3119789483129 
   21    -0.031858737423089*%i-0.27222015896945 
   31    0.20481495858656-0.093918476739506*%i 
   41    0.16377581778393*%i-0.10738879820159 
   51    0.0024786070162412-0.17554486214757*%i 
   61    0.13641169853138*%i+0.084590710242742 
   71    -0.064015208750726*%i-0.13429138505077 
   81    0.13815292462619-0.017659551649835*%i 
   91    0.084630268914206*%i-0.10051430087896 

(%i26) f(integrate(bessel_k(2,%i*y),y,1,100));

(%o26) 1.736293756153568-2.1367387806763253E+42*%i
----------------------------------------------------------------
Where does the 10^42 come from??

(Continue reading)

Edwin Woollett | 1 Jan 09:26 2002
Picon

Re: integrate bessel_k(2,%i*y) surprise

On Oct. 26, 2011, I wrote:
------------------
>(%i26) f(integrate(bessel_k(2,%i*y),y,1,100));
>
>(%o26) 1.736293756153568-2.1367387806763253E+42*%i
>----------------------------------------------------------------
>Where does the 10^42 come from??
----------------------------
I should have included some quad_qag output:

(%i27) quad_qag(imagpart(bessel_k(2.0,%i*y)),y,1,100,3,limit=700);

(%o27) [1.62942476406395,1.1838977094676129E-9,217,0]

Ted
Edwin Woollett | 1 Jan 10:35 2002
Picon

Re: numeric approximation

On Feb. 7, 2012, Rupert Swarbrick wrote:
-----------------------------
>Notice that there are several numerical integration routines and that
>(as far as I know) clever numerical analyst types haven't found a way to
>automatically determine which one is most appropriate, so it's not clear
>how one would implement the "gimme_a_number" flag. Also, numerical
>routines also return a guess at their error, which you're throwing away
>in your code above. How should that be treated by the "gimme_a_number"
>flag?
-----------------------------------
A "wrapper" for the quadpack routines is something I have been
diligently working on, as part of a software package for
an updated version of Ch. 8, Numerical Integration,
(Maxima by Example).

The one-dimensional quadrature code is basically done,
and I am now trying to integrate two dimensional
code into the package.

This code first tries integrate to get an exact
symbolic answer (if possible). Due to some bugs
in the special functions area, there are automatic
filters to bypass integrate and go directly to
quadpack (this also includes multiple valued
functions case).

The code then chooses between the quadpack routines
using a variety of methods and option flags.

If quadpack returns a serious error code, the function
(Continue reading)

Edwin Woollett | 1 Jan 11:29 2002
Picon

Re: numeric approximation

On Feb. 7, 2012, I wrote:
------------------------------------
(%i1) load(nint);
(%o1) "c:/work2/nint.mac"

(%i2) g : 1/(sin(x)^2 + log(x))$

(%i3) f(t) := nint(g,x,1,t)$

(%i4) f(5);

(%o4) 2.729767710374242

(%i5) goutL;

(%o5) [[qag,2.729767710374242,3.3769724748722686E-12,93,0]]
--------------------------------
then one can make a table, for example:

(%i6) ntable (func, x0, dx, xf):= 
block([nL, fL, nfL, jj, ii],
   nL : makelist (zz, zz, x0, xf, dx),
   fL : map ('func, nL),
   nfL : makelist ( ["  ", nL[jj], "  ", fL[jj] ], jj, length(nL)),
   for ii thru length (nfL) do apply ('print, nfL[ii]))$

(%i7) ntable ( 'f, 2, 2, 10)$

   2    0.80395247146645 
   4    2.307860373006648 
(Continue reading)

Edwin Woollett | 1 Jan 13:31 2002
Picon

Re: display patch, was: get name of variable as string in function?

On Feb. 6, 2012, Robert Dodier wrote:
------------------------
....No need to invent a work-around, it's easier to fix the bug.
Here is a patch.....
-----------------------------
Thanks for the patch.

Ted
Edwin Woollett | 1 Jan 15:07 2002
Picon

Re: numeric approximation

On Feb. 7, 2012, I wrote
-----------------------------------------
>(%i1) load(nint);
>(%o1) "c:/work2/nint.mac"

>(%i2) g : 1/(sin(x)^2 + log(x))$

>(%i3) f(t) := nint(g,x,1,t)$

>(%i4) f(5);

>(%o4) 2.729767710374242

>(%i5) goutL;

>(%o5) [[qag,2.729767710374242,3.3769724748722686E-12,93,0]]
>--------------------------------
>then one can make a table, for example:

>(%i6) ntable (func, x0, dx, xf):= 
>block([nL, fL, nfL, jj, ii],
>   nL : makelist (zz, zz, x0, xf, dx),
>   fL : map ('func, nL),
>   nfL : makelist ( ["  ", nL[jj], "  ", fL[jj] ], jj, length(nL)),
>   for ii thru length (nfL) do apply ('print, nfL[ii]))$

>(%i7) ntable ( 'f, 2, 2, 10)$

>   2    0.80395247146645 
>   4    2.307860373006648 
(Continue reading)


Gmane