1138.sagews

def f(x,y):
    return x*y^2
    
def c(x):
    return x^2
    
def d(x):
    return x^(.5)
    
def smpsn_dbl_int(a,b,n,m):
    h = (b-a)/m
    J = [0,0,0]
    for i in range(0,n,1):
        x=a+i*h
        HX = (d(x)-c(x))/m
        K = [  (f( x,c(x) ) + f( x,d(x) )),0,0]
        for j in range(1,m-1,1):
            y = c(x)+j*HX
            Q = f(x,y)
            if is_even(j):
                K[1] = K[1] + Q
            else:
                K[2] = K[2]+Q
        L = ((K[0]+2*K[1]+4*K[2])*HX)/3
        if i == 0:
            J[0] = J[0] + L
        elif i == n:
            J[0] = J[0] + L
        elif is_even(i):
            J[1] = J[1]+L
        else:
            J[2] = J[2]+L
    J = (h*(J[0]+2*J[1]+4*J[2]))/3
    return J

smpsn_dbl_int(0,1,14,14)

0.0459255061521566

def f(x,y):
    return exp(-(x+y))
    
def cee(x):
    return x^2
    
def dee(x):
    return x^(.5)
    
def gssn_dbl_int(a,b,m,n):
    h = [(b-a)/2,(b+a)/2]
    J=0
    r = [[1],[.5773502692,-.5773502692],[.7745966692,0,.7745966692],[.8611363116,.3399810436,-.3399810436,-.8611363116]]
    c = [[1],[1,1],[.55555555556,.88888888889,.55555555556],[.3478548451,.6521451549,.6521451549,.3478548451]]
    for i in range(1,m,1):
        JX = 0
        z = h[0]*r[m-1][i-1]+h[1]
        de = dee(z)
        ce = cee(z)
        k = [(de-ce)/2,(de+ce)/2]
        for j in range(1,n,1):
            v = k[0]*r[n-1][j-1]+k[1]
            Q = f(z,v)
            JX = JX + c[n-1][j-1]*Q
        J = J + c[m-1][i-1]*k[0]*JX
    J = h[0]*J
    return J

gssn_dbl_int(0,1,4,4)

0.0902286146371129