Building Big Matrices: sparse matrix

Let's consider the matrices associated to a Poisson problem:

in a square domain

Using finite differences you obtain a linear system of equations with a matrix A of the form:

Build by hand the blockMatrix

One possibility is to build the matrix step by step:

c=2;

a=-1;

b=-1;

N=3;

Amd=2*c*eye(N); %main diagonal block matrix

Amd(2,1)=a; Amd(1,2)=a;

Amd(3,2)=a; Amd(2,3)=a;

Aupp = b*eye(N);

Alow = b*eye(N);

Azer= zeros(N);

A=[Amd, Aupp, Azer;

Alow, Amd, Aupp;

Azer, Alow,Amd];

spy(A)

Better using idiag function (not in Matlab distribution)

Amd=zeros(N);

Amd(idiag(size(Amd),0))=4;

Amd(idiag(size(Amd),-1))=a;

Amd(idiag(size(Amd),1))=a;

Aupp=zeros(N);

Aupp(idiag(size(Aupp),0))=b;

Alow=zeros(N);

Alow(idiag(size(Alow),0))=b;

Azer= zeros(N);

A=[Amd, Aupp, Azer;

Alow, Amd, Aupp;

Azer, Alow,Amd];

spy(A) %notice that they are sparse matrices

Sparse Matrices

A=zeros(10);

A(1,2)=3;

A(2,3)=4;

A(3,1)=1;

A(6,9)=1;

A(7,1)=2;

A(9,8)=2;

As=sparse(A);

spy(As);

Automatic blockMatrix using blktridiag function (not in Matlab distribution)

a=-1;

b=-1;

c=4;

N=2^3;

% define the matrix blocks

Amd = blktridiag(c,a,a,N);

%Alow = b*eye(N);

Alow =blktridiag(b,0,0,N);

%Aupp = b*eye(N);

Aupp =blktridiag(b,0,0,N);

A = blktridiag(Amd,Alow,Aupp,N);

spy(A)

Poisson's equation: Finitte Differences solution

Let's the Poisson's problem:

in a square domain

with homogeneous Dirichlet boundary conditions:

and

The analytical solution is

In this case:

Linear system

with

, where

is the domain grid size.

In this exemple the number of subdivisions is N+1 (here 8), so the total grid, including the boundaries, has (N+2)x(N+2) points (here 9x9). The number of unknowns per row is N (here 7), which is the size of each bock matrix and finally, the global matrix A has size

(here 49x49).

% define the parameters 
a = -1;
b = -1;
c = 4;
N = 7; %dimension of the block matrices (correspond to the number of unknowns in each row)
%
% define the matrix blocks
%
Amd = blktridiag(c,a,a,N);
Alow = blktridiag(b,0,0,N);
Aupp = blktridiag(b,0,0,N);
A = blktridiag(Amd,Alow,Aupp,N);
% Solve the system
f = @(x,y) -2*pi^2*(cos(2*pi*x).*sin(pi*y).^2+cos(2*pi*y).*sin(pi*x).^2);
u = @(x,y) (sin(pi*x).^2).*sin(pi*y).^2;
h = 1/(N+1); % corresponds to N+1 subdivisions 
x = 0:h:1;
y = x;
[X,Y] = meshgrid(x,y);
rhs = f(X(2:end-1,2:end-1),Y(2:end-1,2:end-1));
solExact = u(X(2:end-1,2:end-1),Y(2:end-1,2:end-1));
solExact = solExact(:);
bb = h^2*rhs(:); %independent vector
xini = zeros(size(A,1),1);
sol = conjgrad(A,bb,xini,1.e-14);
norm(bb-A*sol)
ans = 
   2.4023e-15
norm(sol-solExact,inf)
ans = 
   5.3029e-02

Comparing efficiency with the solution using full matrices

a = -1;
b = -1;
c = 4;
N = 128;  %bigger dimensions
%
% define the matrix blocks
%
tic
Amd = blktridiag(c,a,a,N);
Alow = blktridiag(b,0,0,N);
Aupp = blktridiag(b,0,0,N);
A = blktridiag(Amd,Alow,Aupp,N);
h = 1/(N+1); % corresponds to N+1 subdivisions 
x = 0:h:1;
y = x;
[X,Y] = meshgrid(x,y);
% substitute all grid points but the boundary ones
rhs = f(X(2:end-1,2:end-1),Y(2:end-1,2:end-1));
bb = h^2*rhs(:); %independent vector
tic
solSparse = A\bb; %solution using sparse matrices
toc
Elapsed time is 0.020963 seconds.
solExact = u(X(2:end-1,2:end-1),Y(2:end-1,2:end-1));
solExact = solExact(:);
norm(solSparse-solExact,inf)
ans = 
   1.9766e-04
Cf = full(A);
tic
solFull = Cf\bb; %solution using full matrices
toc
Elapsed time is 8.254292 seconds.