function U = MyHeatEqSolver(p, t, e, InitialData, EqData, BdryData, FinalTime) % Solves the heat equation cdu/dt-D(aDu)=f in Omega, with heat capacity c and % heat conductivity a, and with bdry conditions dudn = k(g-u) and initial % condition u = u0. % The data a, c and f are specified in EqData, the bdry data k and g are % specified in BdryData, and u0 is specified in InitialData (see the example % code below). The initial condition corresponds to time = 0 and the solution % is computed up to FinalTime. The geometry and spatial elements are specified % in p, t and e (as in MyPoissonSolver). The solution at time j*k is stored in % U(:,j+1). time = 0; u = feval(InitialData, p); if size(u,1)==1 u=u'; end U = u; k = 0.01; h = waitbar(0, 'Just a minute .. '); while time < FinalTime time = time + k/2; [S, M, b] = StiffMassLoad(p, t, e, EqData, BdryData, time); u = % to be completed!!!!!!!!!!!!! U=[U u]; time = time + k/2; waitbar(time/FinalTime,h) end % visualization figure(1) set(gcf,'renderer','zbuffer') h = plot(p,U(:,1)); set(gca,'XLim',[0 1],'YLim',[min(min(U)) max(max(U))]) for i=2:size(U,2) set(h,'YData',U(:,i)) drawnow end % subroutines function u0=InitialData(x) u0=zeros(size(x)); function [S, M, b] = StiffMassLoad(p, t, e, EqData, BdryData, time) % Returns the global stiffness matrix S = (s_ij), mass matrix M = (m_ij) % and load vector b = (b_j) with elements % s_ji = int_Omega a Dphi_i Dphi_j, m_ji = int_Omega c phi_i phi_j and % b_j = int_Omega f phi_j, respectively. % initializing n = size(p, 2); S = zeros(n, n); M = zeros(n, n); b = zeros(n, 1); S = sparse(S); M = sparse(M); % collecting element contributions to S, M and b for el = 1:size(t, 2) [dS, dM, db] = ElementContributions(p, t, EqData, time, el); nn = t(1:2, el); S(nn, nn) = S(nn, nn) + dS; M(nn, nn) = M(nn, nn) + dM; b(nn) = b(nn) + db; end % collecting boundary contributions to S, M and b for bel = 1:size(e, 2) [dS, db] = BoundaryContributions(p, e, BdryData, time, bel); nn = e(1, bel); S(nn, nn) = S(nn, nn) + dS; b(nn) = b(nn) + db; end function [dS, dM, db] = ElementContributions(p, t, EqData, time, el) % returns the element stiffness matrix dA and element load % vector db for element number el defined by t(1:3,el) and % p(t(1:2,el), and the equation data returned by EqData, % using element midpoint quadrature. NodeCoords = p(t(1:2, el)); x = sum(NodeCoords)/2; % element midpoint! [a, c, f] = EqData(x, time, t(3, el)); dx = diff(NodeCoords); % element length phi = [.5 .5]; Dphi = [-1/dx 1/dx]; dS = [a*Dphi(1)*Dphi(1)*dx a*Dphi(2)*Dphi(1)*dx;... a*Dphi(1)*Dphi(2)*dx a*Dphi(2)*Dphi(2)*dx]; dM = [c*phi(1)*phi(1)*dx c*phi(2)*phi(1)*dx;... c*phi(1)*phi(2)*dx c*phi(2)*phi(2)*dx]; db = [f*phi(1)*dx; f*phi(2)*dx]; function [dA, db] = BoundaryContributions(p, e, BdryData, time, bel) % returns the contributions dA and db from boundary element/node % bel to the stiffness matrix and load vector. NodeCoord = p(e(1, bel)); x = NodeCoord; [k, g] = BdryData(x, time, e(2, bel)); phi = 1; dA = k*phi*phi; db = k*g*phi; function [a, c, f] = EqData(x, t, tag) % returns the values at x and time t of the data a=a(x,t), c=c(x,t) % and f=f(x,t) in the differential equation cdu/dt-D(aDu)=f a = 1; c = 1; f = 0; function [k, g] = BdryData(x, t, tag) % returns the values at x of the data k=k(x,t) and g=g(x,t) if tag==2 k = 10000; g = sin(5*t); elseif tag==1 k = 10000; g = 0; end