What is wrong with my implementation of numerical integration?

Question

I think similar questions have answered here and here. I am also having similar problems with a code that performs the trapezoidal and simpson's 1by3 rule. I have ran the code on a desktop which runs fine, but using a different machine, the code fails. I am sharing the code and the error message.

Code

import numpy as np

exp = 5.33333
ul = 10
ll = -10
n = 50

def f(x):                           # The function
    return np.e**x**2

def trapezoidal(ll, ul, n):           #Setting up Trapezoidal rule
    traint = (f(ll) + f(ul)) / 2
    h = (ul - ll) / 2
    for i in range(1,n):
        k = ll + i * h
        traint = traint + f(k)
    traint = h * traint
    return traint

def simpson1by3(ll, ul, n):           #Setting up Simpson 1by3 rule
    simint = (f(ll) + f(ul)) / 2
    h = (ul - ll) / 2
    for i in range(1,n,2):
        simint = simint + (2 * f(ll + i * h)) + (4 * f(ll + ((i + 1) * h)))
    simint = simint * h / 3
    return simint

res_sim = simpson1by3(ll, ul, n)
res_trap = trapezoidal(ll, ul, n)

e_rel_trap = (res_trap - exp) / exp
e_rel_simp = (res_simp - exp) / exp

f1 = open("trap.txt", "w+")
f2 = open("simp.txt", "w+")

for n in range(1,50):
    h = (ul - ll) / 2
    res_trap = trapezoidal(ll, ul, n)
    res_simp = simpson1by3(ll, ul, n)
    f1.write(str(n) + " " + str(res_trap) + '\n')
    f2.write(str(n) + " " + str(res_simp) + '\n')

f1.close()
f2.close()

The error message is,

Traceback (most recent call last):
  File "simpson1by3+trapezoidal.py", line 30, in <module>
    res_sim = simpson1by3(ll, ul, n)
  File "simpson1by3+trapezoidal.py", line 26, in simpson1by3
    simint = simint + (2 * f(ll + i * h)) + (4 * f(ll + ((i + 1) * h)))
  File "simpson1by3+trapezoidal.py", line 11, in f
    return np.e**x**2
OverflowError: (34, 'Numerical result out of range')

Your function `f` increases very rapidly. Just to confirm, is `np.e**x**2` really what you want? — cmbfast, Aug 17 '21 at 07:21
@cmbfast Good point. Perhaps the OP means to integrate a normal distribution-like function, where it would be `math.e**(-x**2)` (using `math` instead of `numpy`, since NumPy is unnecessary here). — 9769953, Aug 17 '21 at 07:25
@cmbfast To my knowledge, I do not know the difference that makes between np.e**x**2 or math.e**x**2. — Prasanta Bandyopadhyay, Aug 17 '21 at 08:41
I can use math but indeed want to evaluate ```e**x**2``` not ```e**(-x**2)```. — Prasanta Bandyopadhyay, Aug 17 '21 at 08:46
@PrasantaBandyopadhyay power operator has right-associativity so `e**x**2` is treated as `e**(x**2)`. And **never** use `e**x` to get `e^x`, always use `math.exp(x)` — phuclv, Aug 17 '21 at 11:23
There is no different between `np.e` and `math.e`, except that the former involves installing and loading NumPy. Not a great bother, but an extra nonetheless. — 9769953, Aug 17 '21 at 12:20
What phyclv likely misses to state, is that `math.exp(x)` (or the equivalent `np.exp`) is very likely optimised, compared to using `math.e**x`. — 9769953, Aug 17 '21 at 12:22

score 1 · Accepted Answer · answered Aug 17 '21 at 07:47

1

Seems your integration interval is wrongly set. If I didn't remembered wrongly the interval shall be like (b-a)/2n, i.e., ul-ll/n or ul-ll/2n depends on your setting but not ul-ll/2. Setting the interval to ul-ll/2 makes your actual integral upper limit too large, and for a function like e**x**2, it overflows out of max Pythonic float range (around 1.8e308 == 2^1024 == FP256 MAX on my machine)

answered Aug 17 '21 at 07:47

KazamiMichiru

103
7

It is true about the replacement of n with 2. However, if I may ask, How to check the pythonic float range? Does it vary with the version or computer architechture? – Prasanta Bandyopadhyay Aug 17 '21 at 08:35
2

To the best of my knowledge, it might be the biggest CPU supporting bitwidth of floating point numbers. use `sys.float_info.max` to get the max on your machine – KazamiMichiru Aug 17 '21 at 08:38

Duffy · Answer 2 · 2021-08-17T07:30:59.403

0

This mainly due to the limitation of your local machine. You see, python allocates only some amount of space for integers and floats Its usually 4 for int and 8 for float.

in you code, you have set the data type as float with the f

You computer cannot handle such calculations that is taking 16 bytes of memory, If you are sure that your calculations or input would not have any decimals in it, you can change the data type to integers.. for the inputs, If it still doesnt work, consider also changing the answers's data type/

Hope it is helpfull

edited Aug 17 '21 at 07:30

answered Aug 17 '21 at 07:17

Duffy

123
7

8 bytes for a float: 64 bits equals 8 bytes. Not 16 bytes. – 9769953 Aug 17 '21 at 07:20
Suggesting to change the data type to integers here is a bad idea, because that loses far too much precision for this calculation. And wouldn't work, since various calculations would turn the integer(s) into floats anyway (such as the exponentiation). – 9769953 Aug 17 '21 at 07:23
its only when he is sure that the calculations does not involve floats, – Duffy Aug 17 '21 at 07:32
Just look at the code: there's an exponent involved with base `e`, that is already a float. – 9769953 Aug 17 '21 at 12:19

What is wrong with my implementation of numerical integration?

2 Answers2