non-deterministic find_roots function (was strange quad_qags behavior?)

Sorry, I posted the big float version which requires access to bffind_root() 
which I didn't post.  Here is attached the float version which just uses 
Maxima's find_root function. It probably has bugs.

-----Original Message----- 
From: Richard Hennessy
Sent: Saturday, December 29, 2012 2:55 PM
To: Barton Willis ; maxima at math.utexas.edu
Subject: non-deterministic find_roots function (was strange 
quad_qags behavior?)

Hi,

I decided to try my hand at a root finder that can find most or in some
cases (with a little luck) all roots in a specified interval of a function
f(x).  I have written it so that it partitions the interval into n parts
where n is specified by the user.  Then in each interval, it recursively
partitions the interval into even smaller intervals, which are chosen
partially to divide in half but with a random jiggle . The maximum number of
recursive calls is also specified by the user in a global variable, so it
will not recursively go on forever and generate stack overflow errors. It
works pretty well but I want to improve efficiency.  When a function has a
relatively high number of roots in a subinterval like sin(1/x) between [0,
.1] then the program needs to be able to detect that the function is likely
to have a high number of roots there and try more samples inside that
subinterval.  How can you determine that a subinterval like [.0, .01] has
more roots than [.01, .02] for any function, not just sin(1/x) and any
interval [a,b] not just [0,.01]?  I have been using the ratio of the
standard deviation of the list to the mean of the list below.

[f(a+randnum1),f(b+randnum2),f((a+randnum1+b+randnum2)/2].  Randnum1 and
randnum2 are small random numbers (50% possibility of being negative). Is
there a better way to anticipate where the function is likely to have more
roots?

Here is a sample session.

(%i36) display2d:false;
(out36) false
(%i37) fpprec:16;
(out37) 16
(%i38) find_roots(exp(x)-5*x^2+3,x,-1,10,10);
(out38) [-8.290321491121718b-1,1.094253308924241b0,4.659015373364428b0]
(%i39) fpprec:32;
(out39) 32
(%i40) find_roots(exp(x)-5*x^2+3,x,-1,10,10);
(out40)
[-8.2903214911217184394046498484114b-1,1.0942533089242409104622697911055b0,4.6590153733644280966688664369581b0]
(%i41) find_roots(sin(1/x),x,-1,1,10);
(out41)
[-3.1830988618379067153776752674503b-1,-1.5915494309189533576888376337251b-1,-1.0610329539459689051258917558168b-1,
         -7.9577471545947667884441881686257b-2,-6.3661977236758134307553505349006b-2,-5.3051647697298445256294587790838b-2,
         3.1830988618379067153776752674503b-1]
(%i42) find_roots(sin(1/x),x,-1,1,10);
(out42)
[-3.1830988618379067153776752674503b-1,1.5915494309189533576888376337251b-1,3.1830988618379067153776752674503b-1]
(%i43) find_roots(sin(1/x),x,-1,1,10);
(out43)
[-3.1830988618379067153776752674503b-1,3.1830988618379067153776752674503b-1]
(%i44) find_roots(sin(1/x),x,-1,1,10);
(out44)
[-3.1830988618379067153776752674503b-1,3.1830988618379067153776752674503b-1]
(%i45) find_roots(sin(1/x),x,-1,1,10);
(out45)
[-3.1830988618379067153776752674503b-1,-1.5915494309189533576888376337251b-1,-1.0610329539459689051258917558168b-1,
         1.8724110951987686561045148632061b-2,6.3661977236758134307553505349006b-2,7.9577471545947667884441881686257b-2,
         1.5915494309189533576888376337251b-1,3.1830988618379067153776752674503b-1]

Notice that the number of roots found is non-deterministic for sin(1/x).

Rich

-----Original Message----- 
From: Richard Hennessy
Sent: Sunday, December 23, 2012 6:42 PM
To: Barton Willis ; maxima at math.utexas.edu
Subject: Re: [Maxima] strange quad_qags behavior?


"%solve misses many (most) important cases; for example %solve knows nothing
about bessel functions. Consider:

  (%i24) quad_qags (1/bessel_j(0,x), x, 2, 10, 'epsrel=1d-10);
  (%o24) [19.18306454272454,49.39600889427612,7077,5]

Oops... bessel_j(0,x) has three zeros in [2, 10] ; sometimes a visual
analysis wins:

  (%i25) wxplot2d(bessel_j(0,x),[x,2,10]);

--Barton"

I agree that visual analysis may be the best way.  But if expr is univariate
then you can use find_root(1/expr,x,a,b).  You just have to figure out a and
b, which may be easier that finding the analytical solution, if there is
one.  My TI-89 calculator can find a and b in most cases.  It has a nsolve()
command which does not even require a and b. I assume it samples the real
number line according to some algorithm and finds a and b on its own.  I
does accept help in cases by allowing the user to enter a guess for where
the root is.  Maybe someone knows what algorithm TI-89 uses.  I think you
could for integrate(1/bessel_j(7,x),x,c,d) just try (c,(c+d)/2) and
((c+d)/2, d), and then rinse, repeat.  That would always work.

Rich






_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima 
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: findroots.txt
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20121229/13d96d92/attachment-0001.txt>;