MATH 5620 NUMERICAL ANALYSIS IIOPTIONAL HOMEWORK 6, DUE WED APRIL 21 2010The goal of this assignment is to use the finite element method with triangular piecewise linear (P1) finiteelements to solve the 2D BVP(1)(−∆u = f(x, y) in Ωu = 0 on ∂ΩHere the domain is the unit square Ω = [0, 1]2.1. Solve the problem with f(x, y) = 2(x(1 − x) + y(1 − y)) with true solution being utrue(x, y) = x(1 −x)y(1 − y), using P 1 triangular finite elements. Include plots with N ≡ Nx= Ny= 10 and N = 30 (seebelow for construction of the grid)2. Do a convergence study by plotting the error in L2(Ω) and H1(Ω) norms versus h in a log-log scale (youcan take for example N ∈ {10, 20, 30, 40, 50} and h = 1/N). Estimate the order of convergence. Youhave to convince me that your numerics agree with the theory i.e. that ku − uhkH1(Ω)= O(h) (first orderaccurate) and ku−uhkL2(Ω)= O(h2) (second order accurate). See below for how to estimate these norms.For comparison when N = 10 I get: ku − uhkL2(Ω)= 5.294373 · 10−4and ku − uhkH1(Ω)= 0.005981.1. Implementation details• The easiest way to obtain the grid and triangulation of Ω is as follows:x=lin space ( 0 , 1 ,Nx+1) ; % c r e a t e s d i s c r e t i z a t i o n in x d i r e c t i o ny=lin space ( 0 , 1 ,Ny+1) ; % c r e a t e s d i s c r e t i z a t i o n in y d i r e c t i o n[X,Y]=meshgrid ( x , y ) ; % c r e a t e s (Nx+1) ∗(Ny+1) e q u a l l y spac ed p o i n t spx=X( : ) ; py=Y( : ) ; % co n v e r t X and Y i n t o v e c t o r st r i = d elaunay ( px , py ) ; % do a c t u a l t r i a n g u l a t i o nNnodes = length ( px ) ; % number o f nodes ( d e g r e e s of freedom )Ne lt = si z e ( t r i , 1 ) ; % number o f e l e m e n t s ( t r i a n g l e s )• To assemble the stiffness matrix A and the right hand side F , see p106-108 inhttp://www.math.utah.edu/~fguevara/math5620_s10/na015.pdf.Do this without worrying about the boundary conditions (for now). Use sparse matrices other-wise solving the system will take a long time. You only need to initialize A=sparse(Nnodes,Nnodes);• In order to deal with the boundary conditions we first need to find all nodes that are on the boundary.The array ib below contains the indices of such nodes (of course this could be done directly, but thismethod would work even if the grid were unstructured)% f i n d th e nodes t h a t ar e on t h e boundaryi b=find ( abs ( px−0)<1e −10 | abs ( px−1)<1 e −10 . . .| abs ( py−0)<1 e −10 | abs ( py−1)<1e−10 ) ;• Now to enforce that U is zero at the boundary nodes we can hardwire the Dirichlet boundaryconditions in the linear system by modifying the matrix and right hand side as follows:fo r i =1: length ( i b ) , A( i b ( i ) , : ) = 0 ; A( i b ( i ) , i b ( i ) ) =1; end ;F( i b ) = 0 ;• Solve the system AU = F using Matlab’s backslash.• In order to plot your solution vector U you can use the following command: trimesh(tri,px,py,U).1• A good approximation of ku − uhkL2(Ω)is to look at the L2error between the interpolant of utrueon the grid and the computed uh. Thus we need a way of computing the L2norm of some functionvh∈ Vh(piecewise P 1) and then we will let vh= u − uh. If the elements cover Ω we have:kvhk2L2(Ω)=ZΩ(vh(x))2dx =XeZKe(vh(x))2dx.Then on each element we know that vh∈ P 1 therefore vh(x) =P3i=1vh(xtri(e,i))φei(x), where theφeiare the local basis functions. Thus the contribution of element e to the L2(Ω) norm is:ZKev2h(x)dx =ZKe 3Xi=1vh(xtri(e,i))φei(x)!2dx = vTlocMlocvloc,where Mlocis the same matrix given in p107 of the notes and vloc= [vh(xtri(e,1)), vh(xtri(e,2)), vh(xtri(e,3))]T.Do not forget to take square root after having summed the contributions of all elements.• Use a similar approximation for kvhkH1(Ω):kvhk2H1(Ω)=ZΩ(vh(x))2+ k∇vh(x)k2dx =XeZKe(vh(x))2+ k∇vh(x)k2dx.It is clear that the first term in the integrand has been computed already by the L2norm procedureabove. To compute the derivative term we do a similar procedure:|vh|2H1(Ke)=ZKek∇vh(x)k2dx = vTlocAlocvloc,where Alocis the same local stiffness matrix that is used in p106 of the notes to construct the stiffnessmatrix. Again: do not forget to take square root after having summed the contributions ofall elements. This shows that FEM gives solutions that are meaningful even after differentiation (ofcourse this depends on the polynomial space that is