Source code (MATLAB)

% Solve the problem -\Delta u + c*u = f with Dirichlet boundary conditions.% The domain is the disk centered at the origin of radius 1.% Its triangulation is read from the end of this file.functionmain()c=0;% a parameter in the equation, see abovewhite=0.99*[1,1,1];blue=[0,129,205]/256;green=[0,200,70]/256;% load the triangulation from the end of the file.dummy_arg=0;P=get_points(dummy_arg);T=get_triangles(dummy_arg);% find the number of points and the number of triangles[Np,k]=size(P);[Nt,k]=size(T);Nb=30;% first Nb points in P are on the boundary% plot the triangulationlw=1.4;figure(1);clf;holdon;axisequal;axisoff;forl=1:Nti=T(l,1);j=T(l,2);k=T(l,3);% plot the three edges in each triangleplot([P(i,1),P(j,1)],[P(i,2),P(j,2)],'linewidth',lw,'color',blue)plot([P(i,1),P(k,1)],[P(i,2),P(k,2)],'linewidth',lw,'color',blue)plot([P(j,1),P(k,1)],[P(j,2),P(k,2)],'linewidth',lw,'color',blue)end% a hack to deal with bounding box issuess=1.05;plot(-s,-s,'*','color',white);plot(s,s,'*','color',white);% save as eps and svg (needs the plot2svg function)saveas(gcf,'triangulation.eps','psc2')%  plot2svg('Finite_element_triangulation.svg');% right-hand side%  f=inline ('5-x.^2-y.^2', 'x', 'y');f=inline('4+0*x.^2+0*y.^2','x','y');% the values of f at the nodes in the triangulationF=f(P(:,1),P(:,2));RHS=0*F;RHS=RHS((Nb+1):Np);% will solve A*U=RHS% an empty sparse matrix of size Np by NpA=sparse(Np-Nb,Np-Nb);% iterate through trianglesforl=1:Nt%     fill in the matrixi=T(l,1);j=T(l,2);k=T(l,3);[pii,pij,pik,pjj,pjk,pkk,gii,gij,gik,gjj,gjk,gkk]=calc_elems(i,j,k,P(i,:),P(j,:),P(k,:));%     One has to consider the cases when a given vertex is on the boundary, or not.%     If yes, it can't be included in the matrixifi>NbA(i-Nb,i-Nb)=A(i-Nb,i-Nb)+gii+c*pii;endifi>Nb&j>NbA(i-Nb,j-Nb)=A(i-Nb,j-Nb)+gij+c*pij;A(j-Nb,i-Nb)=A(i-Nb,j-Nb);endifi>Nb&k>NbA(i-Nb,k-Nb)=A(i-Nb,k-Nb)+gik+c*pik;A(k-Nb,i-Nb)=A(i-Nb,k-Nb);endifj>NbA(j-Nb,j-Nb)=A(j-Nb,j-Nb)+gjj+c*pjj;endifj>Nb&k>NbA(j-Nb,k-Nb)=A(j-Nb,k-Nb)+gjk+c*pjk;A(k-Nb,j-Nb)=A(j-Nb,k-Nb);endifk>NbA(k-Nb,k-Nb)=A(k-Nb,k-Nb)+gkk+c*pkk;end%     add the appropriate contributions to the right-hand termsifi>NbRHS(i-Nb)=RHS(i-Nb)+F(i)*pii+F(j)*pij+F(k)*pik;endifj>NbRHS(j-Nb)=RHS(j-Nb)+F(i)*pij+F(j)*pjj+F(k)*pjk;endifk>NbRHS(k-Nb)=RHS(k-Nb)+F(i)*pik+F(j)*pjk+F(k)*pkk;endend%  plot the sparse matrix (more exactly, its sign, then it is easier to see the pattern)figure(6);imagesc(sign(abs(A)));colormap(1-gray);axisij;axisequal;axisoff;%  save as eps and svg (needs the plot2svg function)saveas(gcf,'sparse_dirichlet.eps','psc2')%  plot2svg('Finite_element_sparse_matrix.svg');%  calculate U, then add zeros for the points on the boundaryU_calc=A\RHS;U_calc=[0*(1:Nb)U_calc']';%  exact solutionu=inline('1-x.^2-y.^2','x','y');U_exact=u(P(:,1),P(:,2));%  plot the computed solutionfigure(2);clf;plot_solution(T,P,U_calc,lw,green)%  save as eps and svg (needs the plot2svg function)saveas(gcf,'computed_dirichlet.eps','psc2')%  plot2svg('Finite_element_solution.svg');%   plot the exact solution%   figure(3); clf;%   plot_solution(T, P, U_exact)%   title('Exact Dirichlet solution')%   saveas(gcf, 'exact_dirichlet.eps', 'psc2')disp(sprintf('Error between exact solution and computed solution is %0.9g',max(abs(U_exact-U_calc))))% given the l-th triangle, calculate the integrals% pij = int_K(lambda_i*lambda_j) and  gij = int_K(grad lambda_i dot grad lambda_j)function[pii, pij, pik, pjj, pjk, pkk, gii, gij, gik, gjj, gjk, gkk] = calc_elems(i, j, k, P1, P2, P3)%  extract the x and y coordinates of the pointsx1=P1(1);y1=P1(2);x2=P2(1);y2=P2(2);x3=P3(1);y3=P3(2);%  triangle areaAT=abs((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1))/2;%  the integrals of lambda_i*lambda_jpii=AT/6;pij=AT/12;pik=AT/12;pjj=AT/6;pjk=AT/12;pkk=AT/6;%  the integrals of grad lambda_i dot grad lambda_jgii=((x2-x3)*(x2-x3)+(y2-y3)*(y2-y3))/(4*AT);gjj=((x1-x3)*(x1-x3)+(y1-y3)*(y1-y3))/(4*AT);gkk=((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))/(4*AT);gij=-((x3-x1)*(x3-x2)+(y3-y1)*(y3-y2))/(4*AT);gik=-((x2-x1)*(x2-x3)+(y2-y1)*(y2-y3))/(4*AT);gjk=-((x1-x2)*(x1-x3)+(y1-y2)*(y1-y3))/(4*AT);functionplot_solution(T, P, U, lw, color, fs)[Np,k]=size(P);[Nt,k]=size(T);% create a path in 3D that will trace the outline of the solutionX=zeros(5*Nt,1);Y=zeros(5*Nt,1);Z=zeros(5*Nt,1);forl=1:Nti=T(l,1);j=T(l,2);k=T(l,3);m=5*l-4;X(m)=P(i,1);Y(m)=P(i,2);Z(m)=U(i);X(m+1)=P(j,1);Y(m+1)=P(j,2);Z(m+1)=U(j);X(m+2)=P(k,1);Y(m+2)=P(k,2);Z(m+2)=U(k);X(m+3)=P(i,1);Y(m+3)=P(i,2);Z(m+3)=U(i);X(m+4)=NaN;Y(m+4)=NaN;Z(m+4)=NaN;end%  plot the solutionplot3(X,Y,Z,'linewidth',lw,'color',color)set(gca,'fontsize',15)functionP=get_points(dummy_arg)P=[100.977999999999999980.207999999999999990.914000000000000030.406999999999999970.809000000000000050.587999999999999970.669000000000000040.742999999999999990.50.865999999999999990.3090.950999999999999960.1050.995-0.1050.995-0.3090.95099999999999996-0.50.86599999999999999-0.669000000000000040.74299999999999999-0.809000000000000050.58799999999999997-0.914000000000000030.40699999999999997-0.977999999999999980.20799999999999999-15.6700000000000004e-16-0.97799999999999998-0.20799999999999999-0.91400000000000003-0.40699999999999997-0.80900000000000005-0.58799999999999997-0.66900000000000004-0.74299999999999999-0.5-0.86599999999999999-0.309-0.95099999999999996-0.105-0.9950.105-0.9950.309-0.950999999999999960.5-0.865999999999999990.66900000000000004-0.742999999999999990.80900000000000005-0.587999999999999970.91400000000000003-0.406999999999999970.97799999999999998-0.207999999999999992.6291902682773483e-17-0.00065384615384583118-0.29530436727611131-0.407149903005388670.20468717553848803-0.459508829739425810.5000339709144086-0.0522824392313316210.20444165346679380.45895712720185405-0.373897837866523920.33573030516976354-0.56520453690073769-0.0591754798645009040.551409205813950140.3176130751959379-0.131386952017287120.62243041389833154-0.1444609125190702-0.683046048952053010.51896266786078527-0.46675660322909635-0.676779624251944490.22091656257348974-0.61723328415727241-0.355798314013887280.15396059620050356-0.72709003692960750.739171092208902420.077575788599018480.70634528325694623-0.230425116223339530.153903864563417180.726827008430388940.436562908690756420.60123684202225869-0.436389255152277280.60099824488402309-0.44644144213994863-0.6148097863548887-0.78986314125549695-0.0829374476327926770.48306021940714738-0.66512339089274719-0.717116495909639770.41374161364923850.71966081565749940.41521760024330068-0.172393981826836850.812553915742607380.83205471084583948-0.0874000174933099450.34031764454514540.76427847186026787-0.551646683990411390.34772453168234385-0.581198562733519110.5140347159468972-4.3977950621443274e-18-0.837382469604903370.72531999275454517-0.41850054828302335-0.34213203269662834-0.76835550876536474-0.62483523734501689-0.562383440182595780.564917319789501660.488146785922156210.389172534322699770.41806537561109980.237175459794329020.16912185462803336-0.025962281685754960.286169353657331540.39382883333149710.244348311961457930.54267064766876960.130311242657956930.380050968746640960.072009078853487340.24522905811277632-0.178196962665235210.22329825233727255-0.0049527105560499022-0.018125285606249857-0.28382039345137922-0.25412870038086932-0.1265506189486573-0.220012442533861740.124017969188261430.705945829365216990.24504049798687536-0.289239329419072540.68786790036032985-0.276731169405739050.49861546727758305-0.195804115204457540.33991543187901008-0.064806697817597680.458647757989873510.0398112094972438990.59660089503236968-0.0124525776233673070.75748185020432790.36004239628137363-0.54042597822871430.3198317338457855-0.703124164916731240.26175293701598690.318663722485640730.113649028007965720.34001459209079515-0.80612573707151369-0.26251531242500947-0.843587814987160180.088619865046719065-0.449043852146326190.458022307876681520.35347968445219013-0.0764834061563727660.43529169673866364-0.272014120989298390.33966312039894614-0.38995864108359496-0.441600509938721460.124737089769562650.22042251667383625-0.601090686607012480.01773990267719644-0.551951399569512670.577419672784286960.631721212662533450.836162963310296050.18384745602550973-0.421947736522365220.74460017792571787-0.15191200145175998-0.50793571011947636-0.28468886832575607-0.60129916940148143-0.15211412155416426-0.355825083579045190.0913285322740825940.853961344045157440.25268219760086502-0.82436443072803334-0.565221767943702820.65147861937265739-0.17960003622540516-0.845963804317787620.64262958611925292-0.57845575520448644-0.829776700511539620.27012165242582176-0.75613161550480157-0.436374694077371010.84682787874323973-0.275605448440036930.66460522715024917-0.0702268963894364250.271351871101858490.61028308052190072-0.5291687129171514-0.72847571124389077-0.46589500975993081-0.44182174455360729-0.44717423342713147-0.23275286852251509-0.91658877313709719-0.0963411202356542350.0712381705483739750.15038100559175938-0.075375181679099110.148115108318487010.14411045609387388-0.307201356571660030.046251548820153726-0.411673517544628780.080240926729569728-0.1479533766093776-0.090696855217268765-0.13701177299039563-0.14601986575596898-0.0018723983808613188-0.65587042533541751-0.19614678101972682-0.682228298443264380.0564770016520156070.5820334614227094-0.324312452746667420.56939947784766054-0.17974699975143729-0.330543338366073470.20546019299713006-0.542524154632108350.21357850028719924-0.397822414790862-0.038177513011169929-0.0015128273166917471-0.69402163173758880.12814058587287039-0.86580177941165415-0.30513569295646459-0.26238649633237193];functionT=get_triangles(dummy_arg)T=[518712335476246023206319211122036895810522233712311410818192011263449495251031313010929115881662501128815161415107531314271065247811119712171151613535911873119104981175931295284266911045275226349091291096128106275630124131603394831041112798078386869954013014556141075348657781028210296865859105569659664887511156548111101741324543769731171166738656810050621249312890707212691125296128829554782816757486564774978498978212262501001131718873793124835241528384585342535859596476543811981039818211056451095646747571115747128584236788949104595910413130449510560401256146106614110562221006240113635010863433854646454485866638646592338368667080397874121122856535687069110693476694572719034697072706611812073116726612072311011217410032113117311161144311375117799375127697638977645395577775510797836397778793612779678086806735111813980818180861028247398255924191839484441038484103268685356668858635818667116108871812387431245137107884249598958895991907170903471118929141125419283719291365812812993379433954484941193311811999954864964966419745397249779810987740959911910199409910032100991017312132991017102479102844601312610325981044913104126010523621054061106285210641881071553107426310819871084356109306110946126110345611046651113557111486211221631125011432132631134332114113129114748711517881155172116318611666791176711775122921183312011871101119733311995721207112112031120121731221213111712231122751291141234351123374288124511248861125414611012691126344612612579127751281273612412842128931271141293712974122406013044130602513124441311031011323274114132];