Figure 2.7:

Function f(x)=H(x) and truncated Legendre-Fourier series approximations with n=10,50,100.

Code for Figure 2.7

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
% Legendre series of the step function
% (after JBR solution)
%MDG 11/25/12
nxs = 301; % number of points
x = linspace(-1, 1, nxs); %points for plotting
nterms = [10, 50, 100]; % highest polynomial power for sums
u = zeros(length(nterms), nxs);
ntermsmax=nterms(length(nterms));
phi=zeros(ntermsmax,nxs);
for m=0:ntermsmax
    P=legendre(m,x); % this generates associated legendres P_n^m from n=0 to m
    phi(m+1,:)=P(1,:)*sqrt((2*m+1)/2); % this picks out the m= 0 case
                                       %(the legendre polynomial P_n)
                                       % and also applies the normalization
end
% n=1;
% plot(x,phi(n+1,:)) % plot nth basis fn (recalling that indices start at
% % 0
flist=(sign(x)+1)/2; % step function
for i = 1:length(nterms)
  k=0;  %k=0 is the only even value of k that contributes
  un0 = sqrt(pi)/(2*gamma(1-k/2)*gamma((3+k)/2))*sqrt((2*k+1)/2);
  u(i,:) = u(i,:)+ un0*phi(k+1,:);
  for k = 1:2:nterms(i); % sum over odd values of k (odd Legendre polys)
      %integrals from Mathematica -- gamma Functions!
        un0 = sqrt(pi)/(2*gamma(1-k/2)*gamma((3+k)/2))*sqrt((2*k+1)/2);
        u(i,:) = u(i,:)+ un0*phi(k+1,:);
  end
end
plot(x,flist,'-',x,u(1,:),'--',x,u(2,:),':',x,u(3,:),'-.');
legend('H(x)','n=10','n=50','n=100');

%save results for plotting with gnuplot
table = [x' flist' u'];
save LegendreStep.dat table