Assignment 6 Scripts
P1
%P1
% Explore accuracy of AB2 method with fixed step size on an
% orbit problem. The orbiting body has period 2pi. So let's
% explore whether or not we get back to the original position
% at T = 2pi.
close all
abserr = .000001;
for e = [.1 .5 .9]
% e = eccentricity. 0<= e < 1. e=0 corresponds to a circle. e = .9 is "cigarlike".
x0 = 1-e; xp0 = 0; y0 = 0; yp0 = sqrt((1+e)/(1-e));
Tfinal = 2*pi;
% Find out how many f-evals are required by ODE23 to get back to
% original position with abs error <= .000001
[xfinal,yfinal,FunEvals] = ODE23Orbit(Tfinal,x0,xp0,y0,yp0,abserr);
figure
for N = [50 100 200 500 1000 2000 4000 8000]
[x,y] = AB2Orbit(Tfinal,N,x0,xp0,y0,yp0);
plot(x,y,':')
axis([-1.2 1.2 -1.2 1.2])
axis equal
title(sprintf('Eccentricity = %5.2f N = %1d ODE23 F-evals = %1d',e,N,FunEvals))
xlabel(sprintf('|x(0) - x(2pi)| = %10.6f |y(0) - y(2pi)| = %10.6f ',abs(x(1)-x(N+1)),abs(y(1)-y(N+1))))
pause
end
end
P2
% P2
% Poisson on an L-shaped region
N = 5;
u = LShape(N,@P2F,@P2G);
clc
k = 0;
for i=1:3*N-1
if i<=2*N
jmax = N-1;
else
jmax = 2*N-1;
end
disp(sprintf('%7.2f ',u(k+1:k+jmax)))
k = k+jmax;
end
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
function z = P2F(x,y)
z = exp(-(x^2+y^2);
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
function z = P2G(x,y)
z = (1+cos(2*x))*sin(x^2+1)*sin(x*y)^2
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%P2Hints
clc
format short
%%%%%%%%%%%%%%%%%%
% The Schur Decomposition of the 1-2-1 matrices.
% The starting point is the matrix E_m, which is m-by-m,
% zeros everywhere except ones along the sub and superdiagonal.
m = 6;
E_m = diag(ones(m-1,1),-1)+diag(ones(m-1,1),1)
% Let Q_m be the m-by-m matrix with (i,j) entry = sin(i*j*pi/(m+1)) / sqrt(2/(m+1))
% This matrix is orthogonal an diagonalizes E_m.Indeed,
% Q_m' * E_m * Q_m = diag(g1,...,gm) gj = 2*cos(j*pi/(m+1))
% and so since S_m = 2eye(m,m) - E_m
% Q_m' * S_m * Q_m = diag(d1,...,dm) dj = 2(1 - cos(j*pi/(m+1))
S_m = 2*eye(m,m) - E_m
Q_m = sqrt(2/(m+1)) * sin( (1:m)'*((1:m)*(pi/(m+1))) )
D_m = Q_m' * S_m * Q_m
disp(' ')
disp('The matrix S_m is positive definite.')
disp(' ')
disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>')
disp(' ')
disp('Fast and slow Poisson equation solving.')
m = 3; n = 4;
disp(' ')
disp('Solve L_m,n u = b, L_m,n = kron(eye(m,m),S_n) + kron(S_m,eye(n,n))')
b = randn(m*n,1);
uF = FastPS(b,m,n);
uS = SlowPS(b,m,n);
format long
Sols = [uF uS]
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Display the matrix of coeffients A associated with the solving of Poisson's equation
% on an L-shaped region
N = 5;
A = MakeL_matrix(N);
% There are P unknowns (= number of interior meshpoints) where
P = 4*N^2 - 5*N + 1;
% There are M unknowns above the line segment from (0,1) to (1,1).
% There are M unknowns below the line segment from (0,1) to (1,1).
M = (N-1)*(2*N-1);
% Display A ...
fill([.5 P+.5 P+.5 .5 .5],[.5 .5 P+.5 P+.5 .5],'y')
axis ij equal off
hold on
for i=1:P
for j=1:P
if A(i,j)
plot(i,j,'.')
end
end
end
% Partition A into 3-by-3. The (1,1) and (3,3) blocks are M-by-M and are
% discrete Laplacian matrices.
plot([.5 P+.5],[M+.5 M+.5],'r')
plot([.5 P+.5],[M+N-.5 M+N-.5],'r')
plot([M+.5 M+.5],[.5 P+.5],'r')
plot([M+N-.5 M+N-.5],[.5 P+.5],'r')
text(-5,M/2,'M'); text(-5,M+N/2,'N-1'); text(-5,M+N+M/2,'M')
title(sprintf('The L-system matrix: N = %1d, M = %1d, P = %1d',N,M,P))
hold off
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
function A = MakeL_Matrix(N)
% Sets up the matrix of coefficients when Poisson's equation is solved
% on the L-shaped region defined by (0,0), (2,0), (2,1), (1,1), (1,3), (0,3)
% There are P unkowns (= number of interior meshpoints) where
P = 4*N^2 - 5*N + 1;
% Let S be the line segment from (0,1) to (1,1)
% There are M meshpoints above S, N-1 on S, and M meshpoints below S where
M = (N-1)*(2*N-1);
nT = N-1; % # meshpoints per horizontal line in the top rectangle.
nB = 2*N-1; % # meshpoints per horizontal line in the bottom rectangle.
% Set up A
A = 4*eye(P,P);
k=0;
% Top Rectangle
for i = 1:2*N
for j=1:N-1
k = k+1;
% Finish setting up the k-th row
if j>1
A(k,k-1) = -1; % There is a west neighbor
end
if j<N-1
A(k,k+1) = -1; % There is an east neighbor
end
if i>1
A(k-nT,k) = -1; % There is a north neighbor
end
A(k+nT,k) = -1; % There always is a south neighbor.
end
end
% Set up rows associated with the meshpoints on S
for j=1:2*N-1
k = k+1;
if j>1
A(k,k-1) = -1; % There is a west neighbor
end
if j<2*N-1
A(k,k+1) = -1; % There is an east neighbor
end
if j<=N-1
A(k-nT,k) = -1; % There is a north neighbor
end
A(k+nB,k) = -1; % There is a south neighbor
end
% The bottom rectangle
for i=2:N-1
for j=1:2*N-1
k = k+1;
if j>1
A(k,k-1) = -1; % There is a west neighbor
end
if j<2*N-1
A(k,k+1) = -1; % There is an east neighbor
end
A(k-nB,k) = -1; % There is a north neighbor
if i<N-1
A(k+nB,k) = -1; % There is a south neighbor
end
end
end
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
function u = FastPS(b,m,n)
% Solves L_mn u = b
% The eigenvectors and values of S_m
Q_m = sqrt(2/(m+1)) * sin( (1:m)'*((1:m)*(pi/(m+1))) );
D_m = 2*(1- cos((1:m)'*pi/(m+1)));
% The eigenvectors and values of S_n
Q_n = sqrt(2/(n+1)) * sin( (1:n)'*((1:n)*(pi/(n+1))) );
D_n = 2*(1 - cos((1:n)'*pi/(n+1)));
% L_mn u = b transforms to Dy = c where D = kron(eye(m,m),diag(D_n)) + kron(diag(D_m),eye(n,n))
% y = Q'*u, and c = Q'*b where Q = kron(Q_m,Q_n).
c = reshape(Q_n'*reshape(b,n,m)*Q_m ,m*n,1); % (n-by-n) times (n-by-m) times (m-by-m)
D = kron(ones(m,1),D_n) + kron(D_m,ones(n,1));
y = c./D;
u = reshape( Q_n*reshape(y,n,m)*Q_m' , m*n,1 ); % (n-by-n) times (n-by-m) times (m-by-m)
% PS Kronecker products times vectors can be done fast as above.
% Total flops = 4(nm^2 + mn^2) flops
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
function u = SlowPS(b,m,n)
% Solves L_mn u = b
S_m = diag(-ones(m-1,1),-1) + 2*eye(m,m) + diag(-ones(m-1,1),1);
S_n = diag(-ones(n-1,1),-1) + 2*eye(n,n) + diag(-ones(n-1,1),1);
A = kron(eye(m,m),S_n) + kron(S_m,eye(n,n));
u = A\b;
% O( (mn)^3 ) flops