1_5

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)

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)

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)

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