ASSIGNMENT 2 SOLUTIONS

ASSIGNMENT 2 SOLUTIONS
Nick Trefethen, 4 November 2014
1. CONDITION NUMBERS
We set up anonymous functions to construct the matrix A as indicated:
clear all, close all
D = @(n) sparse(toeplitz([2 -1 zeros(1,n-3)]));
I = @(n) speye(n-1);
A = @(n) kron(I(n),kron(I(n),D(n))) + kron(I(n),kron(D(n),I(n))) + ...
kron(D(n),kron(I(n),I(n)));
Here is a table of condition numbers:
nn = 3:10;
for n = nn
c(n-2) = cond(full(A(n)));
fprintf(’n = %2d (n-1)^3 = %3d
end
n
n
n
n
n
n
n
n
= 3
= 4
= 5
= 6
= 7
= 8
= 9
= 10
(n-1)^3
(n-1)^3
(n-1)^3
(n-1)^3
(n-1)^3
(n-1)^3
(n-1)^3
(n-1)^3
=
=
=
=
=
=
=
=
8
27
64
125
216
343
512
729
cond(A(n))
cond(A(n))
cond(A(n))
cond(A(n))
cond(A(n))
cond(A(n))
cond(A(n))
cond(A(n))
cond(A(n)) = %6.3f\n’,n,(n-1)^3,c(n-2))
=
=
=
=
=
=
=
=
3.000
5.828
9.472
13.928
19.196
25.274
32.163
39.863
The results for n = 5 and n = 10 suggest that the condition number is approximately 0.4n2 , and this is confirmed by a plot. (The actual proportionality
constant as n → ∞ is known to be 4/π 2 .)
LW = ’linewidth’; MS = ’markersize’; FS = ’fontsize’;
hold off, loglog(nn,c,’.-’,LW,1.2,MS,16), grid on
title(’Condition numbers’,FS,16)
xlabel(’n’,FS,16)
hold on, loglog(nn,.4*nn.^2,’--r’,LW,2)
text(5.1,23,’0.4 n^2’,’color’,’r’,FS,16)
axis([3 10 2 50])
set(gca,’xtick’,3:10,’ytick’,[5 10:10:50])
1
Condition numbers
50
40
30
0.4 n2
20
10
5
3
4
5
6
7
8
9
10
n
2. SOLUTION BY DIRECT METHODS
To my eye the time seems to grow in proportion to n4 . I do not know if this
is the ”right” answer, for I do not know the details of the algorithms from
UMFPACK used by Matlab’s backslash command.
nn = [5:40]; tt = [];
for n = nn
N = (n-1)^3;
b = ones(N,1);
AA = A(n);
tic, AA\b; tt = [tt toc];
end
hold off, loglog(nn,tt,’.-’,LW,1.2,MS,16), grid on
title(’Times for direct solution’,FS,16)
xlabel(’n’,FS,16)
hold on, loglog(nn,5e-7*nn.^4,’--r’,LW,2)
text(11,0.2,’5e-7 n^4’,’color’,’r’,FS,16)
axis([5 50 1e-4 30]), set(gca,’xtick’,[5 10:10:50])
Times for direct solution
0
10
4
5e−7 n
−2
10
−4
10
5
10
20
30
40
50
n
The dimension of this matrix is n3 , so the time would be proportional to n9 ,
which is hugely worse, if backslash did not take advantage of sparsity.
2
3. SOLUTION BY CG
Our CG experiment is, please note, CG without a preconditioner, which is far
from the best iterative method available. Nevertheless it is faster than backslash.
Again the time seems to scale in proportion to n4 , but the constant is better.
Note that the dashed red line is at the same height as before.
hold on, nn = 5:50; tt = [];
for n = nn
N = (n-1)^3;
b = ones(N,1);
AA = A(n);
tic, pcg(AA,b,1e-10,500); tt = [tt toc];
end
hold off, loglog(nn,tt,’.-’,LW,1.2,MS,16), grid on
title(’Times for CG solution’,FS,16)
xlabel(’n’,FS,16)
hold on, loglog(nn,5e-7*nn.^4,’--r’,LW,2)
text(21,.003,’5e-7 n^4’,’color’,’r’,FS,16)
axis([5 50 1e-4 30]), set(gca,’xtick’,[5 10:10:50])
pcg
pcg
pcg
pcg
pcg
pcg
...
pcg
pcg
pcg
pcg
converged
converged
converged
converged
converged
converged
at
at
at
at
at
at
iteration
iteration
iteration
iteration
iteration
iteration
4 to a solution with relative residual 4.2e-16.
7 to a solution with relative residual 9.1e-16.
10 to a solution with relative residual 1.5e-15.
16 to a solution with relative residual 1.8e-15.
20 to a solution with relative residual 2.6e-15.
24 to a solution with relative residual 8.4e-11.
converged
converged
converged
converged
at
at
at
at
iteration
iteration
iteration
iteration
131
134
136
140
to
to
to
to
a
a
a
a
solution
solution
solution
solution
with
with
with
with
relative
relative
relative
relative
residual
residual
residual
residual
Times for CG solution
0
10
−2
10
5e−7 n4
−4
10
5
10
20
30
40
50
n
We can explain the n4 timing as follows. From problem 1, we know that the
condition number is proportional to n2 . Thus CG theory tells us that the
3
7.5e-11.
8.3e-11.
9.5e-11.
7.6e-11.
number of steps for convergence will scale at worst linearly with n, and this
prediction is confirmed by the numbers printed by pcg during execution. Each
step, however, requires work proportional to n3 , since the biggest part of the
work is a matrix-vector multiplication and the matrix has dimension O(n3 ) with
a constant number of nonzeros on each row.
4. A MATRIX OF DIMENSION 40000
There are many ways to solve this. A very simple one is to use preconditioned
CG and take the diagonal matrix as the preconditioner. Convergence takes
place in just a few iterations.
N = 40000;
p = primes(10+round(1.5*N*log(N)));
% more than enough!
d = sparse(p(1:N));
A = diag(d);
I = 2.^(0:floor(log2(N)));
O = ones(N,2*length(I));
A = spdiags(O,[I -I],A);
b = (1:N)’; tol = 1e-14; maxit=100;
tic, x = pcg(A,b,tol,maxit,diag(d)); toc
x10000 = x(N/2)
pcg converged at iteration 12 to a solution with relative residual 7.5e-15.
Elapsed time is 0.163716 seconds.
x10000 =
0.088980964585111
4

Download Report