% 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. % Number of elements Ne = 10; x = linspace(-1,1,Ne+1); % 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); % Add forcing term to kn1 weighted residual b(kn1) = b(kn1) + COMPLETE THIS LINE; % Add forcing term to kn2 weighted residual b(kn2) = b(kn2) + COMPLETE THIS LINE; 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 plot(x,v,'*-'); xlabel('x'); ylabel('T'); % Plot exact solution xe = linspace(-1,1,2001); Te = -50*exp(xe) + 50*xe*sinh(1) + 100 + 50*cosh(1); hold on; plot(xe,Te); hold off;