python - 我正在尝试将我的 matlab 代码转换为 python 以获得自适应高斯求积

Question

Matlab 代码运行良好

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 代码

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))

但是当我测试函数的 python 代码时，出现错误。出了什么问题？想要返回积分值，我更改了一些索引，但它一直返回积分值 0，这是我初始化积分的值。现在它告诉我错误是

IndexError Traceback(最近一次调用最后一次) 在 () 45返回积分 46---> 47 打印高斯(0,2,f,10**(-3)) 48

高斯分布(a、b、f、公差) 41如果j>n: 42 继续---> 43 a_j=abt[j,0] 44 b_j=abt[j,1] 45返回积分

C:\Users\Brandon Tran\AppData\Local\Enthought\Canopy\User\lib\site-packages\numpy\matrixlib\defmatrix.pyc in getitem(self, index) 314 第315章--> 316 out = N.ndarray.getitem(self, 索引) 第317章 最后: 第318章

索引错误:索引 1 超出尺寸为 1 的轴 0 的范围

Answer 1

无需自己实现。 quadpy (我创作的)具有自适应正交。

安装

pip install quadpy

并尝试

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

代码通过 Gauss-Kronrod 实现自适应性。