% FEM solver for d2T/dx2 + q = 0 where q = 50 exp(x) % % BC's: T(-1) = 100 and T(1) = 100. % % Note: the finite element degrees of freedom are % stored in the vector, v. % Set number of Gauss points (only used for forcing term in this example) NGf = 3; if (NGf == 3), xiGf = [-sqrt(3/5); 0; sqrt(3/5)]; aGf = [ 5/9; 8/9; 5/9]; elseif (NGf == 2), xiGf = [-1/sqrt(3); +1/sqrt(3)]; aGf = [ 1.0; 1.0]; else, NGf = 1; xiGf = [0.0]; aGf = [2.0]; end % Number of elements Ne = 5; x = linspace(-1,1,Ne+1); % Zero stiffness matrix K = zeros(2*Ne+1, 2*Ne+1); b = zeros(2*Ne+1, 1); % Loop over all elements and calculate stiffness and residuals for ii = 1:Ne, kn1 = 1 + 2*(ii-1); kn2 = 2 + 2*(ii-1); kn3 = 3 + 2*(ii-1); x1 = x(ii); x3 = x(ii+1); dx = x3 - x1; dxidx = 2/dx; dxdxi = 1/dxidx; % Add contribution to kn1 weighted residual K(kn1, kn1) = K(kn1, kn1) - dxidx*( 1/2); K(kn1, kn2) = K(kn1, kn2) - dxidx*( 0); K(kn1, kn3) = K(kn1, kn3) - dxidx*(-1/2); % Add contribution to kn2 weighted residual K(kn2, kn1) = K(kn2, kn1) - dxidx*( 0); K(kn2, kn2) = K(kn2, kn2) - dxidx*( 2/3); K(kn2, kn3) = K(kn2, kn3) - dxidx*( 0); % Add contribution to kn3 weighted residual K(kn3, kn1) = K(kn3, kn1) - dxidx*(-1/2); K(kn3, kn2) = K(kn3, kn2) - dxidx*( 0); K(kn3, kn3) = K(kn3, kn3) - dxidx*( 1/2); % Use Gaussian quadrature to evaluate forcing term integral for nn = 1:NGf, % Get xi for Gauss point xiG = xiGf(nn); % Find N1, N2 and N3 (i.e. weighting/intepolants) at xiG N1 = 0.5*(1-xiG); N2 = 0.5*(xiG^2 - 1); N3 = 0.5*(1+xiG); % Find x for Gauss point xG = 0.5*(1-xiG)*x1 + 0.5*(1+xiG)*x3; % Find f for Gauss point fG = -50*exp(xG); % Evaluate integrand at Gauss point for weight functions at nodes gG1 = N1*fG*dxdxi; gG2 = N2*fG*dxdxi; gG3 = N3*fG*dxdxi; % Send to correct right-hand side term b(kn1) = b(kn1) + aGf(nn)*gG1; b(kn2) = b(kn2) + aGf(nn)*gG2; b(kn3) = b(kn3) + aGf(nn)*gG3; end end % Set Dirichlet conditions at x=0 kn1 = 1; K(kn1,:) = zeros(size(1,2*Ne+1)); K(kn1, kn1) = 1.0; b(kn1) = 100.0; % Set Dirichlet conditions at x=1 kn1 = 2*Ne+1; K(kn1,:) = zeros(size(1,2*Ne+1)); K(kn1, kn1) = 1.0; b(kn1) = 100.0; % Solve for solution v = K\b; % Plot it and compare. Note: since even the finite element % solution varies more than linearly across an element, we need to % subdivide each element, evaluate the basis functions, and plot % the FEM solution to see the higher order variations. Nplot = 20; % Number of points per element to plot nnn = 0; for ii = 1:Ne, kn1 = 1 + 2*(ii-1); kn2 = 2 + 2*(ii-1); kn3 = 3 + 2*(ii-1); x1 = x(ii); x3 = x(ii+1); v1 = v(kn1); v2 = v(kn2); v3 = v(kn3); for nn = 1:Nplot, % Get xi for plot point xiG = -1 + 2*(nn-1)/(Nplot-1); % Find N1, N2 and N3 (i.e. weighting/intepolants) at xiG N1 = 0.5*(1-xiG); N2 = 0.5*(xiG^2 - 1); N3 = 0.5*(1+xiG); % Find x and v for plot point xG = 0.5*(1-xiG)*x1 + 0.5*(1+xiG)*x3; vG = v1*N1 + v2*N2 + v3*N3; nnn = nnn + 1; xp(nnn) = xG; vp(nnn) = vG; up(nnn) = -50*exp(xG) + 50*xG*sinh(1) + 100 + 50*cosh(1); end end figure(1); plot(xp,vp,'r');hold on; plot(xp,up); hold off; xlabel('x'); ylabel('u'); figure(2); plot(xp,vp-up); xlabel('x'); ylabel('Error');