#Nathan Fox #Homework 25 #I give permission for this work to be posted online #Read procedures from class read(`C25.txt`): Help:=proc(): print(` Tgh(N,h) , FEMh(p,q,r,f,g,x,y,N,h) , CheckFEM(u,p,q,r,x,y,pt,H:=[1/4,1/8/,1/12,1/16]) `): end: ##PROBLEM 2## #Tgh(N,h): triangulates into right-angled triangles #of edge-side h the square [0,N]x[0,N] #h is the reciprocal of a positive integer Tgh:=proc(N, h) local i,j: #{seq(seq({[i,j],[i+1,j],[i,j+1],[i+1,j+1]},i=0..N-1),j=0..N-1)}: return {seq(seq({[i*h,j*h],[(i+1)*h,j*h],[i*h,(j+1)*h]},i=0..N/h-1),j=0..N/h-1)} union {seq(seq({[(i+1)*h,j*h],[i*h,(j+1)*h],[(i+1)*h,(j+1)*h]},i=0..N/h-1),j=0..N/h-1)}: end: ##PROBLEM 3## #FEMh(p,q,r,f,g,x,y,N,h): the finite element approximate solution #to the elliptic PDE Dirichlet BVP #(12.41), u=g on the boundary using the #triangluation Tgh(N,h) FEMh:=proc(p, q, r, f, g, x, y, N, h) local w, c, var, B, T, i, b, N1244, eq, v1: T:=Tgh(N,h): B:={seq([0, i*h], i=0..N/h), seq([i*h, 0], i=0..N/h), seq([N, i*h], i=0..N/h), seq([i*h, N], i=0..N/h)}: w, var:=TestF(T,c,x,y): for b in B do w:=subs(c[b]=subs({x=b[1], y=b[2]}, g), w): var:=var minus {c[b]}: od: N1244:=Iw(w,p,q,r,f,x,y): eq:={seq(diff(N1244, v1), v1 in var)}: var:=solve(eq, var): if var = NULL then return FAIL: fi: return subs(var, w): end: ##PROBLEM 4## #CheckFEM(u,p,q,r,x,y,pt,H:=[1/4,1/8,1/12,1/16]): #starts with a solution, u #(an expression of x,y), expressions p,q,r as in (12.41) of the #book, computes f, and then computes #FEMh(p,q,r,f,u,x,y,1,h) #(note that g is u, since we already know the solution and it #better be equal to g on the boundary) at #h as each element of H #and plugs it in at the point p. The output is a list of nops(H) #numbers, followed by the EXACT value of u at (x,y)=pt. #NOTE: H defaults to [1/4,1/8,1/12,1/16] CheckFEM:=proc(u, p, q, r, x, y, pt, H:=[1/4,1/8,1/12,1/16]) local h, ret: ret:=[]: for h in H do ret:=[op(ret), EvalF1(FEMh(p,q,r,f,u,x,y,1,h), x, y, pt)]: od: return ret, subs({x=pt[1], y=pt[2]}, u): end: #CheckFEM(x+3*y+10,x+y+3,-3-x^2*y,x^2+y,x,y,[0.32,0.52]); returned [11.88, 11.88, 11.88, 11.88], 11.88 #CheckFEM(x^2*y+7*x*y+2,x^2+y+5,-4*x-y,x^3+3*y,x,y,[0.32,0.52]); rereturned [2, 2, 2, 2], 3.218048