1137.sagews

def f(x,y):
    return x*y^2
    
def c(x):
    return 1.2
    
def d(x):
    return 1.4
    
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

%latex
This is the 4.4 Algorithm answer
\sage{smpsn_dbl_int(2.1,2.5,4,4)}

var('x')
var('y')
actual = integrate(integrate(x*y^2,y,1.2,1.4),x,2.1,2.4)

def prob2():
    var('x')
    var('y')
    actual = integrate(integrate(x*y^2,y,1.2,1.4),x,2.1,2.4)
    H = 1
    i = 1
    answer = []
    while H == 1:
        error  = abs(smpsn_dbl_int(2.1,2.5,i,i)-actual)
        if error <=10^(-6):
            H = 0
            print answer
            print 'n=m= ', i
            break
        else:
            answer.append(error)
            i = i+1

answer = []
answer.append(3)