Fixing the singularity of a function

Question

Assume you have a function like

F = lambda x: sin(x)/x

Evaluating F(0.0) would result in a divide by zero warning, and would not give the expected result of 1.0. Is it possible to write another function fix_singularity that would give the desired result when applied to the above function, so that

fix_singularity(F)(0.0) == 1.0

Or formally fix_singularity should pass the following tests:

import numpy as np

def test_fix_singularity():

    F = lambda x: np.sin(x)/x

    x = np.array((0.0, pi))

    np.testing.assert_array_almost_equal( F(x), [nan, 0] )

    np.testing.assert_array_almost_equal( fix_singularity(F)(x), [1, 0] )

One possible implementation is

def fix_singularity(F):
    """ Fix the singularity of function F(x) """

    def L(x):
        f = F(x)
        i = np.isnan(f)
        f[i] = F(x[i] + 1e-16)
        return f

    return L

Are there better ways of doing this?

EDIT: Also how can I suppress the warning:

Warning: invalid value encountered in divide

lafras · Accepted Answer · 2011-04-05T19:38:06.080

7

numpy has a sinc() function, which is the normalised form of your function, i.e.

F = lambda x: sin(pi*x) / (pi*x)

It handles the case for x == 0.0 correctly,

In [16]: x = numpy.linspace(-1,1,11)

In [17]: print x
[-1.  -0.8 -0.6 -0.4 -0.2  0.   0.2  0.4  0.6  0.8  1. ]

To "unnormalize" do,

In [22]: s = numpy.sinc(x/numpy.pi)

In [23]: print s.round(2)
[ 0.84  0.9   0.94  0.97  0.99  1.    0.99  0.97  0.94  0.9   0.84]

edited Apr 05 '11 at 19:38

answered Apr 05 '11 at 19:08

lafras

8,712
4
29
28

+1: I actually think this is the best answer for the OP's specific case. I hadn't come across this method before. Good work. – JoshAdel Apr 05 '11 at 19:34

JoshAdel · Answer 2 · 2011-04-05T19:22:41.047

3

If you are already using numpy then:

a = np.linspace(0.0,2*np.pi,100)
b = np.sin(a)/a

Will calculate without error leaving a NaN value in b[0]. You could then just replace it by the following if that's how you want to handle it:

b[np.isnan(b)] = 1.0

Update To suppress the warning, try:

np.seterr(divide='ignore') # Or possibly np.seterr(invalid='ignore')

edited Apr 05 '11 at 19:22

answered Apr 05 '11 at 18:52

JoshAdel

66,734
27
141
140

I don't understand your reference to fallacies. OP is referring to the fact that sin(x)/x -> 1 as x -> 0 – David Heffernan Apr 05 '11 at 19:15
@David: Yeah, should have applied L'Hospital before I answered that. Corrected. – JoshAdel Apr 05 '11 at 19:24

score 2 · Answer 3 · answered Apr 05 '11 at 18:59

In general you can't write a simple fix decorator as you might imagine. For example, a general function need not have a finite limiting value at a singularity as this particular example does.

Normal practice is to implement special handling on a case by case basis.

score 1 · Answer 4 · answered Apr 05 '11 at 19:04

1

I'll try this

>>> def fix_singularity(F):
...     def L(x):
...         x1 = max(x,1e-16) if x >=0 else min(x,-1e-16)
...         return F(x1)
...     return L
...
>>> FS = fix_singularity(F)
>>> FS(0.0)
1.0
>>> FS(-1e-17)
1.0

answered Apr 05 '11 at 19:04

Wai Yip Tung

18,106
10
43
47

score 0 · Answer 5 · answered Apr 05 '11 at 18:59

0

I don't know if this would work for your exact purposes, but there's a python library called sage that can handle quite a bit of Calculus-type situations.

answered Apr 05 '11 at 18:59

Matt Howell

15,750
7
49
56

score 0 · Answer 6 · answered Apr 05 '11 at 19:07

0

I believe sympy (symbolic python) can do limits, which is what you are really asking (that solution is only true as a limit). Regardless, you should check it out.

answered Apr 05 '11 at 19:07

Garrett Berg

2,585
1
22
21

Fixing the singularity of a function

6 Answers6