function blatt7_1D pmax = 4; mmax = 7; %create convergence plots for different values of p and uniform meshes exact_energy = pi; energy = zeros(mmax,pmax); h = zeros(mmax,pmax); % for p=1:pmax for m=1:mmax M = 2^m; energy(m,p) = fem_1D(M,p) h(m,p) = 1.0/M; h end; end; error = sqrt(abs(exact_energy - energy)); figure(1) loglog(h(:,1),error(:,1),'*-k',... h(:,2),error(:,2),'o--b',... h(:,3),error(:,3),'d:r',... h(:,4),error(:,4),'+-.m','LineWidth',2) set(gca,'FontSize',14) axis([1/200 1 10^(-7) 1]) xlabel('mesh size h') ylabel('energy norm error') title('-u^{\prime\prime} + u =f; solution u(x) = sin(x); uniform meshes') legend('p=1','p=2','p=3','p=4','Location','SouthEast') figure(2) semilogy([1:pmax],,error(1,:),'*-k',LineWidth',2) set(gca,'FontSize',14) xlabel('polynomial degree p') ylabel('energy norm error') title('-u^{\prime\prime} + u =f; solution u(x) = sin(x); uniform meshes') return end %------------------------------------------------------------------------------- function energy = fem_1D(m,p) %m = parameter for the mesh (here: uniform mesh) %p = polynomial degree % %numbering of the unknowns is such that the first and the last %unknown correspond to the left and right endpts % g = -1; %value of the Neumann bc at the right endpt %define the mesh xleft = 0; xright = pi; h = (xright-xleft)/m; mesh =[xleft:h:xright]; %compute the stiffness matrix and the load vector without boundary conditions [B,l] = assemble_neumann_problem(mesh,p); N = size(B,1); %get size of B % %ENFORCE NEUMANN BOUNDARY CONDITIONS AT RIGHT ENDPOINT %COMPUTE FEM SOLUTION %COMPUTE ENERGY % % return end %------------------------------------------------------------------------------- function [B,l] = assemble_neumann_problem(mesh,p); %input: a mesh, i.e., a (sorted) list of nodes x1, x2,...,xM %input: polynomial degree p % % [T,N] = local_to_global_maps(mesh,p); % B = zeros(N,N); l = zeros(N,1); % %LOOP OVER ALL ELEMENTS AND ASSEMBLE % return end %------------------------------------------------------------------------------- function [T,N] = local_to_global_maps(mesh,p) %input: the mesh and the polynomial degree %output: matrix T, where T(K,i) gives the number of the global DOF corresponding % to the function N_i on element K % N = total number of DOF % recall: N_1, N_2 are linears, N_i for i \ge 3 are bubbles % global numbering: for each element, first the left endpt, then all bubbles, then % the right endpt % %SET UP THE MATRIX T AND COMPUTE N return end %------------------------------------------------------------------------------- function B = element_stiffness_matrix(vertex_coords,p) % vertex_coords(1:2) contains the two endpoints of the element % p >= 1 polynomial degree to be employed % h = vertex_coords(2) - vertex_coords(1); K = zeros(p+1,p+1); M = zeros(p+1,p+1); q = p+1; %number of Gaussian quadrature points to be employed [x,w] = gauleg(q); %x = quadrature points on (-1,1); w = quadrature weights x = x'; w = w'; %gauleg returns column vectors, we want row vectors PHI = N(x,p+1); %evaluate the shape fcts N_i at the quadrature points GRAD_PHI = grad_N(x,p+1); %evaluate N_i' at the quadrature points % %COMPUTE THE ELEMENT STIFFNESS MATRIX % return end %------------------------------------------------------------------------------- function l = element_load_vector(vertex_coords,p) % vertex_coords(1:2) contains the two endpoints of the element % p >= 1 polynomial degree to be employed % h = vertex_coords(2) - vertex_coords(1); l = zeros(p+1,1); q = p+1; %number of Gaussian quadrature points to be employed [x,w] = gauleg(q); %x = quadrature points on (-1,1); w = quadrature weights x = x'; w = w'; %gauleg returns column vectors, we want row vectors PHI = N(x,p+1); %evaluate the shape fcts N_i at the quadrature points % %COMPUTE THE ELEMENT LOAD VECTOR (WITHOUT THE NEUMANN BC) % return end %------------------------------------------------------------------------------- function [U] =N(x,n) % % Calculates the first n trial functions on [-1,1]; % we assume implicitly that the vector x is a row vector % % N_1 = (1-x)/2 % N_2 = (1+x)/2 % N_n = \sqrt{(2n-3)/2}\int_{-1}^x L_{n-2}(t) dt if n \ge 3 % U(1,:) = (1-x)/2; U(2,:) = (1+x)/2; if n >2 [P] = legtable(x,n); for i = 3:n ip=i-2; U(i,:) = sqrt((2*i-3)/2)*1/(2*ip+1)*(P(ip+2,:)-P(ip,:)); end; end return end %------------------------------------------------------------------------------- function [U] =grad_N(x,n) % % Calculates the derivatives of the first n trial functions on [-1,1]; % we assume implicitly that the vector x is a row vector % % N_1 = (1-x)/2 % N_2 = (1+x)/2 % N_n = \sqrt{(2n-3)/2}\int_{-1}^x L_{n-2}(t) dt if n \ge 3 % U(1,:) = -1/2*ones(size(x)); U(2,:) = 1/2*ones(size(x)); if n >2 [P] = legtable(x,n); for i = 3:n ip=i-2; U(i,:) = sqrt((2*i-3)/2)*P(ip+1,:); end; end return end %------------------------------------------------------------------------------- function [P] = legtable(x,m) % % input: row vector x \subset (-1,1) of points % m the maximal order of the Legendre polynomials % output: a matrix P=P(m+1,length(x)) containing the values of the % Legendre polynomials at the points x % l=length(x); P=ones(m+1,l); P(2,:)=x; for i=2:m P(i+1,:)=((2*i-1)*x.*P(i,:)-(i-1)*P(i-1,:))/i; end; return end %------------------------------------------------------------------------------- function [x,w]=gauleg(n) x=zeros(n,1); w=zeros(n,1); m=(n+1)/2; xm=0.0; xl=1.0; for i=1:m z=cos(pi*(i-0.25)/(n+0.5)); while 1 p1=1.0; p2=0.0; for j=1:n p3=p2; p2=p1; p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j; end pp=n*(z*p1-p2)/(z*z-1.0); z1=z; z=z1-p1/pp; if (abs(z-z1)