numerical integration methods question python

Question

I am trying to integrate the following formula:

Below is my attempt to perform this integration using a homemade scheme for f(a) = sin(a).

def func(x):
    return math.sin(x)

def integration (f, n, r, a, dtheta ):

    summation = 0
    theta = 0
    while theta <= 2*np.pi:

        f_arg = a + r*np.exp(1j*theta)
        second = np.exp(-1j*theta*n)

        summation += f(f_arg) * second * dtheta
        theta += dtheta

    return math.factorial(n)*summation / (2*np.pi*r**n)

integration(func, n=1, r=1, a=0, dtheta=2*np.pi/10000)

The first derivative (n=1) of f(a) = sin(a) is f'(a) = cos(a). When evaluated at a = 0, this should give cos(0) = 1, however, it does not. Where am I going wrong?

return math.factorial(n) / (2*np.pi*r**n) doesn't seem to involve f, theta, dtheta, or anything you did in the while loop. — user2357112
– user2357112, Commented Nov 12, 2018 at 20:00
@user2357112, thank you for attention, updating question (as that was a typo), still didn't fix. — RedPen
– RedPen, Commented Nov 12, 2018 at 20:01

Matthias Ossadnik · Accepted Answer · 2018-11-12 20:29:52Z

3

It seems that your problem is the math.sin function, which doesn't support complex arguments:

i = np.exp(.5j * np.pi)
math.sin(i), np.sin(i)
(6.123233995736766e-17, (9.44864380126377e-17+1.1752011936438014j))

It also throws a warning (but not an error...):

ComplexWarning: Casting complex values to real discards the imaginary part

Using np.sin instead fixes the problem.

In general, the implementation might be simpler to express (and easier to debug) with more use of numpy, like

def integration(func, a, n, r, n_steps):
    z = r * np.exp(2j * np.pi * np.arange(0, 1, 1. / n_steps))
    return math.factorial(n) * np.mean(func(a + z) / z**n)

np.allclose(1., integration(np.sin, a=0., n=1, r=1., n_steps=100))

True

edited Nov 12, 2018 at 20:29

answered Nov 12, 2018 at 20:22

Matthias Ossadnik

9114 silver badges8 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

user2357112 Over a year ago

Replacing math.sin with np.sin isn't enough to produce a result of 1 when I try it.

Matthias Ossadnik Over a year ago

@user2357112 in the original example there is also n=2 used instead of n=1 in the function call. If you change this it is OK. Have edited question.

Collectives™ on Stack Overflow

numerical integration methods question python

1 Answer 1

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related