Plots weird problem plotting a piecewise function at discontinuity points

Specific Domains Visualization
,
1 
using Plots;
g_k =0.0:0.01:10.0;
A1=1.0; A2=2.0; α=0.33;
κ2 =3.0;
κ̂2 = κ2 / (1.0 - ((A1/A2)^(1/α))); #3.41842
function f3p(k)
    if (k < κ̂2)
        return α*A1*(k^(α - 1.0))
    elseif (k > κ̂2)
        return α*A2*(( max(0.0, k - κ2) )^(α- 1.0))
    end    
end
plot(g_k, f3p)

Gives

image
image1108×740 7.69 KB

The weird looking vertical line at the discontinuity at 3.41842 should not be there.
How can I remove it?

2

My only hack is to do it manually, but that’s not a good solution for more complicated problems

plot!(0.0:0.01:κ̂2, f3p, lab="f'(k)", c=1, lw=1)
plot!(κ̂2:0.01:10, f3p, lab="f'(k)", c=1, lw=1)

image
image1140×738 19.4 KB

3

Well I don’t think there is an easy solution for exactly what you propose. plot essentially just takes 2 arrays of x and y coordinates and draws lines between them. How should it know which line to omit?

Some other ideas:

  1. Instead of using two plot calls you can mark the jump with the value NaN. IIRC NaNs break the line. So your picewise function could return NaN exactly at the jump (or a small neighbourhood) and if your are plotting enoigh points such that one of them is inside that window, then it should work out.
  2. Instead of plot/plot! you could use scatter/scatter! which just draws dots at the given location. If you play a bit with markersize and number of points you probably can get something reasonable as well.
  3. Maybe there is a package for piecewise functions that can help you?
2 Likes
4

Scatter() looks way too weird.
Here is a more minimal/parsimonious MWE

using Plots;
grid= 0.0:0.0001:1.0
g(x)= (x<0.5) ? (x) : (x-0.5)
plot(grid, g)
scatter!(grid, g)

image
image1138×736 35.3 KB

I’m hoping for some kind of options to ignore discontinuities.

5

Yes you’ll need to play with some setting to make scatter look better. Did you try e.g. with color=:blue, markersize=0.25 or something like that? I am on mobile and can’t play with these settings myself.

Edit: You probably want to disable the border (black line) around the markers. You can try markerstrokewidth=0 for that.

6

It would probably make sense to use the information about the intervals that you have, by incorporating them into a type for which you can make a recipe which then makes gaps at the correct places. If you just have a blackbox function outputting numbers, you cannot really know whether two adjacent points are separated by a singularity or not, could always be continuous on a very small scale I guess?

7

As @jules mention, since this is just some blackbox you would need to provide it yourself with your own recipe, or break it into NaNs. scatter is probably okay too. Here’s one other approach you could take, identifying singularities by large slopes:

using Plots
function fix_singularities(y, x; slope_tol=100)
    _x = similar(x, 0)
    _y = similar(y, 0)
    for i in firstindex(y):(lastindex(y)-1)
        yᵢ, yᵢ₊₁ = y[i], y[i+1]
        xᵢ, xᵢ₊₁ = x[i], x[i+1]
        slope = (yᵢ₊₁ - yᵢ) / (xᵢ₊₁ - xᵢ)
        if abs(slope) < slope_tol
            push!(_x, x[i])
            push!(_y, y[i])
        else
            push!(_x, x[i], (x[i] + x[i+1]) / 2)
            push!(_y, y[i], NaN)
        end
    end
    return _y, _x
end
g(x) = (x < 0.5) ? (x) : (x - 0.5)
x = 0.0:0.0001:1.0
y = g.(x)
_y, _x = fix_singularities(y, x)
plot(_x, _y)

01cb465c1a053dbbbe6591bcb369ebda
01cb465c1a053dbbbe6591bcb369ebda717×479 7.79 KB

Maybe this is robust enough for your use case, if you are willing to recompute the x and y like this (you can of course define your own method for providing a function and a grid still).

1 Like
8

Thanks.
I was wrong to call it a singularity.
I meant discontinuity.

And I STRONGLY believe a good plotting package should have a fixdiscontinuity option.
Just b/c others don’t do it, doesn’t mean we shouldn’t try.

9

Random thought: Maybe Autodiff could identify true discontinuities? Just computing a local gradient does of course not contain this information. I am not too familiar with how all of these different packages work, but maybe one of the approaches could detect when there is a switch to a different code branch or just know due to code analysis when a jump occurs.

1 Like
10

Mathematica plots discontinuous functions correctly

image
image1125×2436 742 KB

11

Yes but this quite different from plotting values from a blackbox. In principle you could define a type in Julia that representa a piecewise function and define a plotting recipe that handles the discontinuity.

1 Like
12

No!
I had no idea my original function was discontinuous or where the discontinuity points were

13

Some simple approach based on the derivative and on dense sampling (try also the tan() case):

using Statistics, Plots; gr(dpi=600)

A1 = 1.0; A2 = 2.0; α = 0.33
κ2 = 3.0
κ̂2 = κ2 / (1.0 - ((A1/A2)^(1/α)));  # 3.41842

f3p(k) = (k < κ̂2) ?  α*A1*(k^(α - 1.0)) : α*A2*((max(0.0, k - κ2))^(α- 1.0))
# f3p(k) = tan(k)               # try this too!

N = 1000                        # the higher/denser the sampling the better
k = range(0, 10, N)
g_k = f3p.(k)

δk = [diff(g_k); 0]
dp = abs.(δk[2:end-1])
ϵ = N * median(dp)              # discontinuity threshold
ix = [false; dp .> ϵ; false]
g_k[ix] .= NaN

plot(k, g_k)
2 Likes
14

The said functions in Mathematica is not the same thing as in Plots.jl. Mathematica knows it is defined piecewise (because of its symbolic nature) and does not assume continuity whereas, in Plots.jl, it is always numeric and it almost never knows it’s symbolic structure unless a recipe with a custom type defines piecewise function structure. Hence, like Matlab and several other plotting libraries, discontinuities in the lines of plot are inferred by NaN values.

3 Likes