Integrating over piecewise approxfun

Specific Domains Numerics
,
1

Following up on my earlier post, I now have an ApproxFun.jl Fun defined over a piecewise domain:

using ApproxFun
f_pw = Fun(x -> (x < 1.0 ? exp(x) : 1.), union(Segment(0.,1.), Segment(1.,2.)))
f_pw_i = integrate(f_pw)

The integrated f_pw_i does not quite behave like the integral, as each piece is separately centered to yield a zero constant term.
Is there a way to define an operator that integrates from the leftmost boundary onwards, and over the discontinuity to yield a continuous function on the joint domain of f_pw?

2

Try cumsum? Also I think there’s a PiecewiseSegment.

1 Like
3

Ah yes that’s better.

I was using integrate rather than cumsum because in the actual application in need an integral over one dimension of a two-dimensional Fun. AFAIR cumsum did not work for that when I tried it. Currently I am doing that like this (and that’s where I ran into the issue):

Integral(factor(space(twodim_fun), 1)) ⊗ I

Would cumsum generalize like this?

Also, I recall PiecewiseSegment not combining well into a two-dimensional domain, but I can retry to be sure.

4

Maybe the confusion stems from the fact that integrate gives you the antiderivative. So to get the function that is the integral from the leftmost boundary you need to subtract that:

antiderivative = integrate(f_pw)
integral_from_left = x -> antiderivative(x) - antiderivative(0) # IIUC that 0 is left boundary

I know that is confused me for a short while when starting with ApproxFun.jl :sweat_smile:

1 Like
5

That’s useful, but it covers only the left boundary, not the jump in the middle.

6

here is a non-working example:

d2 = cross(Segment(0.,1), PiecewiseSegment(0.,1.,2.))
f_pw2 = Fun((t, x) -> t * (x < 1.0 ? exp(x) : 1.), d2)

where the product domain with one piecewise factor can be constructed but the Fun cannot.

7

I think you are pushing ApproxFun.jl beyond what it’s been tested on…

It’s possible that PiecewiseOrthogonalPolynomials.jl might be better suited to your needs, eventually, though it’s pretty half-baked at the moment. There’s hidden support for 2D (it allows quasi-optimal complexity hp-FEM in a rectangle ) but it’s not cleaned up in a user-friendly way.

In the future the following would work:

julia> using PiecewiseOrthogonalPolynomials, ClassicalOrthogonalPolynomials

julia> f = expand(ContinuousPolynomial{0}([0,1,2]), x -> x < 1 ? exp(x) : 1.0)
ContinuousPolynomial{0}([0, 1, 2]) * [1.718281828459045, 1.0, 0.8451545146228638, 0.0, 0.13986399606658317, 0.0, 0.013931255854518054, 0.0, 0.0009925875385252538, 0.0  …  ]

julia> cumsum(f)
ERROR: Not implemented
2 Likes
8

too bad… It’s almost working already:

d3 = cross(Segment(0.,1), union(Segment(0.,1.), Segment(1., 2.)))
f_pw3 = Fun((t, x) -> t * (x < 1.0 ? exp(x) : 1.), d3)

is accepted – only acting on it with I ⊗ Integral(factor(d3, 2)) has the jumpy behavior that motivated the topic.

9

After way too much spelunking and trying to coerce operators to combine, I stumbled upon cumsum(::ProductFunction, n) which does integration over one direction in a multidimensional function, with correct handling of piecewise dimensions! I could not imagine it being that easy before, unfortunately.

I think I will get by using that, like so:

cumsum(ProductFun(f_pw3), 2)

Thanks everyone.