## 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 104

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 2,720

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.