I'm trying to convert my matlab code to python for an adaptive gaussian quadrature

Question

Matlab Code-runs fine

function [int, abt]= gadap(a,b,f,p,tol);
% a,b: interval endpoints with a < b
% f: function handle f(x, p) to integrate (p for user parameters)
% tol: User-provided tolerance for integral accuracy
% int: Approximation to the integral
% abt: Endpoints and approximations
a_j=a;
b_j=b;
j=1;
int=0;
n=1;
abt=[[a_j,b_j,0]];
while j<=n
    %%%Evaluation of t%%%
    t_j=(b_j-a_j)/2*(5/9*f((b_j-a_j)/2*-1*sqrt(3/5)+(b_j+a_j)/2)+8/9*f((b_j+a_j)/2)+5/9*f((b_j-a_j)/2*sqrt(3/5)+(b_j+a_j)/2));

    %%%Evaluation of l%%%
    k=j+1;
    a_k=a_j;
    b_k=a_j+(b_j-a_j)/2;
    l_j=(b_k-a_k)/2*(5/9*f((b_k-a_k)/2*-1*sqrt(3/5)+(b_k+a_k)/2)+8/9*f((b_k+a_k)/2)+5/9*f((b_k-a_k)/2*sqrt(3/5)+(b_k+a_k)/2));

    %%%Evaluation of r%%%
    z=j+2;
    a_z=a_j+(b_j-a_j)/2;
    b_z=b_j;
    r_j=(b_z-a_z)/2*(5/9*f((b_z-a_z)/2*-1*sqrt(3/5)+(b_z+a_z)/2)+8/9*f((b_z+a_z)/2)+5/9*f((b_z-a_z)/2*sqrt(3/5)+(b_z+a_z)/2));

    %%%List Generation%%%
    if abs(t_j-(l_j+r_j))>tol*max(abs(t_j), (abs(l_j)+abs(r_j)))
        abt=[abt; [a_k, b_k, l_j]; [a_z, b_z, r_j]];
    else 
        int=int+t_j;
    end

    n=size(abt,1);
    j=j+1;
    if j>n
        continue
    end
    a_j=abt(j,1);
    b_j=abt(j,2);


end

Python Code

import numpy as np
from numpy import *
from numpy import array


f = lambda x:x**2
def gaussian(a,b,f,tolerance):
    # a,b: interval endpoints with a < b
    # f: function f(x) to integrate 
    # tol: tolerance for integral accuracy
    # integral: Approximation to the integral
    # abt: Endpoints and approximations
    a_j=a
    b_j=b
    j=0
    integral=0
    n=0
    abt=matrix([a_j,b_j,0])
    while j<=n:
        #Evaluation of t
        t_j=(b_j-a_j)/2*(5/9*f((b_j-a_j)/2*-1*sqrt(3/5)+(b_j+a_j)/2)+8/9*f((b_j+a_j)/2)+5/9*f((b_j-a_j)/2*sqrt(3/5)+(b_j+a_j)/2))

        #Evaluation of right
        z=j+2;
        a_z=a_j+(b_j-a_j)/2;
        b_z=b_j;
        r_j=(b_z-a_z)/2*(5/9*f((b_z-a_z)/2*-1*sqrt(3/5)+(b_z+a_z)/2)+8/9*f((b_z+a_z)/2)+5/9*f((b_z-a_z)/2*sqrt(3/5)+(b_z+a_z)/2))

        #Evaluation of left
        k=j+1;
        a_k=a_j;
        b_k=a_j+(b_j-a_j)/2;
        l_j=(b_k-a_k)/2*(5/9*f((b_k-a_k)/2*-1*sqrt(3/5)+(b_k+a_k)/2)+8/9*f((b_k+a_k)/2)+5/9*f((b_k-a_k)/2*sqrt(3/5)+(b_k+a_k)/2))

        if abs(t_j-(l_j+r_j))>tolerance*max(abs(t_j), (abs(l_j)+abs(r_j))):
            abt=numpy.vstack([abt, [a_k, b_k, l_j], [a_z, b_z, r_j]])
        else:
            integral=integral+t_j
        n=abt.shape[0]
        j=j+1
        if j>n:
            continue
        a_j=abt[j,0]
        b_j=abt[j,1]
    return integral

print gaussian(0,2,f,10**(-3))

But when I test the python code for my function I get an error. What's wrong?want to return integral value,I changed some indices but it keeps returning the integral value of 0, which was what I initialized the integral with. Now it tells me that the error is

IndexError Traceback (most recent call last) in () 45 return integral 46 ---> 47 print gaussian(0,2,f,10**(-3)) 48

in gaussian(a, b, f, tolerance) 41 if j>n: 42 continue ---> 43 a_j=abt[j,0] 44 b_j=abt[j,1] 45 return integral

C:\Users\Brandon Tran\AppData\Local\Enthought\Canopy\User\lib\site-packages\numpy\matrixlib\defmatrix.pyc in getitem(self, index) 314 315 try: --> 316 out = N.ndarray.getitem(self, index) 317 finally: 318 self._getitem = False

IndexError: index 1 is out of bounds for axis 0 with size 1

Answer 1

No need for implementing it yourself. quadpy (which I authored) has adaptive quadrature.

Install with

pip install quadpy

and try

from numpy import exp
import quadpy

val, err = quadpy.quad(lambda x: exp(0.3 * x), 0.0, 1.0)
print(val)
print(err)

1.1661960252533439
3.170663072723789e-20

The code implements adaptivity via Gauss-Kronrod.

Answer 2

Your problem is that MATLAB and numpy use different ordering for array indices. MATLAB uses the order from the Fortran programming language, while numpy (by default) uses the ordering from the C language. These are opposites. So while your MATLAB matrix is an Nx1 matrix, your numpy matrix is a 1xN matrix. So, for example, a_j=abt[j,0] should be a_j=abt[0,j]. You can verify this with abt.shape, which will gives you (1, 3).

There are other problems as well. In Python 2.x, division of two integers returns an integer. So 5/9, for example, is 0. At the top of your code you should put from __future__ import division to make it behave the way you expect.

Also, you really shouldn't be using numpy matrices to begin with. Stick with numpy arrays, they are much better supported. That way you can avoid your problem entirely by just using a 1D array.

Finally, it is much safer to do import numpy as np rather than from numpy import *, which can produce function name collisions with the built-in versions.