% FEM solver for d2T/dx2 + q = 0 where q = 50 exp(x) % % BC's: T(-1) = 100 and T(1) = 100. % % Gaussian quadrature is used in evaluating the forcing integral. % % Note: the finite element degrees of freedom are % stored in the vector, v. % Number of elements Ne = 10; x = linspace(-1,1,Ne+1); % Set quadrature rule Nq = 2; if (Nq == 1), alphaq(1) = 2.0; xiq(1) = 0.0; elseif (Nq == 2), alphaq(1) = 1.0; xiq(1) = -1/sqrt(3); alphaq(2) = 1.0; xiq(2) = 1/sqrt(3); else fprintf('Error: Unknown quadrature rule (Nq = %i)\n',Nq); return; end % Zero stiffness matrix K = zeros(Ne+1, Ne+1); b = zeros(Ne+1, 1); % Loop over all elements and calculate stiffness and residuals for ii = 1:Ne, kn1 = ii; kn2 = ii+1; x1 = x(kn1); x2 = x(kn2); dx = x2 - x1; % Add contribution to kn1 weighted residual due to kn1 function K(kn1, kn1) = K(kn1, kn1) - (1/dx); % Add contribution to kn1 weighted residual due to kn2 function K(kn1, kn2) = K(kn1, kn2) + (1/dx); % Add contribution to kn2 weighted residual due to kn1 function K(kn2, kn1) = K(kn2, kn1) + (1/dx); % Add contribution to kn2 weighted residual due to kn2 function K(kn2, kn2) = K(kn2, kn2) - (1/dx); % Evaluate forcing term using quadrature for nn = 1:Nq, % Get xi location of quadrature point xi = xiq(nn); % Calculate x location of quadrature point xq = x1 + 0.5*(1+xi)*dx; % Calculate q qq = 50*exp(xq); % Calculate phi1 and phi2 phi1 = 0.5*(1-xi); phi2 = 0.5*(1+xi); % Subtract forcing term to kn1 weighted residual b(kn1) = b(kn1) - alphaq(nn)*0.5*phi1*qq*dx; % Add forcing term to kn2 weighted residual b(kn2) = b(kn2) - alphaq(nn)*0.5*phi2*qq*dx; end end % Set Dirichlet conditions at x=0 kn1 = 1; K(kn1,:) = zeros(size(1,Ne+1)); K(kn1, kn1) = 1.0; b(kn1) = 100.0; % Set Dirichlet conditions at x=1 kn1 = Ne+1; K(kn1,:) = zeros(size(1,Ne+1)); K(kn1, kn1) = 1.0; b(kn1) = 100.0; % Solve for solution v = K\b; % Plot solution figure(1); plot(x,v,'*-'); xlabel('x'); ylabel('T'); % For the exact solution, we need to use finer spacing to plot % it correctly. If we only plot it at the nodes of the FEM mesh, % the exact solution would also look linear between the nodes. To % make sure there is always enough resolution relative to the FEM % nodes, the size of the vector for plotting the exact solution is % set to be 20 times the number of FEM nodes. Npt = 20*Ne+1; xe = linspace(-1,1,Npt); Te = -50*exp(xe) + 50*xe*sinh(1) + 100 + 50*cosh(1); hold on; plot(xe,Te); hold off; % Plot the error. To do this, calculate the error on the same % set of points in which the exact solution was plot. This % requires that the location of the point xx(i) be found in the % FEM mesh to construct the true solution at this point by linearly % interpolating between the two nodes of the FEM mesh. verr(1) = v(1) - Te(1); h = x(2)-x(1); for i = 2:Npt-1, xxi = xe(i); Ti = Te(i); j = floor((xxi-xe(1))/h) + 1; x0 = x(j); x1 = x(j+1); v0 = v(j); v1 = v(j+1); xi = 2*(xxi - x0)/(x1-x0)-1; % This gives xi between +/-1 vi = 0.5*(1-xi)*v0 + 0.5*(1+xi)*v1; verr(i) = vi - Ti; end verr(Npt) = v(Ne+1) - Te(Npt); figure(2); plot(xe,verr); xlabel('x'); ylabel('Error');