Figure 2.6:

Function f(x)=\exp \big (-8x^{2}\big ) and truncated Legendre-Fourier series approximations with n=2,5,10.

Code for Figure 2.6

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
% Legendre series of a Gaussian function
% (after JBR solution)
%MDG 12/3/12

nxs = 1000; % number of points
x = linspace(-1, 1, nxs); %points for plotting
dx=2/(nxs-1);
nterms = [2,5,10]; % highest polynomial power for sums
u = zeros(length(nterms), nxs);
ntermsmax=nterms(length(nterms));
phi=zeros(ntermsmax,nxs);
%generate basis functions
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
% generate function at points x
flist=exp(-8*x.^2);
% plot(x,flist)
for i = 1:length(nterms)
  for k = 0:nterms(i);
      % compute coefficients by simple rectangle rule sum
      % (quick and dirty, but inaccurate for large k, because of the
      % oscillations of the basis functions)

        ck = dot(flist(:),phi(k+1,:))*dx;
        u(i,:) = u(i,:)+ ck*phi(k+1,:);
  end
end
plot(x,flist,'-',x,u(1,:),'--',x,u(2,:),':',x,u(3,:),'-.');
legend('exp(-8x^2)','n=2','n=5','n=10');

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