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

matlab code-runs fine

function [int, abt]= gadap(a,b,f,p,tol); % a,b: interval endpoints < b % f: function handle f(x, p) integrate (p user parameters) % tol: user-provided tolerance integral accuracy % int: approximation integral % abt: endpoints , 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 np numpy import * numpy import array   f = lambda x:x**2 def gaussian(a,b,f,tolerance):     # a,b: interval endpoints < b     # f: function f(x) integrate      # tol: tolerance integral accuracy     # integral: approximation integral     # abt: endpoints , 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 test python code function error. what's wrong?want return integral value,i changed indices keeps returning integral value of 0, initialized integral with. tells me 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 out of bounds axis 0 size 1

your problem matlab , numpy use different ordering array indices. matlab uses order fortran programming language, while numpy (by default) uses ordering c language. these opposites. while matlab matrix nx1 matrix, numpy matrix 1xn matrix. so, example, a_j=abt[j,0] should a_j=abt[0,j]. can verify abt.shape, gives (1, 3).

there other problems well. in python 2.x, division of 2 integers returns integer. 5/9, example, 0. @ top of code should put from __future__ import division make behave way expect.

also, shouldn't using numpy matrices begin with. stick numpy arrays, better supported. way can avoid problem entirely using 1d array.

finally, safer import numpy np rather from numpy import *, can produce function name collisions built-in versions.