Integrating e^-(x^2)

[SK6114]

Hi,

I need to integrate e^-(x^2) between negative infinity and positive infinity, now the function vanishes to a minute value much before negative and positive infinity.

I need to approximate the integral using the trapezium rule, which I have coding for but I just don't know how to apply it to the function.

Here is my coding for the trapezium rule:

def trap1 (f,a,b,delta , maxtraps=512):
          n=8
          inew= trap0(f,a,b,n)
          iold=- inew 
          while ( fabs(inew - iold)>delta * fabs( inew )):
            iold= inew
            n=2*n
            if n> maxtraps:
              print " Cannot reach requested accuracy with", \
                     maxtraps , " trapezia"
                     return
               inew= trap0 (f,a,b,n)
           return inew

Note sure if I have the indentation right; 'n' is the number of trapeziums or widths I take, delta is just a small number I specify.

I basically need to apply the function f as e^-(x^2) and I need to alter the trapezium coding such that it gives me a good approximate for the integral - so I'm guessing I take some -100 to 100 as my 'a' and 'b' values.

Is the above coding the correct basis for the trapezium rule as well?

Thanks!

[Gribouillis]

Here is your code with the correct indentation

def trap1 (f,a,b,delta , maxtraps=512):
    n=8
    inew= trap0(f,a,b,n)
    iold=- inew 
    while ( fabs(inew - iold)>delta * fabs( inew )):
        iold= inew
        n=2*n
        if n> maxtraps:
            print " Cannot reach requested accuracy with", \
                maxtraps , " trapezia"
            return
        inew= trap0 (f,a,b,n)
    return inew

Indent python code with 4 spaces. Since you didn't post the function trap0(), it's difficult to say if your code will work. You can call the function by defining your own function

from math import exp

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

result = 2 * trap1(myfunc, 0, 6, .001)

I suggest [0, 6] for a and b because when A >= 0.5

int_A^infty exp(-x^2) dx <= int_A^infty 2x exp(-x^2) dx = exp(-A^2)

and

>>> exp(-6*6)
2.3195228302435696e-16

[SK6114]

I am trying to calculate the approximate integral of a function using the trapezium rule; the code below shows what I have so far:

def f(x):
    (x*x*x*x*(1-x)**4)/(1+x*x)

def trap0 (f,a,b,n):
    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

The function f(x) is supposed to be ((x^4)(1-x)^4)/(1+(x^2)), I think I have typed it in correctly. Now I want to calculate an approximate integral using the trapezium rule between the limits 1 and 0, 'n' is just the number of widths I use.

I believe I should state the values of a,b and n somewhere but I have tried that and I keep getting this when I type print trap0 <function trap0 at 0x01DD1330>.

What am I doing wrong?

[TrustyTony]

Why wrong? You have the function defined at line 4.

[SK6114]

I put the trap0 function at line 1 and I still receive the same error, where do I state my values of a,b and n?

def trap0 (f, a, b, n):
    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 f(x):
   return (x*x*x*x*(1-x)**4)/(1+x*x)

a=0.0
b=1.0
n=50.0

print trap0

[TrustyTony]

You are not calling function, you are printing function objects string representation. Give function parameters in ()