# How to test if a function is continuous

**URL:** <https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737>\
**Category:** Numerics\
**Tags:** question\
**Created:** [April 9, 2024, 3:23pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737 "2024-04-09T15:23:37Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 9, 2024, 3:23pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/1 "2024-04-09T15:23:38Z")

</div>

I am writing some tests for a package (not yet public) where I want to test that a function `f(x::Real)::Real` is continuous at some `x`.

Is there a numerical way of doing this?

A solution I thought of is finding the left- and right limits with Richardson extrapolation, as in

```julia
using Richardson

function diff_at_x(f, x, h; kwargs...)
    f_left, f_left_err = extrapolate(f, -h; x0 = x, kwargs...)
    f_right, f_right_err = extrapolate(f, h; x0 = x, kwargs...)
    # roughly, < 1 is OK as limits are conservative
    abs(f_left - f_right) / (abs(f_left_err) + abs(f_right_err) + eps(x))
end

julia> diff_at_x(abs, 0.0, 0.1)
0.0

julia> diff_at_x(x -> x > 0 ? sin(x) : sin(x) - 1, 0.0, 0.1)
2.7536530891901533e12

```

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [April 9, 2024, 4:22pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/2 "2024-04-09T16:22:27Z")

</div>

Can you autodiff through it? A similar idea to what you’re doing now would be to compare `ForwardDiff.derivative(f, x+eps(x))` and `ForwardDiff.derivative(f, x-eps(x))` or something like that.

---

<div class="post-metadata">

**Author:** ![NikoBiele](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikobiele/32/208184_2.png) [@NikoBiele](https://discourse.julialang.org/u/NikoBiele)\
**Post date:** [April 9, 2024, 5:09pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/3 "2024-04-09T17:09:19Z")

</div>

I’m not an expert, but how about using FFT? Shouldn’t discontinuous functions give a sharp peak in a FFT? Due to a sudden change in frequency?

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [April 9, 2024, 5:46pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/4 "2024-04-09T17:46:52Z")

</div>

You need a whole lot of frequencies to approximate a discontinuity, so it wouldn’t be a peak. Explaining it as a sudden change in frequency, whatever that means, seems questionable at best. Which frequency would suddenly change to which if you have `f(x) = sign(x)`?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 9, 2024, 5:53pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/5 "2024-04-09T17:53:23Z")

</div>

> [@Tamas\_Papp](#):
>
> A solution I thought of is finding the left- and right limits with Richardson extrapolation

If the function is sufficiently well behaved, you probably don’t need to extrapolate? i.e. just compare `f(nextfloat(x))` and `f(prevfloat(x))`?

> [@cgeoga](#):
>
> Can you autodiff through it?

Most AD tools silently take one-sided derivatives at discontinuities and kinks, e.g.

```julia
julia> using ForwardDiff

julia> ForwardDiff.derivative(x -> x > 0 ? 1 : 0, 0)
0

julia> ForwardDiff.derivative(abs, 0)
1

```

> [@cgeoga](#):
>
> A similar idea to what you’re doing now would be to compare `ForwardDiff.derivative(f, x+eps(x))` and `ForwardDiff.derivative(f, x-eps(x))` or something like that.

That’s checking whether the (one-sided) derivative is continuous, not the function itself?

---

<div class="post-metadata">

**Author:** ![NikoBiele](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikobiele/32/208184_2.png) [@NikoBiele](https://discourse.julialang.org/u/NikoBiele)\
**Post date:** [April 9, 2024, 6:06pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/6 "2024-04-09T18:06:18Z")

</div>

Fair points, that’s why I said I’m not an expert on FFT. 🙂 You’re right that a sudden change of slope is not a change of frequency. I guess it is more like an impulse, which contains all frequencies. Thanks for the clarification.

I guess I was thinking along the lines of those artifacts (peaks) that occur when you Fourier series a square wave function. But I guess derivative is a better method. Like the derivate of the Heaviside step function can be understood as the dirac delta function at the step.

EDIT: I was just trying to brainstorm but I guess it was not that helpful…

---

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [April 9, 2024, 9:04pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/7 "2024-04-09T21:04:04Z")

</div>

Also, in order to calculate an FFT (DFT) you’d need to sample the function’s outputs at a number of points in the domain around the particular value. A discontinuity anywhere in that sampled region, not just at the desired point, would probably be indistinguishable.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 9, 2024, 9:23pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/8 "2024-04-09T21:23:40Z")

</div>

> [@mike.ingold](#):
>
> Also, in order to calculate an FFT (DFT) you’d need to sample the function’s outputs at a number of points in the domain around the particular value. A discontinuity anywhere in that sampled region, not just at the desired point, would probably be indistinguishable.

Note that you wouldn’t use a DFT or Fourier series directly, because that implies periodic boundaries (so a generic non-periodic function will effectively always have “discontinuities” at the endpoints from the perspective of the DFT). If you have a function on an arbitrary finite interval x \in [a,b] and you want to detect discontinuities and other singularities automatically (e.g. at unknown positions), you could expand in Chebyshev polynomials (e.g. via ApproxFun.jl or FastChebInterp.jl … these use FFTs internally because they are related to a DFT by a change of variables), and look at the convergence rate of the Chebyshev coefficients c\_n — if it is converging as only \sim 1/n asymptotically, that indicates a discontinuity somewhere in the interior (a,b) of the domain. (More generally, slower-than-exponential asymptotic convergence indicates some non-analyticity.)

If you want to locate unknown discontinuities/singularities in an interval [a,b] (if any) there are various recursive search algorithms. e.g. the algorithm used by [Chebfun](http://www.chebfun.org/) with the `splitting on` option is described in the article [“Piecewise smooth chebfuns”](https://www.cs.ox.ac.uk/files/717/NA-08-07.pdf) ([2010](https://doi.org/10.1093/imanum/drp008))

![image](https://global.discourse-cdn.com/julialang/original/3X/4/7/477693f5afeab63eb63d34959300b9792119a45f.png)

Their [`findjump` code](https://github.com/chebfun/chebfun/blob/db207bc9f48278ca4def15bf90591bfa44d0801d/%40fun/detectEdge.m#L177-L258) (free/open-source under [3-clause BSD](https://github.com/chebfun/chebfun/blob/db207bc9f48278ca4def15bf90591bfa44d0801d/LICENSE.txt) is pretty short and wouldn’t be too hard to re-implement in Julia.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [April 10, 2024, 1:30am UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/9 "2024-04-10T01:30:23Z")

</div>

Oh geez, you’re right. I read that too quickly and thought Tamas was asking for a test of differentiability (where the point was to compare the two one-sided derivatives), not continuity. Sorry guys.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 10, 2024, 5:18am UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/10 "2024-04-10T05:18:28Z")

</div>

> [@stevengj](#):
>
> i.e. just compare `f(nextfloat(x))` and `f(prevfloat(x))`?

I thought about this but could not figure out a reasonable cutoff. Eg consider

```julia
function f(x; A)
    z = x - 1
    z > 0 ? A * z : 2 * A * z
end

```

which is definitely continuous at 0, but the difference between f(1\pm\epsilon) will depend on A, which can be arbitrarily large.

I should have emphasized that the locations of the possible discontinuities are _known_. No need to search for them.

Thanks for all the answers. I will probably go with Richardson extrapolation as it is fairly straightforward, and the cost does not matter (again, this is for CI).

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [April 10, 2024, 11:25am UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/11 "2024-04-10T11:25:30Z")

</div>

The way to go is more symbolic. Like in AD a tape of operations of a function on a “dual variable type” will collect the points of potential discontinuities and their left/right/middle values.

TL;DR copy-paste AD logic and apply to this problem, but `prevfloat`/`nextfloat` is a good approx.

If there is a piecewise-linear functions implementation in some package, it could also be a guide (I’m not familiar with one).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 11, 2024, 11:02am UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/12 "2024-04-11T11:02:19Z")

</div>

> [@Dan](#):
>
> TL;DR copy-paste AD logic and apply to this problem

I may be missing something, but I don’t see how that can be done. Is there a paper I can read about this?

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [April 11, 2024, 12:16pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/13 "2024-04-11T12:16:18Z")

</div>

> [@Tamas\_Papp](#):
>
> I don’t see how that can be done. Is there a paper I can read about this?

It isn’t trivial. And I’m not aware of a paper. The idea would go, by extending a Value into a (Value, Dir) where Dir is a “direction”, either + or - or 0. Most calculation would go through to work on Value, except expressions like `<`,`<=`,`==`,`>` which would virtually add an infinitismal Dir. Thus `(3, +) <= 3` would be false.  
With this setup, finding continuity of `f(3)` would be seeing if `f(3+) == f(3) == f(3-)`. Some other discontinuous functions would also need to be defined, such as `round`, `sign`, etc…

ADDED:  
This “etc…” above does a lot of work. For example `-(3,+)` would be `(-3,-)`.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 11, 2024, 12:18pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/14 "2024-04-11T12:18:30Z")

</div>

I would wager that interval arithmetic (IEEE 1788) is the correct tool for this job. It defines the `DAC` decoration which confirms that a function is ‘defined and continuous’ on an interval (which can be restricted to just a single real number).

Also without decorations one could easily imagine a numerical epsilon-delta test for continuity: As you shrink your input interval around a suspect point on your function, the result interval should also shrink at the appropriate rate.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [April 12, 2024, 2:39pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/15 "2024-04-12T14:39:25Z")

</div>

> [@Tamas\_Papp](#):
>
> I should have emphasized that the locations of the possible discontinuities are _known_. No need to search for them.

Is your setting that your function is stitched together of different definitions at certain nodes, like e.g. a spline? And you now want to do a quick plausibility check that everything matches up?

For a “mathematical / extensional” function defined on floating point numbers, i.e. a blackbox, there is not a lot more to do than comparing `f(prevfloat(p))`, `f(p)` and `f(nextfloat(p))` to some cutoff (which may depend on `f(p)`). What cutoff is appropriate depends on your domain.

For a “mathematical / extensional” function defined on real numbers, I have a halting-problem oracle to sell you. Jokes aside, that doesn’t work ☹

For a AD-able / interval-math-able implementation of a function, there might be some clever definitions and tests. These tests may depend on the implementation of the function, not just its input-output behavior, same as all these things (you can implement f(x)=x as a lookup table on Float32, but that won’t AD).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 12, 2024, 3:02pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/16 "2024-04-12T15:02:11Z")

</div>

> [@foobar\_lv2](#):
>
> Is your setting that your function is stitched together of different definitions at certain nodes, like e.g. a spline? And you now want to do a quick plausibility check that everything matches up?

Yes.

> [@foobar\_lv2](#):
>
> For a “mathematical / extensional” function defined on floating point numbers, i.e. a blackbox, there is not a lot more to do than comparing `f(prevfloat(p))`, `f(p)` and `f(nextfloat(p))` to some cutoff (which may depend on `f(p)`). What cutoff is appropriate depends on your domain.

I am finding Richardson extrapolation easier, especially because the cutoff is easier to determine.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 12, 2024, 3:03pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/17 "2024-04-12T15:03:34Z")

</div>

BTW, I think that Richardson extrapolation _is_ the AD equivalent of taking limits numerically and has a broad range of applications.

Consider the following exercise: numerically verify that

\lim\_{x \to \infty} \log(1 + \exp(x)) - x = 0

> **Solution**
>
> julia\> using Richardson, LogExpFunctions
> 
> julia\> to\_inf(x) = 1/x  
> to\_inf (generic function with 1 method)
> 
> julia\> extrapolate((x → log1pexp(x) - x) ∘ to\_inf, 0.01)  
> (0.0, 0.0)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 12, 2024, 3:21pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/18 "2024-04-12T15:21:24Z")

</div>

> [@Tamas\_Papp](#):
>
> I think that Richardson extrapolation _is_ the AD equivalent of taking limits numerically

AD is exact in exact arithmetic (i.e. in the absence of roundoff error). Richardson extrapolation (with a finite number of iterations) is not. So, they seem qualitatively different to me.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [April 12, 2024, 4:54pm UTC](https://discourse.julialang.org/t/how-to-test-if-a-function-is-continuous/112737/19 "2024-04-12T16:54:10Z")

</div>

> [@Tamas\_Papp](#):
>
> BTW, I think that Richardson extrapolation _is_ the AD equivalent of taking limits numerically and has a broad range of applications.

Richardson is blackbox: You apply it to a function, i.e. a relation.

AD applies to something like a circuit that computes a function, and outputs a circuit that computes / approximates the derivative. The same holds for interval math.