## Code for Figure 2.31

### main.m


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57% Legendre Galerkin solution of u''+u=-x,u(0)=u(1)=0
% MDG 1/4/13
N=10;
A=zeros(N,N);
u=zeros(N,1);
b=zeros(N,1);
for i=0:N-3
for j=0:N-1
A(i+1,j+1)=2./(2*j+1)*(i==j);
for k=0:j-1
if (mod(j-k,2)==0)
A(i+1,j+1)=A(i+1,j+1)+4*(j-k)*(j+k+1)*(i==k);
end
end
end
end
for j=0:N-1
A(N-1,j+1)=(-1)^j;
A(N,j+1)=1;
end
for i=0:N-3
b(i+1)=-(i==0)-1/3*(i==1);
end
b(N-1)=0;
b(N)=0;
c = A \ b;
%generate Legendre polynomials
nxs = 1000; % number of points
x = linspace(-1, 1, nxs); %points for plotting
soln=zeros(nxs,1);
dx=2/(nxs-1);
%generate basis functions
for m=0:N-1
P=legendre(m,x); % this generates associated legendres P_n^m from n=0 to m
phi(m+1,:)=P(1,:);                  % this picks out the m= 0 case
%(the legendre polynomial P_n)
end
for j=1:N
for i=1:nxs
soln(i)=soln(i)+c(j)*phi(j,i);
end
end
% plot((x+1)/2,soln) %plot in original coordinate system
% exact solution
Nex=100;
hex=1/(Nex+1);
uex=zeros(Nex+2,1);
xex=0:hex:1;
uex=-xex+sin(xex)*csc(1);
% hold on
% plot(xex,uex,'-')
% hold off

semilogy(abs(c),'o');

t = [[1:N]', abs(c)];
save LegendreGalerkin1.dat  t