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 above
white=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 triangulation
lw=1.4;
figure(1);clf;holdon;axisequal;axisoff;
forl=1:Nt
i=T(l,1);j=T(l,2);k=T(l,3);

% plot the three edges in each triangle
plot([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 issues
s=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 triangulation
F=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 Np
A=sparse(Np-Nb,Np-Nb);

% iterate through triangles 
forl=1:Nt
% fill in the matrix
i=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 matrix
ifi>Nb
A(i-Nb,i-Nb)=A(i-Nb,i-Nb)+gii+c*pii;
end
ifi>Nb&j>Nb
A(i-Nb,j-Nb)=A(i-Nb,j-Nb)+gij+c*pij;
A(j-Nb,i-Nb)=A(i-Nb,j-Nb);
end
ifi>Nb&k>Nb
A(i-Nb,k-Nb)=A(i-Nb,k-Nb)+gik+c*pik;
A(k-Nb,i-Nb)=A(i-Nb,k-Nb);
end
ifj>Nb
A(j-Nb,j-Nb)=A(j-Nb,j-Nb)+gjj+c*pjj;
end
ifj>Nb&k>Nb
A(j-Nb,k-Nb)=A(j-Nb,k-Nb)+gjk+c*pjk;
A(k-Nb,j-Nb)=A(j-Nb,k-Nb);
end
ifk>Nb
A(k-Nb,k-Nb)=A(k-Nb,k-Nb)+gkk+c*pkk;
end
% add the appropriate contributions to the right-hand terms
ifi>Nb
RHS(i-Nb)=RHS(i-Nb)+F(i)*pii+F(j)*pij+F(k)*pik;
end
ifj>Nb
RHS(j-Nb)=RHS(j-Nb)+F(i)*pij+F(j)*pjj+F(k)*pjk;
end
ifk>Nb
RHS(k-Nb)=RHS(k-Nb)+F(i)*pik+F(j)*pjk+F(k)*pkk;
end

end
% 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 boundary
U_calc=A\RHS;
U_calc=[0*(1:Nb)U_calc']';
% exact solution
u=inline('1-x.^2-y.^2','x','y');
U_exact=u(P(:,1),P(:,2));
% plot the computed solution
figure(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 points
x1=P1(1);y1=P1(2);
x2=P2(1);y2=P2(2);
x3=P3(1);y3=P3(2);
% triangle area
AT=abs((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1))/2;
% the integrals of lambda_i*lambda_j
pii=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_j 			 
gii=((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 solution
X=zeros(5*Nt,1);
Y=zeros(5*Nt,1);
Z=zeros(5*Nt,1);

forl=1:Nt
i=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 solution
plot3(X,Y,Z,'linewidth',lw,'color',color)
set(gca,'fontsize',15)

functionP=get_points(dummy_arg)
P=[10
0.977999999999999980.20799999999999999
0.914000000000000030.40699999999999997
0.809000000000000050.58799999999999997
0.669000000000000040.74299999999999999
0.50.86599999999999999
0.3090.95099999999999996
0.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.995
0.105-0.995
0.309-0.95099999999999996
0.5-0.86599999999999999
0.66900000000000004-0.74299999999999999
0.80900000000000005-0.58799999999999997
0.91400000000000003-0.40699999999999997
0.97799999999999998-0.20799999999999999
2.6291902682773483e-17-0.00065384615384583118
-0.29530436727611131-0.40714990300538867
0.20468717553848803-0.45950882973942581
0.5000339709144086-0.052282439231331621
0.20444165346679380.45895712720185405
-0.373897837866523920.33573030516976354
-0.56520453690073769-0.059175479864500904
0.551409205813950140.3176130751959379
-0.131386952017287120.62243041389833154
-0.1444609125190702-0.68304604895205301
0.51896266786078527-0.46675660322909635
-0.676779624251944490.22091656257348974
-0.61723328415727241-0.35579831401388728
0.15396059620050356-0.7270900369296075
0.739171092208902420.07757578859901848
0.70634528325694623-0.23042511622333953
0.153903864563417180.72682700843038894
0.436562908690756420.60123684202225869
-0.436389255152277280.60099824488402309
-0.44644144213994863-0.6148097863548887
-0.78986314125549695-0.082937447632792677
0.48306021940714738-0.66512339089274719
-0.717116495909639770.4137416136492385
0.71966081565749940.41521760024330068
-0.172393981826836850.81255391574260738
0.83205471084583948-0.087400017493309945
0.34031764454514540.76427847186026787
-0.551646683990411390.34772453168234385
-0.581198562733519110.5140347159468972
-4.3977950621443274e-18-0.83738246960490337
0.72531999275454517-0.41850054828302335
-0.34213203269662834-0.76835550876536474
-0.62483523734501689-0.56238344018259578
0.564917319789501660.48814678592215621
0.389172534322699770.4180653756110998
0.237175459794329020.16912185462803336
-0.025962281685754960.28616935365733154
0.39382883333149710.24434831196145793
0.54267064766876960.13031124265795693
0.380050968746640960.07200907885348734
0.24522905811277632-0.17819696266523521
0.22329825233727255-0.0049527105560499022
-0.018125285606249857-0.28382039345137922
-0.25412870038086932-0.1265506189486573
-0.220012442533861740.12401796918826143
0.705945829365216990.24504049798687536
-0.289239329419072540.68786790036032985
-0.276731169405739050.49861546727758305
-0.195804115204457540.33991543187901008
-0.064806697817597680.45864775798987351
0.0398112094972438990.59660089503236968
-0.0124525776233673070.7574818502043279
0.36004239628137363-0.5404259782287143
0.3198317338457855-0.70312416491673124
0.26175293701598690.31866372248564073
0.113649028007965720.34001459209079515
-0.80612573707151369-0.26251531242500947
-0.843587814987160180.088619865046719065
-0.449043852146326190.45802230787668152
0.35347968445219013-0.076483406156372766
0.43529169673866364-0.27201412098929839
0.33966312039894614-0.38995864108359496
-0.441600509938721460.12473708976956265
0.22042251667383625-0.60109068660701248
0.01773990267719644-0.55195139956951267
0.577419672784286960.63172121266253345
0.836162963310296050.18384745602550973
-0.421947736522365220.74460017792571787
-0.15191200145175998-0.50793571011947636
-0.28468886832575607-0.60129916940148143
-0.15211412155416426-0.35582508357904519
0.0913285322740825940.85396134404515744
0.25268219760086502-0.82436443072803334
-0.565221767943702820.65147861937265739
-0.17960003622540516-0.84596380431778762
0.64262958611925292-0.57845575520448644
-0.829776700511539620.27012165242582176
-0.75613161550480157-0.43637469407737101
0.84682787874323973-0.27560544844003693
0.66460522715024917-0.070226896389436425
0.271351871101858490.61028308052190072
-0.5291687129171514-0.72847571124389077
-0.46589500975993081-0.44182174455360729
-0.44717423342713147-0.23275286852251509
-0.91658877313709719-0.096341120235654235
0.0712381705483739750.15038100559175938
-0.075375181679099110.14811510831848701
0.14411045609387388-0.30720135657166003
0.046251548820153726-0.41167351754462878
0.080240926729569728-0.1479533766093776
-0.090696855217268765-0.13701177299039563
-0.14601986575596898-0.0018723983808613188
-0.65587042533541751-0.19614678101972682
-0.682228298443264380.056477001652015607
0.5820334614227094-0.32431245274666742
0.56939947784766054-0.17974699975143729
-0.330543338366073470.20546019299713006
-0.542524154632108350.21357850028719924
-0.397822414790862-0.038177513011169929
-0.0015128273166917471-0.6940216317375888
0.12814058587287039-0.86580177941165415
-0.30513569295646459-0.26238649633237193];
functionT=get_triangles (dummy_arg)
T=[5187123
35476
246023
206319
2111220
368958
1052223
37123114
1081819
2011263
449495
25103131
3010929
1158816
6250112
881516
1415107
531314
2710652
4781111
9712
1711516
135359
11873119
1049811
7593129
528426
6911045
275226
349091
2910961
2810627
56301
2413160
339483
1041112
798078
386869
9540130
14556
1410753
48657
78102
821029
686585
91055
6965
96648
8751115
6548111
10174132
4543
76973
11711667
386568
1005062
12493128
907072
12691125
296128
82955
478281
6757
486564
774978
498978
212262
50100113
171887
3793124
835241
528384
585342
535859
5964
765438
119810
398182
1105645
1095646
74757
1115747
1285842
367889
4910459
5910413
1304495
1056040
1256146
1066141
1056222
1006240
1136350
1086343
385464
64544
858666
386465
923383
686670
803978
74121122
856535
687069
1106934
766945
727190
346970
727066
11812073
1167266
1207231
10112174
10032113
11731116
11443113
7511779
9375127
697638
977645
395577
775510
797836
397778
7936127
796780
868067
3511181
398081
818086
1028247
398255
924191
839484
4410384
8410326
868535
666885
863581
8667116
1088718
1238743
1245137
1078842
495989
588959
919071
709034
7111892
9141125
419283
719291
3658128
1299337
943395
448494
11933118
1199995
486496
49664
19745
3972
497798
109877
409599
11910199
4099100
3210099
10173121
3299101
710247
91028
4460131
2610325
9810449
1310412
6010523
6210540
6110628
5210641
8810715
5310742
6310819
8710843
5610930
6110946
12611034
5611046
6511135
5711148
6211221
6311250
11432132
6311343
32114113
12911474
8711517
8811551
7211631
8611666
7911767
11775122
9211833
12011871
10111973
3311995
7212071
12112031
12012173
12212131
11712231
12275129
11412343
5112337
4288124
5112488
6112541
46110126
9112634
46126125
7912775
12812736
12412842
12893127
11412937
12974122
4060130
4413060
2513124
44131103
10113232
74114132];