Developed by Randy LeVeque for a course on Approximation Theory and Spectral Methods at the University of Washington.
See http://faculty.washington.edu/rjl/classes/am590a2013/codes.html for more IPython Notebook examples.
This Notebook presents some examples from Chapter 9 of "Spectral Methods in Matlab" by L. N. Trefethen. See http://people.maths.ox.ac.uk/trefethen/spectral.html for other m-files.
p22.m is also rewritten using Chebfun to avoid solving ill-conditioned Vandermonde matrices.
First we need to set things up to use Matlab and set the path to find chebfun.
You must have pymatbridge installed and chebfun available for this to work.
import pymatbridge
ip = get_ipython()
pymatbridge.load_ipython_extension(ip)
Chebfun is required -- to run this notebook you will need to download and add the correct directory to the path here, along with the directory of programs from "Spectral Methods in Matlab"
%%matlab
path(path,'/Users/rjl/chebfun/chebfun_v4.2.2889/chebfun')
path(path,'/Users/rjl/SMM/all_M_files')
%%matlab
% p21.m - eigenvalues of Mathieu operator -u_xx + 2qcos(2x)u
% (compare p8.m and p. 724 of Abramowitz & Stegun)
N = 42; h = 2*pi/N; x = h*(1:N);
D2 = toeplitz([-pi^2/(3*h^2)-1/6 ...
-.5*(-1).^(1:N-1)./sin(h*(1:N-1)/2).^2]);
qq = 0:.2:15; data = [];
for q = qq;
e = sort(eig(-D2 + 2*q*diag(cos(2*x))))';
data = [data; e(1:11)];
end
clf, subplot(1,2,1)
set(gca,'colororder',[0 0 1],'linestyleorder','-|--'), hold on
plot(qq,data), xlabel q, ylabel \lambda
axis([0 15 -24 32]), set(gca,'ytick',-24:4:32)
%%matlab
% p22.m - 5th eigenvector of Airy equation u_xx = lambda*x*u
for N = 12:12:48
[D,x] = cheb(N); D2 = D^2; D2 = D2(2:N,2:N);
[V,Lam] = eig(D2,diag(x(2:N))); % generalized ev problem
Lam = diag(Lam); ii = find(Lam>0);
V = V(:,ii); Lam = Lam(ii);
[foo,ii] = sort(Lam); ii = ii(5); lambda = Lam(ii);
v = [0;V(:,ii);0]; v = v/v(N/2+1)*airy(0);
xx = -1:.01:1; vv = polyval(polyfit(x,v,N),xx);
subplot(2,2,N/12), plot(xx,vv), grid on
title(sprintf('N = %d eig = %15.10f',N,lambda))
end
[Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.] [> In polyfit at 76 In web_eval at 44 In webserver at 241] [Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.] [> In polyfit at 76 In web_eval at 44 In webserver at 241]
Note the Matlab warnings from using polyfit (which uses the ill-conditioned Vandermonde matrix).
To avoid this, we rewrite Program 22 to use Chebfun in defining the Chebyshev points and D2 matrix, and to interpolate the resulting eigenfunction and create a smooth plot.
%%matlab
for N = 12:12:48
% Use Chebfun:
x = chebpts(N+1);
D2op = chebop(@(u) diff(u,2));
D2 = D2op(N+1);
D2 = D2(2:N,2:N);
[V,Lam] = eig(D2,diag(x(2:N))); % generalized ev problem
Lam = diag(Lam); ii = find(Lam>0);
V = V(:,ii); Lam = Lam(ii);
[foo,ii] = sort(Lam); ii = ii(5); lambda = Lam(ii);
v = [0;V(:,ii);0]; v = v/v(N/2+1)*airy(0);
% Use Chebfun to interpolate and plot:
d = domain(-1,1);
vpoly = interp1(x,v,d);
subplot(2,2,N/12), plot(vpoly), grid on
title(sprintf('N = %d eig = %15.10f',N,lambda))
end
%%matlab
% p23.m - eigenvalues of perturbed Laplacian on [-1,1]x[-1,1]
% (compare p16.m)
% Set up tensor product Laplacian and compute 4 eigenmodes:
N = 16; [D,x] = cheb(N); y = x;
[xx,yy] = meshgrid(x(2:N),y(2:N)); xx = xx(:); yy = yy(:);
D2 = D^2; D2 = D2(2:N,2:N); I = eye(N-1);
L = -kron(I,D2) - kron(D2,I); % Laplacian
L = L + diag(exp(20*(yy-xx-1))); % + perturbation
[V,D] = eig(L); D = diag(D);
[D,ii] = sort(D); ii = ii(1:4); V = V(:,ii);
% Reshape them to 2D grid, interpolate to finer grid, and plot:
[xx,yy] = meshgrid(x,y);
fine = -1:.02:1; [xxx,yyy] = meshgrid(fine,fine);
uu = zeros(N+1,N+1);
[ay,ax] = meshgrid([.56 .04],[.1 .5]); clf
for i = 1:4
uu(2:N,2:N) = reshape(V(:,i),N-1,N-1);
uu = uu/norm(uu(:),inf);
uuu = interp2(xx,yy,uu,xxx,yyy,'spline'); %%% changed from 'cubic' to 'spline' to avoid matlab error
subplot('position',[ax(i) ay(i) .38 .38])
contour(fine,fine,uuu,-.9:.2:.9)
colormap(1e-6*[1 1 1]); axis square
title(['eig = ' num2str(D(i)/(pi^2/4),'%18.12f') '\pi^2/4'])
end
%%matlab
% p24.m - pseudospectra of Davies's complex harmonic oscillator
% (For finer, slower plot, change 0:2 to 0:.5.)
% Eigenvalues:
N = 70; [D,x] = cheb(N); x = x(2:N);
L = 6; x = L*x; D = D/L; % rescale to [-L,L]
A = -D^2; A = A(2:N,2:N) + (1+3i)*diag(x.^2);
clf, plot(eig(A),'.','markersize',14)
axis([0 50 0 40]), drawnow, hold on
% Pseudospectra:
x = 0:2:50; y = 0:2:40; [xx,yy] = meshgrid(x,y); zz = xx+1i*yy;
I = eye(N-1); sigmin = zeros(length(y),length(x));
%h = waitbar(0,'please wait...'); %%% removed waitbar commands
for j = 1:length(x),
for i = 1:length(y), sigmin(i,j) = min(svd(zz(i,j)*I-A)); end
end
contour(x,y,sigmin,10.^(-4:.5:-.5)), colormap(1e-6*[1 1 1]);