|
From: | David Bateman |
Subject: | Re: quadl crashes on evaluation |
Date: | Sat, 20 Sep 2008 15:43:03 +0200 |
User-agent: | Thunderbird 2.0.0.16 (X11/20080725) |
rmcd wrote:
Ben Abbott wrote:When you say "crash" what do you get, exactly? I get the following errors repeated many times error: evaluating binary operator `+' near line 169, column 3 ... If your function is modified f(x) = x .* sin (1 ./ (x+eps)) .* sqrt (abs (1 - x)); and an integration tolerance is added v = quadl ('f' ,0, 3,. 0001) v = 1.9818 If any event, please describe what you mean by "crash".Sorry for being unclear. By "crash" I mean that when I invoke quadl, octave exits with no error message and leaves me at the bash prompt (this is what I tried to describe at the start of my message). Presumably this should not happen. Under Ubuntu with 3.0.0 and 3.0.1 I just get the error messages, sothis seems like a cygwin-specific bug.Your suggestions to add eps to the function and a tolerance to quadl work even under cygwin, thanks very much.
It shouldn't crash as the issue here is only that the recursion limit of Octave is met. I suspect that cygwin doesn't like a 256 deep recursion as set by default with Octave and you should try something like
max_recursion_depth (128) quadl (@(x) x .* sin (1 ./ x) .* sqrt (abs (1 - x)), 0, 3)and see if it still crashes. If it doesn't, I'd suggest that the cygwin maximum recursion depth be reduced.
That being said, quadl is not really adapted to this type of singular integral. The QUADPACK based "quad" function is but is a little slow. The best for this type of singularity on the edge integrand is probably the quadgk function in the development branch of Octave as it explicitly extracts these singularities before treating the integrand and also uses vectorized calls to the integrand.
You can just take quadgk directly from the development tree and use it with Octave 3.0.x as it is not dependent on any other functions and neither is it dependent on any 3.1.x features.
Regards David
[Prev in Thread] | Current Thread | [Next in Thread] |