Integration with Riemann Sum Python

Question

I have been trying to solve integration with riemann sum. My function has 3 arguments a,b,d so a is lower limit b is higher limit and d is the part where a +(n-1)*d < b. This is my code so far but. My output is 28.652667999999572 what I should get is 28.666650000000388. Also if the input b is lower than a it has to calculate but I have solved that problem already.

def integral(a, b, d):
    if a > b:
        a,b = b,a
    delta_x = float((b-a)/1000)
    j = abs((b-a)/delta_x)
    i = int(j)
    n = s = 0
    x = a
    while n < i:
        delta_A = (x**2+3*x+4) * delta_x
        x += delta_x
        s += delta_A

        n += 1

    return abs(s)

print(integral(1,3,0.01))

Your implementation does not use the d variable. Is this intentional? — MrHug
– MrHug, Commented Mar 28, 2015 at 15:15
3.4 Python version. Yes i don't know how to use d correctly so i did first without it. — Hybr1d
– Hybr1d, Commented Mar 28, 2015 at 15:25
Did you try to divide by more than 1000? Maybe 1e6 where you define delta_x. — Daniel Thaagaard Andreasen
– Daniel Thaagaard Andreasen, Commented Mar 28, 2015 at 15:53
@Daniel Thaagaard Andreasen, yes it's closer than... but still not what it should be if i give more 1e8 or more then the function takes too long. =/ — Hybr1d
– Hybr1d, Commented Mar 28, 2015 at 16:01

mbomb007 · Accepted Answer · 2017-02-06 22:33:35Z

2

There is no fault here, neither with the algorithm nor with your code (or python). The Riemann sum is an approximation of the integral and per se not "exact". You approximate the area of a (small) stripe of width dx, say between x and x+dx, and f(x) with the area of an rectangle of the same width and the height of f(x) as it's left upper corner. If the function changes it's value when you go from x to x+dx then the area of the rectangle deviates from the true integral.
As you have noticed, you can make the approximation closer by making thinner and thinner slices, at the cost of more computational effort and time. In your example, the function is f(x) = x^2 + 3*x + 4, and it's exact integral over x in [1.0,3.0) is 28 2/3 or 28.66666...

The approximation by rectangles is a crude one, you cannot change that. But what you could change is the time it takes for your code to evaluate, say, 10^8 steps instead of 10^3. Look at this code:

def riemann(a, b, dx):
    if a > b:
        a,b = b,a
    # dx = (b-a)/n
    n = int((b - a) / dx)
    s = 0.0
    x = a
    for i in xrange(n):
        f_i = (x + 3.0) * x + 4.0
        s += f_i
        x += dx
    return s * dx

Here, I've used 3 tricks for speedup, and one for greater precision. First, if you write a loop and you know the number of repetions in advance then use a for-loop instead of a while-loop. It's faster. (BTW, loop variables conventionally are i, j, k ... whereas a limit or final value is n). Secondly, using xrange instead of range is faster for users of python 2.x. Thirdly, factorize polynoms when calculating them often. You should see from the code what I mean here. This way, the result is numerically stable. Last trick: operations within the loop which do not depend on the loop variable can be extracted and applied after the loop has ended. Here, the final multiplication with dx.

edited Feb 6, 2017 at 22:33

mbomb007

4,3913 gold badges51 silver badges79 bronze badges

answered Mar 28, 2015 at 22:23

user1016274

4,2992 gold badges25 silver badges20 bronze badges

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

5 Comments

Hybr1d Over a year ago

Tnx, you so much :) and nice explanation :).

user1016274 Over a year ago

You're welcome! I hope you go on fom here after a while and use trapezes instead of rectangles to approximate the integral (called Simpson's rule). I started doing num-int on a calculator and found it very rewarding. Have fun!

mbomb007 Over a year ago

I think the answer is still wrong somehow. You need to multiply by dx twice. You divided by it earlier, and so the first time will undo that, but you still have to multiply each value in the sum by dx to get the tiny rectangle. tio.run/nexus/python2#fY/… . To fix it, do s += f_i * dx

user1016274 Over a year ago

No. We divide the interval over which we integrate into n equal pieces. As we're given dx instead of n, we calculate n from dx (not the other way around). This is the division by dx. Still, the sum is over f_i * dx where dx is a constant so we can extract it from the sum and multiply it at the end.

user1016274 Over a year ago

Actually, looking at the code you'll notice that at each iteration (of n) a constant (4.0) is added. We can leave out the n additions and add n * 4.0 after the loop has finished, slicing off yet another microsecond. Just to demonstrate the principle of avoiding operations with constants in loops if possible.

Collectives™ on Stack Overflow

Integration with Riemann Sum Python

1 Answer 1

5 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related