Trapezium rule for exp(-x*x)

[John_47]

Hello, i'm currently trying to write some code to find the integral of exp(-x*x).

I have been given code for one trapezium rule (trap0) which is imperfect, and a second set of code for a new trapezium rule (trap1), which uses the old code and some new conditions to 'perfect' it. Here is my attempt.

import numpy 
import matplotlib 
import pylab as P
import math
from math import exp

def f(x):
    return exp(-x*x)

def trap0 (f,a,b,n):
    """Basic trapezium rule. Integrate f(x) over the 30 4 Numerical Integration interval from a to b using n strips."""
    h= float (b-a)/n
    s =0.5*( f(a)+f(b))
    for i in range (1,n):
        s=s+f(a+i*h)
        return s*h


def trap1 (f,a,b,delta , maxtraps=512):
    """ Improved trapezium rule. Integrate f(x) over interval from a to b, trying to get relative accuracy delta. Optional last argument is maximum allowed number of trapezia."""
    n=8
    inew= trap0(f,100,-100,n)
    iold=- inew # make sure we don’t terminate immediately
    while ( fabs(inew - iold)>delta * fabs( inew )):
        iold= inew
        n=2*n
        if n> maxtraps:
            print " Cannot reach requested accuracy with"

The integral rapidly becomes negligible if x is large so i understand you would simply use the limits of -5 to 5 to calculate this integral. However, i'm still unsure of how to do this.

I am using python 2.7

Thanks for any help

[BearofNH]

The integral rapidly becomes negligible if x is large so i understand you would simply use the limits of -5 to 5 to calculate this integral. However, i'm still unsure of how to do this.

Unless I'm mistaken, I think you can fix this by inserting between line 11 and 12:

# Assume a<b since it's our code
if a >  5: return 0
if a < -5: a = -5 
if b < -5: return 0
if b >  5: b =  5

LFTR: for completeness, you may also want to handle the case of a>b. See if you can do that in one line.

I don't see how your fabs() calls work. They should be math.fabs(), or you can include fabs with exp on the import line.

I notice you don't change inew inside the while loop. That doesn't look right, but maybe it's just the name. E.g., if it represents the "best" integral value it should be called ibest. Using names iold and inew sort of implies you are changing inew each time thru the while loop, and saving the prior value in iold, which isn't happening.

[ddanbe]

In this case, you could take advantage of the fact that the curve is symetric around the y-axis.

[John_47]

Thanks for the help, your additions worked BearofNH. Indentations were slightly wrong as well.