Difficulty in detecting true discontinuities

[Tim Mitchell]

As I recently mentioned on he...@chebfun.org, I want to find discontinuities of certain functions.  Unfortunately, distinguishing between a true discontinuity and a breakpoint can be difficult, since current versions of Chebfun no longer guarantee exact agreement of function values at breakpoints (where the function in continuous).  So we have to look at how large these gaps are.  If the gaps are small, then it is usually easy to sort out the actual discontinuities (which have large gaps) from all the other breakpoints.  However, for the functions I’ve been using, even breakpoints where the function is still continuous can sometimes have relatively large gaps, i.e., 1e-7.  This makes it hard to reliably detect which breakpoints are the actual discontinuities in an automatic way.

I’ve created a simple demo, with just two discontinuities (at x=1 and x=3), which demonstrates this problem:

>> breakpointsDemo(false)  % doesn’t use merge 

>> breakpointsDemo(true)   % Uses merge but merge isn’t helpful

 

For now, to get around this problem, I am now using a slightly modified Chebfun which uses a global variable to stores all the edge values found by findJump routine inside @fun/detectEdge.  This works rather well, despite being a rather ugly kludge.  The odd thing is that findJump repeatedly finds the same two discontinuities, over and over again.  So in the end, I get a vector of 10 edges, with x=1 and x=3, each with multiplicity five (modulo rounding error).

Moving forward, it would great if edges detected by findJump could be added as a new field of the resulting chebfun in a future release.  That way this info would be immediately available to the user.  Lastly, if possible, it would also be nice if findJump routine didn’t have find the same discontinuities more than once, as this would be more efficient (particularly since findJump evaluates the function at value at a time in its bisection iteration so there’s no vectorization here).  Perhaps this too is possible?

Many thanks,

Tim

[Nick Hale]

However, for the functions I’ve been using, even breakpoints where the function is still continuous can sometimes have relatively large gaps, i.e., 1e-7. 

I can't solve your problem, but I can at least explain why this occurs:

Your function has square root branch singularities in addition to the jumps, which invokes 

the adaptive splitting. The error control in the adaptive splitting is not 'L_\infty' in the sense 

that it allows the newly created funs to have a larger error if they live on smaller intervals. 

The norm this then defines I don't know (if it even is a norm) but it at least (usually) allows 

the construction to terminate. 

In particular, if you look at the 6th fun in the constructed chebfun (the one with length 8, 

closest to the singularity / jump) ones finds

chebcoeffs(f.funs{6})

     9.999099322018328e-01

     5.994850992857182e-05

     1.193195459526655e-05

     5.072226977573609e-06

     2.786094459815282e-06

     1.746665523485859e-06

     1.187486566134302e-06

     8.516762607987371e-07

Suggesting this fun only has a local accuracy of around 1e-7 (since it lives on a domain

of length 1e-8, it's 'global' accuracy around 1e-15). Hence, when you evaluate this fun at

it's right and/or left limit, you can expect an error of around 1e-7.

Storing the jumps located in fundjump would be a big design change, and one I don't 

think would find much traction. In particular, those jumps are properties of the original 

function, not the constructed chebfun. No other properties of the construction function

are stored (for example, one could keep hold of the function handle used to construct

in case that was useful later, but we don't).

Hope this helps

Nick 

[Tim Mitchell]

Dear Nick,

Thanks for the additional clarifications and sorry for my *extremely* tardy response.  

In the end, I just added a global variable inside the findJump routine so that I could more reliably obtain the locations of any jump discontinuities.  It's a bit of a kludge but it suffices for my purposes now.

I understand your point that incorporating this info into the output would be a big design change and so is unlikely to happen.  Perhaps it would be easier and be more relevant to users if chebfun.m could optionally return a second output argument containing various metadata?  For example, I personally would also find it helpful to have a status code that tells me if chebfun.m converged or not, as I often use chebfun as a subroutine.  The current warnings are difficult to catch in a program, i.e., I think one has to resort to undocumented features in MATLAB.  In any case, I understand this too may still be a big design change and there might not be sufficient interest.

Thanks and best regards,

Tim

[Nick Hale]

Hi Tim - I will match your tardy response with a prompt one!

The ishappy function will actually check for convergence (which we call 'happiness'). 

Here's an example:

>> ishappy(chebfun(@(x) sqrt(1-x)))

ans =

  logical

   0

>> ishappy(chebfun(@(x) sin(x)))

ans =

  logical

   1

Hopefully that will solve some of your issues.

I agree that at times some other outputs/metadata/intermediate computations might be useful if they were user-accessible, but as you say, that would be quite a big design change.

Best regards

Nick

[Tim Mitchell]

Hi Nick, 

Yes very helpful!  When I was perusing the source code last year, I did notice ishappy() and its name definitely made me smile when I saw it.  But for some reason, it never occurred to me to call it myself!  I'll definitely take advantage of this, thanks! 

Best,

Tim