3420.sagews

ttr1 = (0, 0); ttr2 = (0, 1); ttr3 = (1, 1); ttr = (ttr1, ttr2, ttr3)

Polynoms = PolynomialRing(QQ, 'x, y')
(x, y) = Polynoms.gens()

def S (P, cubatur, (tr1, tr2, tr3)):
    answer = 0;
    for (alpha1, alpha2, alpha3, w) in cubatur:
        answer = answer + (P(alpha1*tr1[0] + alpha2*tr2[0]+alpha3*tr3[0], alpha1*tr1[1] + alpha2*tr2[1]+alpha3*tr3[1]) + P(alpha1*tr1[0] + alpha3*tr2[0]+alpha2*tr3[0], alpha1*tr1[1]+alpha3*tr2[1]+alpha2*tr3[1]) + P(alpha2*tr1[0] + alpha1*tr2[0]+alpha3*tr3[0], alpha2*tr1[1]+alpha1*tr2[1]+alpha3*tr3[1]) + P(alpha3*tr1[0] + alpha1*tr2[0]+alpha2*tr3[0], alpha3*tr1[1]+alpha1*tr2[1]+alpha2*tr3[1]) + P(alpha2*tr1[0] + alpha3*tr2[0]+alpha1*tr3[0], alpha2*tr1[1]+alpha3*tr2[1]+alpha1*tr3[1]) + P(alpha3*tr1[0] + alpha2*tr2[0]+alpha1*tr3[0], alpha3*tr1[1]+alpha2*tr2[1]+alpha1*tr3[1]))*w
    return answer*(-(tr2[0]-tr1[0])*(tr3[1]-tr1[0])+(tr2[1]-tr1[1])*(tr3[0]-tr1[0]))/2

def integrr (P, p1, p2):
    pint = integr(P(p1[0] + x*(p2[0]-p1[0]), p1[1] + x*(p2[1]-p1[1]))*(p2[1]-p1[1]))
    if pint == 0:
        return 0
    else:
        return (pint(1,0) - pint(0,0))

def integr (P):
    if P == 0:
        return 0;
    if P.is_monomial():
        return P*x/(P.degree(x)+1)
    else:
        answer = 0
        for monom in P.monomials():
            answer = answer +  P.monomial_coefficient(monom)*monom*x/(monom.degree(x)+1)
        return answer

def myI (P, (tr1, tr2, tr3)):
    Q = integr(P)
    return integrr (Q, tr2, tr1) + integrr (Q, tr1, tr3) + integrr (Q, tr3, tr2)

def error (P, cubator, tr):
    return abs(S(P, cubator, tr)-myI(P,tr))

def split ((tr1, tr2, tr3), ord):
    tr11 = ((tr2[0]+tr3[0])/2, (tr2[1]+tr3[1])/2)
    tr21 = ((tr1[0]+tr3[0])/2, (tr1[1]+tr3[1])/2)
    tr31 = ((tr2[0]+tr1[0])/2, (tr2[1]+tr1[1])/2)
    if ord == 1:
        return ((tr1, tr31, tr21), (tr2, tr11, tr31), (tr3, tr21, tr11), (tr11, tr21, tr31))
    else:
        return split((tr1, tr31, tr21), ord-1) + split((tr2, tr11, tr31), ord-1) + split((tr3, tr21, tr11), ord-1) + split((tr11, tr21, tr31), ord-1)

eps=1e-10
def test (cubatur, ord):
    minorder = 100;
    maxmod = 0;
    mincord = (100, 0)
    for i in [0..ord]:
        for j in [0..ord]:
            locerr = error (x^i*y^j, cubatur, ttr)
            if (locerr > eps) & (i+j < minorder):
                minorder = i+j; mincord = (i, j); maxmod = locerr
            if (locerr > eps) & (i+j == minorder) & (locerr > maxmod):
                mincord = (i, j); maxmod = locerr
            print locerr, "&"
        print ("\\\\")
    print ("\n\n\n\n")
    P = x^mincord[0]*y^mincord[1]
    corrI = myI(P, ttr)
    print maxmod
    preveps = maxmod
    for i in [1..4]:
        curreps = abs(reduce (lambda u, v: u+v, map(lambda tr: S(P, cubatur, tr), split(ttr, i))) - corrI)
        print curreps/preveps, "&"
        preveps = curreps

test([(0.666666666666666, 0.166666666666667, 0.166666666666667, 0.166666666666667)], 3)

9.99200722162641e-16 & 6.66133814775094e-16 & 3.33066907387547e-16 & 0.000925925925925802 & \\ 3.33066907387547e-16 & 1.66533453693773e-16 & 0.000462962962962901 & 0.00192901234567899 & \\ 5.55111512312578e-17 & 0.000462962962962943 & 0.000385802469135846 & 0.00190145502645508 & \\ 0.000925925925925872 & 0.000385802469135742 & 0.000799162257495670 & 0.00219693072702340 & \\ 0.000925925925925872 3.37500000000014 & 6.40653935185195 & 0.711105810487330 & 1.24322491897020 &

test([(0.736712498969535, 0.237932366475534, 0.025355134554931, 0.104166666666667), (1/3, 1/3, 1/3, 1/16)], 5)

9.99200722162641e-16 & 7.21644966006352e-16 & 3.00925950824649e-13 & 6.60915766559356e-13 & 1.04083408558608e-12 & 0.000352733687494333 & \\ 3.60822483003176e-16 & 1.50462975412324e-13 & 3.30457883279678e-13 & 5.20417042793042e-13 & 0.000176366843747153 & 0.000594650206666575 & \\ 3.00579006129453e-13 & 2.70658495615805e-13 & 3.20646287299553e-13 & 0.0000352733690099569 & 0.000176954733011055 & 0.000438957476599000 & \\ 2.40751862889965e-13 & 2.20760909552808e-13 & 0.0000352733683586445 & 0.0000318930038167090 & 0.0000466392321854282 & 0.000200914495170918 & \\ 2.00776895109556e-13 & 0.000176366842846096 & 0.000316872427779918 & 0.000384087791258755 & 0.000372382258523395 & 0.000293228027552623 & \\ 0.000352733685892968 & 0.000639917695310840 & 0.000834019204218248 & 0.000933516232088955 & 0.000951582214241440 & 0.000906437026621859 & \\ 0.000352733687494333 338.554686129883 & 0.640941574362544 & 1.06482967125409 & 0.873514756600137 &