function [u, A, b] = MyPoissonSolver(p, t, e, EqData, BdryData) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% Purpose: Solves D(-aDu)+cu=f in a domain D with %% boundary conditions anDu=k(g-u) where nDu %% nDu denotes the derivative out of D. %% %% p 1 x n vector of node coordinates. %% %% t 3 x m matrix of element node numbers t(1,el)=nn1 %% and t(2,el)=nn2, 1 <= el <= m, pointing to %% corresponding node coordinates p(nn1) and p(nn2), %% and subdomain reference tags t(3,el), useful %% if data are given by different expresssions in %% different parts of the domain. %% %% e 2 x q matrix of boundary node numbers and %% reference tags, with node number e(1,bel)=nn pointing %% to the boundary node coordinate p(nn), and e(2,bel) %% is a reference tag (see example below). %% %% EqData matlab function returning [a,c,f] given %% the coordinate x and subdomain reference tag. %% %% BdryData matlab function returning [k,g] given x and %% a boundary reference tag. %% %% A n x n global stiffness matrix %% %% b n x 1 global load vector %% %% u n x 1 solution at node points p %% %% %% Example: To solve D(-(1+x)Du)=1/x in [0,1] with boundary %% conditions Du(0)=u(0) and u(1)=7 (k = oo), define %% p = 0:0.01:1, t = [1:100; 2:101; ones(1, 100)], %% e = [1 101; 1 2] and functions %% %% function [a, c, f] = EqData(x, tag) %% a = 1+x; c = 0; f = 1/x; %% %% and %% %% function [k, g] = BdryData(x, tag) %% if tag==1 k = 1; g = 0; end %% if tag==2 k = 10000; g = 7; end %% %% To call the solver, type %% [u, A, b] = MyPoissonSolver(p, t, e, 'EqData', 'BdryData') %% %% To plot the solution, type plot(p, u) %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % initializing n = size(p, 2); A = zeros(n, n); b = zeros(n, 1); % collecting element contributions to A and b for el = 1:size(t, 2) [dA, db] = ElementContributions(p, t, EqData, el); nn = t(1:2, el); A(nn, nn) = A(nn, nn) + dA; b(nn) = b(nn) + db; end % collecting boundary contributions to A and b for bel = 1:size(e, 2) [dA, db] = BoundaryContributions(p, e, BdryData, bel); nn = e(1, bel); A(nn, nn) = A(nn, nn) + dA; b(nn) = b(nn) + db; end % solving resulting system u = A\b; % subroutines function [dA, db] = ElementContributions(p, t, EqData, 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, t(3, el)); dx = phi = Dphi = dA = % the 2x2 element stiffness matrix % int_el a Dphi(I) Dphi(J) + c phi(I) phi(J) dx % I,J = 1,2 db = % the 2x1 element load vector % int_el f phi(J) dx, J=1,2 function [dA, db] = BoundaryContributions(p, e, BdryData, 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, e(2, bel)); ds = 1; phi = 1; dA = k*phi*phi*ds; % the 1x1 bdry element contribution % int_bel k phi(I) phi(J) ds, I,J = 1,2 db = k*g*phi*ds; % the 1x1 bdry element contribution % int_bel k g phi(J) ds, J=1 function [a, c, f] = EqData(x, tag) % returns the values at x of the data a=a(x), c=c(x) and f=f(x) in the % differential equation D(-aDu)+cu=f a = 1; c = 0; f = 1/x; function [k, g] = BdryData(x, tag) % returns the values at x of the data k=k(x) and g=g(x) if tag==1 k = 10000; g = 0; elseif tag==2 k = 1; g = 1; end