3838.sagews

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

defaultcubatur = [(0.736712498969535, 0.237932366475534, 0.025355134554931, 0.104166666666667), (0.333333333333333, 0.333333333333333, 0.333333333333333, 0.0625)]

def S0 (P, (tr1, tr2, tr3), f=lambda x,y:x , cubatur=defaultcubatur):
    answer = 0;
    for (alpha1, alpha2, alpha3, w) in cubatur:
       answer = answer + (        return c0
    f(P(alpha1*tr1[0]+alpha2*tr2[0]+alpha3*tr3[0], alpha1*tr1[1]+alpha2*tr2[1]+alpha3*tr3[1]), (alpha1, alpha2, alpha3)) +
    f(P(alpha1*tr1[0]+alpha3*tr2[0]+alpha2*tr3[0], alpha1*tr1[1]+alpha3*tr2[1]+alpha2*tr3[1]), (alpha1, alpha3, alpha2)) + 
    f(P(alpha2*tr1[0]+alpha1*tr2[0]+alpha3*tr3[0], alpha2*tr1[1]+alpha1*tr2[1]+alpha3*tr3[1]), (alpha2, alpha1, alpha3)) +
    f(P(alpha3*tr1[0]+alpha1*tr2[0]+alpha2*tr3[0], alpha3*tr1[1]+alpha1*tr2[1]+alpha2*tr3[1]), (alpha3, alpha1, alpha2)) +
    f(P(alpha2*tr1[0]+alpha3*tr2[0]+alpha1*tr3[0], alpha2*tr1[1]+alpha3*tr2[1]+alpha1*tr3[1]), (alpha2, alpha3, alpha1)) +
    f(P(alpha3*tr1[0]+alpha2*tr2[0]+alpha1*tr3[0], alpha3*tr1[1]+alpha2*tr2[1]+alpha1*tr3[1]), (alpha3, alpha2, alpha1)))*w
    return answer

def isvalidpoint ((i,j), size):
    return (i>=0)&(i<size)&(j>=0)&(j<=i)
    
def weight ((i0,j0), (i1,j1), size): 
    if (i0,j0) == (i1,j1):
        return len(neartr((i0,j0),size))*2
    if (j0==0)&(j1==0)|((i0-j0)==0)&((i1-j1)==0)|(i0==size-1)&(i1==size-1):
        return 1
    return 2

        
def isvalidtr (tr, size):
    return reduce (lambda b, point: b & isvalidpoint (point, size), tr, True)

def neartr ((i,j), size):
    return filter (lambda tr: isvalidtr(tr, size), [((i,j), (i,j-1), (i-1,j-1)), ((i,j), (i-1,j-1), (i-1,j)), ((i,j), (i-1,j), (i,j+1)), ((i,j), (i,j+1), (i+1,j+1)), ((i,j), (i+1, j+1), (i+1,j)), ((i,j), (i+1,j), (i,j-1))])

def weights ((i,j),size):
    return map(lambda (i0,j0): (weight((i,j),(i0,j0), size)*1/12, pointtonom((i0,j0))),
        filter (lambda point: isvalidpoint(point, size), [(i,j), (i,j-1), (i-1,j-1), (i-1,j), (i,j+1), (i+1, j+1), (i+1,j)]))


def pointtonom ((i,j)):
    return i*(i+1)/2+j

def pointtocord ((i,j), size):
    return (ttr1[0]+i*(ttr2[0]-ttr1[0])/(size-1)+ j*(ttr3[0]-ttr2[0])/(size-1),
            ttr1[1]+i*(ttr2[1]-ttr1[1])/(size-1)+ j*(ttr3[1]-ttr2[1])/(size-1))
            
def trtocord (tr, size):
    return map (lambda point:pointtocord(point,size), tr)

def multbyg (f, point, size):
    return reduce(lambda x,y: x+y, 
                  map(lambda tr: S0(f, tr, lambda x,(alpha1,alpha2,alpha3):x*alpha1), 
                      map(lambda tr:trtocord(tr,size),
                          neartr(point,size))))

def rightside (f, size):
    result=[]
    for i0 in [0..(size-1)]:
        for j0 in [0..i0]:
            result.append((multbyg (f, (i0,j0), size), i0, j0))
    return result

