6

I'm trying to use Gaussian quadrature to approximate the integral of a function. (More info here: http://austingwalters.com/gaussian-quadrature/). The first function is on the interval [-1,1]. The second function is generalized to [a,b] by change of variable. The problem is that I keep getting the error "'numpy.ndarray' object is not callable". I assume (please correct me if I'm wrong) this means I've tried to call the arrays w and x as functions, but I'm not sure how to fix this.

This is the code

from __future__ import division
from pylab import *
from scipy.special.orthogonal import p_roots

def gauss1(f,n):
    [x,w] = p_roots(n+1)
    f = (1-x**2)**0.5
    for i in range(n+1):
        G = sum(w[i]*f(x[i]))
    return G

def gauss(f,a,b,n):
    [x,w] = p_roots(n+1)
    f = (1-x**2)**0.5
    for i in range(n+1):
        G = 0.5*(b-a)*sum(w[i]*f(0.5*(b-a)*x[i]+ 0.5*(b+a)))
    return G

These are the respective error messages

gauss1(f,4)
Traceback (most recent call last):

  File "<ipython-input-82-43c8ecf7334a>", line 1, in <module>
    gauss1(f,4)

  File "C:/Users/Me/Desktop/hw8.py", line 16, in gauss1
    G = sum(w[i]*f(x[i]))

TypeError: 'numpy.ndarray' object is not callable

gauss(f,0,1,4)
Traceback (most recent call last):

  File "<ipython-input-83-5603d51e9206>", line 1, in <module>
    gauss(f,0,1,4)

  File "C:/Users/Me/Desktop/hw8.py", line 23, in gauss
    G = 0.5*(b-a)*sum(w[i]*f(0.5*(b-a)*x[i]+ 0.5*(b+a)))

TypeError: 'numpy.ndarray' object is not callable

4 Answers 4

5

As Will says you're getting confused between arrays and functions.

You need to define the function you want to integrate separately and pass it into gauss.

E.g.

def my_f(x):
    return 2*x**2 - 3*x +15 

gauss(m_f,2,1,-1)

You also don't need to loop as numpy arrays are vectorized objects.

def gauss1(f,n):
    [x,w] = p_roots(n+1)
    G=sum(w*f(x))
    return G

def gauss(f,n,a,b):
    [x,w] = p_roots(n+1)
    G=0.5*(b-a)*sum(w*f(0.5*(b-a)*x+0.5*(b+a)))
    return G
1

quadpy, a small project of mine, might help:

import numpy
import quadpy


def f(x):
    return numpy.exp(x)


scheme = quadpy.line_segment.gauss_legendre(10)
val = scheme.integrate(f, [0.0, 1.0])
print(val)
1.7182818284590464

There are many other quadrature schemes for 1D.

0

In gauss1(f,n), you are treating f as if it's a function when it is an array, since you are reassigning it;

def gauss1(f,n):
  [x,w] = p_roots(n+1)
  f = (1-x**2)**0.5 # This line is your problem.
  for i in range(n+1):
    G = sum(w[i]*f(x[i]))
  return G

You are doing something similar in the second function.

2
  • If I change it to: def test(f,n): [x,w] = p_roots(n+1) for i in range(n+1): f[i] = (1-x[i]**2)**0.5 G = sum(w[i]*f([i])) return G Then I get the error: File "C:/Users/Me/Desktop/hw8.py", line 29, in test f[i] = (1-x[i]**2)**0.5: TypeError: 'builtin_function_or_method' object does not support item assignment: (sorry about the lack of line breaks....) Nov 24, 2014 at 23:30
  • yes, that's because you're passing a function to gauss1(f,n), with the aim of creating a Gauss-Legendre quadrature approximating it. You should then by calling f(x) inside the function to create the quadrature
    – will
    Nov 24, 2014 at 23:49
0

Example: solving using gaussian integral with n = 2 for integral 2+sinX with b = pi/2 and a = 0

import numpy as np

E = np.array([-0.774597, 0.000000, 0.774597])
A = np.array([0.555556, 0.888889, 0.555556])

def gauss(f, a, b, E, A):
    x = np.zeros(3)
    for i in range(3):
        x[i] = (b+a)/2 + (b-a)/2 *E[i]

    return (b-a)/2 * (A[0]*f(x[0]) + A[1]*f(x[1]) + A[2]*f(x[2]))


f = lambda x: 2 + np.sin(x)
a = 0.0; b = np.pi/2

areaGau = gauss(f, a, b, E, A)
print("Gaussian integral: ", areaGau)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.