I'm trying to use Gaussian quadrature to approximate the integral of a function. (More info here: ). 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
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 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.
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.
2Example: 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 = areaGau = gauss(f, a, b, E, A) print("Gaussian integral: ", areaGau) 
