Numerical Evaluation of Discontinuous Integrals

[Nathan Woods]

Hi everyone,

As part of some work I was doing, I had to develop a method for efficiently evaluating integrals of multidimensional, piece-wise-continuous functions based on a priori knowledge of where the jumps would happen. The resulting code takes a Sympy expression and processes it into a SciPy numerical integration call. The interesting bit is that the resulting integration call exactly fits the discontinuity(s), so the integrator only has to deal with the smooth regions, rather than trying to progressively refine a polynomial fit to a step function. Accomplishing that relies heavily on symbolic math manipulations and the code generation tools in Sympy.

Anyway, I would like to package this up in a way that would be publicly useful, but I'm not sure where it fits. Sympy seemed a likely guess, but the SciPy dependence is problematic. Alternatively, the discontinuity-processing could be decoupled and used with any iterated integrator that supports manually specified points of discontinuity. If Sympy isn't a good fit for this, I would appreciate suggestions about other places that might be better.

Thanks for the help,

Nathan Woods

[Joachim Durchholz]

Am 24.02.2015 um 17:03 schrieb Nathan Woods:
> Anyway, I would like to package this up in a way that would be publicly
> useful, but I'm not sure where it fits. Sympy seemed a likely guess, but
> the SciPy dependence is problematic.

Out of the box, a separate project with a dependency on both SymPy and
SciPy would probably fit best.

> Alternatively, the
> discontinuity-processing could be decoupled and used with any iterated
> integrator that supports manually specified points of discontinuity.

Decoupling and putting the modules into their respective projects would
work, too.
It's much more work though.

> If
> Sympy isn't a good fit for this, I would appreciate suggestions about other
> places that might be better.

I did not fully understand your description of what your code does, so I
can't be very specific.
In general, anything that analyzes mathematical objects for properties
of interest that go beyond a one-shot task would be a worthy addition to
SymPy.

[Nathan Woods]

The decoupling is actually pretty easy. Unfortunately, the Sympy-specific stuff isn't especially useful without the context of a numerical integrator. 

Maybe an example will help. Say you want to numerically integrate a step function H(x) over the interval [-1, 1]. Most numerical integrators will try to fit a polynomial to this, but polynomials can't handle the jump at x=0 gracefully, so you end up wasting a lot of processing time trying to refine it. The fix is to break up the integral, and do it in two pieces, one from [-1, 0), and one from [0, 1]. The integrator has no trouble with this, so you get a very fast evaluation. Many integrators, and in particular scipy.integrate.quad, have this capability.

What I've done is figure out a way to get the same benefits for multidimensional functions. So, you could integrate H(x**2 + y**2 - z) over a volume, and get the same kind of efficiency. Without this, not only is evaluation of the integral very slow, but the results can sometimes be very inaccurate too. 

Anyway, from what I can tell, a similar effect could be achieved using mpmath, rather than SciPy. Maybe that project is a better venue. 

Thanks again!

N

[Jason Moore]

Nathan,

Can you show us the code? It may help us understand what you are doing.

Also, we have code in sympy that optional depends on scipy as do we other packages: cython, numpy, theano, matplotlib, etc. The code generation and tightly coupled symboli/numeric code is in that blurry zone about what we'd add.

But if you make a working separate package it can be a good start. For example the PyDy project depends on sympy but deals with generating specific numerical codes. If you make the package and we want to include it, then great. If not, you have an easily installable package that others can use and we can link to it from our docs/website.

Jason
moorepants.info
+01 530-601-9791

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/125f4bc9-7c1d-441f-bf41-5df851379a7d%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Nathan Woods]

I'll need to break it out of a more specialized package first, but this all sounds very promising. It'll be a few weeks before I can really get started (dissertation writing...), but I think it should only take a day or so after that, at least to have something to share. Maybe longer to flesh out the feature set and add a cleaner interface. Mostly I just want it to be available in a place where it'll get use.

N

[Aaron Meurer]

This sounds like it could be useful for lambdify. I'd like to see
lambdification of integrals be improved.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/bed69b8c-fd47-4cba-a3fa-e800e34acea8%40googlegroups.com.

[Nathan Woods]

All right, it took me longer than I had anticipated, but I have this project assembled in a more complete form now. To recap, the idea is to simplify allow for accurate evaluation of integrals of discontinuous functions, for which the form of the discontinuity is known beforehand. There is a paper describing the concept available here: http://authors.elsevier.com/a/1QsPV508HRs4u (this link expires June 2). I have also broken out the integration software, which is available in my GitHub: https://github.com/woodscn/sympy_discontinuous_integration

As it stands right now, this package uses symbolic processing to convert multidimensional integrals into the correct form for SciPy's nquad integration routine. For discontinuous integrand functions, it also fits the multidimensional discontinuity, which eliminates the error that would normally be introduced into the integral. I see no reason why this idea could not be adapted to use the built-in Sympy integration routines; that's just not how I developed it. Please let me know what you think, and whether this would be a good fit for Sympy.

Nathan Woods