Figure 2.31:

Dependence of \left | c(j) \right | on j for the Legendre-Galerkin approximation of \eqref {eq:MWRexample} with n=10.

Code for Figure 2.31

Text of the GNU GPL.

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