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.
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 =
areaGau = gauss(f, a, b, E, A)
print("Gaussian integral: ", areaGau)