def zero(size, a=0):
    result = []
    for i0 in [0..(size-1)]:
        for j0 in [0..i0]:
            result.append((a, i0, j0))
    return result

def multbyA (x, size):
#    print 'start'
    result = []
    for i0 in [0..(size-1)]:
        for j0 in [0..i0]:
            result.append((reduce (lambda tmp, (a,n): tmp+a*x[n][0], weights((i0,j0),size), 0), i0, j0))
#            result.append((weights((i0,j0),size), i0, j0))
#    print 'end'
    return result

def l2prod (l1, l2):
    return reduce (lambda x, (y1,y2): x+y1[0]*y2[0], zip(l1,l2), 0)
    
def l2norm (l1):
    return numerical_approx(sqrt(reduce(lambda x, y: x+y[0]*y[0], l1, 0)))

def zipwith (l1,l2,f):
    return map(f, zip(l1,l2))

def mult (l1, alpha):
    return map(lambda (x,i,j): (alpha*x,i,j),l1)
    
def minus (l1, l2):
    return zipwith (l1,l2, lambda ((x1,i1,j1),(x2,i2,j2)):((x1-x2),i1,j1))

def solveBiCGStab (size, F, eps):
    c0 = zero(size)
    r0 = minus(F, multbyA(c0, size))
    p0 = r0
    z  = r0
    beta = l2norm (r0)
    if (beta < eps):
        return c0
#    j=0  
    while (l2norm(r0)/beta > eps):
#        j=j+1
#        print 'c0'
#        print c0
        alpha0 = l2prod(r0,z)/l2prod(multbyA(p0, size), z)
#        print 'r0'
#        print r0
        s0 = minus(r0, mult(multbyA(p0, size), alpha0))
#        print 'l2norm'
#        print l2prod(multbyA(s0, size), multbyA(s0, size))
        w0 = l2prod(multbyA(s0, size),s0)/l2prod(multbyA(s0, size), multbyA(s0, size))
        c1 = minus(minus(c0, mult(p0, -alpha0)), mult(s0, -w0))
        r1 = minus(s0, mult(multbyA(s0, size), w0))
        beta0 = numerical_approx(l2prod(r1,z)/l2prod(r0,z)*alpha0/w0)
        p1 = minus(r1,mult(minus(p0,mult(multbyA(p0,size),w0)),-beta0))
        p0 = p1
        c0 = c1
        r0 = r1
#    print j
    return c0

def errortr (f, tr, resh, size):
    return S0(f, (pointtocord(tr[0], size), pointtocord(tr[1], size), pointtocord(tr[2], size)), lambda x,(alpha1,alpha2,alpha3):(x - alpha1*resh[pointtonom(tr[0])][0] - alpha2*resh[pointtonom(tr[1])][0] - alpha3*resh[pointtonom(tr[2])][0])^2)

def split (size):
    split = []
    for i in [1..size-1]:
        for j in [0..i-1]:
            split.append (((i,j),(i-1,j),(i,j+1)))
    for i in [1..size-2]:
        for j in [0..i-1]:
            split.append (((i,j),(i,j+1),(i+1,j)))
    return split

def error (f, answ, size):
    return reduce (lambda a,b:a+b,
                   map(lambda tr: errortr(f,tr, answ, size), split(size)))

def test (f, size, eps):
    return error (f, solveBiCGStab(size, rightside (f, size), eps), size)

def tests (f, eps):
    t1 = test (f, 9, eps)
    print t1
    t2 = test (f, 17, eps)
    print t2
    t3 = test (f, 33, eps)
    print t3
    t4 = test (f, 65, eps)
    return (t1, t2, t3, t4)

tests (lambda x,y:1, 1e-6)

7.17985244868109e-10 9.65101085758785e-11 3.06473588508056e-9 (7.17985244868109e-10, 9.65101085758785e-11, 3.06473588508056e-9, 7.60323318058235e-9)

tests (lambda x,y:x+y, 1e-6)

2.56698576860977e-10 2.02940885577074e-10 2.07636958156847e-9 (2.56698576860977e-10, 2.02940885577074e-10, 2.07636958156847e-9, 5.27501943114878e-9)

tests (lambda x,y: y*y, 1e-6)

0.0000488034800092498 0.0000122057584881344 3.05195630030148e-6 (0.0000488034800092498, 0.0000122057584881344, 3.05195630030148e-6, 7.63460414187703e-7)