Iterative Methods

Functions

Here, we recall the Laplace finite difference operators from Chapter 4.

DXmat[Nx_, Ny_, dx_, dy_, left_, right_] := Module[{g,retmat,bside, index, indey,i,j,k},
g[n_,m_] := (m - 1)*Nx + n;

retmat = Table[0.0,{i,1,Nx Ny},{j,1,Nx Ny}];
bside = Table[0.0,{i,1,Nx Ny}];

For[indey = 1, indey ≤ Ny, indey = indey+1,
j = g[1,indey];
k = g[2,indey];

retmat[[j,j ]] += -2.0/dx^2;
retmat[[j,k ]] += 1.0/dx^2;
bside[[ j ]] += -left[indey dy]/dx^2;

For[index = 2, index ≤ Nx-1, index = index+1,
i = g[index-1,indey];
j = g[index,indey];
k = g[index + 1,indey];
retmat[[j,i]] += 1.0/dx^2;
retmat[[j,j]] += -2.0/dx^2;
retmat[[j,k]] += 1.0/dx^2;
];
i = g[Nx-1,indey];
j = g[Nx,indey];
retmat[[j,i ]] += 1.0/dx^2;
retmat[[ j,j]] += -2.0/dx^2;
bside[[ j]] += -right[indey dy]/dx^2;
];
Return[{retmat, bside}];
]

DYmat[Nx_, Ny_, dx_, dy_, top_, bottom_] := Module[{g,retmat,bside, index, indey,i,j,k},
g[n_,m_] := (m - 1)*Nx + n;

retmat = Table[0.0,{i,1,Nx Ny},{j,1,Nx Ny}];
bside = Table[0.0,{i,1,Nx Ny}];

For[index = 1, index ≤ Nx, index = index+1,
j = g[index,1];
k = g[index, 2];

retmat[[j,j ]] += -2.0/dy^2;
retmat[[j,k ]] += 1.0/dy^2;
bside[[ j ]] += -bottom[index dx]/dy^2;

For[indey = 2, indey ≤ Ny-1, indey = indey+1,
i = g[index,indey-1];
j = g[index,indey];
k = g[index ,indey+1];
retmat[[j,i]] += 1.0/dy^2;
retmat[[j,j]] += -2.0/dy^2;
retmat[[j,k]] += 1.0/dy^2;
];
i = g[index, Ny-1];
j = g[index, Ny];
retmat[[j,i ]] += 1.0/dy^2;
retmat[[ j,j]] += -2.0/dy^2;
bside[[ j]] += -top[indey dy]/dx^2;
];
Return[{retmat, bside}];
]

MakeSource[ Nx_, Ny_, dx_, dy_, source_] := Module[{svec, g,j, index, indey},
g[n_,m_] := (m - 1)*Nx + n;

svec = Table[0.0,{i,1,Nx Ny}];
For[index = 1, index ≤ Nx, index=index+1,
For[indey = 1, indey ≤ Ny, indey = indey + 1,
j = g[index,indey];
svec[[j ]] = source[ index dx, indey dy];
];
];
Return[svec];
]

PlotPoisson[invec_, Nx_,dx_, Ny_, dy_] := Module[{g, Datatable,retplot},
g[n_,m_] := (m - 1)*Nx + n;
Datatable= Table[{ii dx, jj dy, invec[[ g[ii,jj] ]]},{ii,1, Nx}, {jj,1,Ny}];
Datatable = Partition[Flatten[Datatable],3];
 retplot = ListContourPlot[Datatable,Contours-> 20, ColorFunction->"GrayTones"];
Return[retplot];
]

Jacobi Method

The Jacobi method proceeds by the iteration in (11.6), where P is the diagonal of the matrix A and Q is “everything else”, Q = P - A.  The function below takes “Nsteps” steps, in order to keep track of the convergence, but we could also exit when a certain residual requirement was met (like Max[Abs[res]] < epsilon for given epsilon, for example)

JacobiInverse[inA_, inb_,Nsteps_] := Module[{retx, P,Pinv,Q,S,c,index},
P = DiagonalMatrix[Table[inA[[jj,jj]],{jj,1,Length[inA]}]];
Pinv = DiagonalMatrix[Table[1.0/inA[[jj,jj]],{jj,1,Length[inA]}]];
Q = P - inA;
S = Pinv.Q;
c = Pinv.inb;
retx = inb;
For[index = 1, index ≤ Nsteps, index=index+1,
retx = S.retx + c;
];
Return[retx];
]

The Jacobi method works well when the diagonal entries are large - below, we’ll make a random matrix with larger entries along the diagonal:

SeedRandom[11111];
Ns = 10;
A = N[Table[RandomReal[{-1/2,1/2}],{i,1,Ns},{j,1,Ns}]];
A = A + DiagonalMatrix[Table[RandomReal[{-10,10}],{j,1,Ns}]];
b = Table[1.,{i,1,Ns}];

MatrixPlot[A]

The residual, after 10 iterations:

tstx = JacobiInverse[A,b,10];
Print[A.tstx - b];

If we run out further, to 1000 iterations, then we can make the residual small.

tstx = JacobiInverse[A,b,1000];
Print[A.tstx - b];

Not every matrix will lead to a converging sequence for the Jacobi method, one has to check that the S matrix in the method has maximum eigenvalue less than 1 (in absolute value).

Poisson Problem - Figure 11.1

Here, we set up the same parallel plate capacitor with grounded boundaries as in Chapter 4.  This time, we’ll use the JacobiInverse function to perform the matrix inversion (instead of the built-in LinearSolve or Inverse functions).

RHS[x_,y_] := Module[{retval, rho},
rho = 1;
retval = 0;
If[y == Round[Ny/4] dy && (x ≤(Nx -  Round[Nx/5]+1) dx && x ≥ (Round[Nx/5] dx)) , 
retval = -rho;
];
If[y == (Ny - Round[Ny/4]+1) dy && (x ≤(Nx -  Round[Nx/5]+1) dx && x ≥ (Round[Nx/5] dx)) , 
retval = rho;
];
Return[-retval];
]

top[x_] := 0

bottom[x_] := 0

left[y_] := 0

right[y_] := 0

Make a modest grid,

Nx =40;
Ny = 40;
dx = 1/(Nx+ 1);
dy = 1/(Ny+1);

Construct the matrix, boundary, and source vectors:

{DX, bx} = DXmat[Nx, Ny, dx, dy, left, right] ;
{DY, by} = DYmat[Nx, Ny, dx, dy, top, bottom];
bs = MakeSource[ Nx, Ny, dx, dy,RHS] ;
rhsvec = bx + by + bs;
D2mat = DX + DY;

Here’s the previous solution, using LinearSolve

fxy = LinearSolve[D2mat, rhsvec];

And next, a table of Jacobi solutions at different levels of iteration.

Jxytabs = Table[JacobiInverse[D2mat, rhsvec,jj],{jj,200,1000,200}];

Gtab = Table[PlotPoisson[Jxytabs[[jj]],Nx,dx,Ny,dy],{jj,1,Length[Jxytabs]}];

Ga =PlotPoisson[fxy,Nx,dx,Ny,dy]

Grout = GraphicsGrid[{{Gtab[[1]],Gtab[[2]],Gtab[[3]]},{Gtab[[4]],Gtab[[5]],Ga}}]

Successive Over Relaxation

Below, we implement the matrix decomposition and iteration associated with (11.18) - this has a parameter w in it, and that determines whether we are 
iterating via under-relaxation (w < 1), Gauss-Seidel (w = 1) or over-relaxation (w > 1).  Here, we make the matrices as usual, and then convert them to sparse matrices.  Then the lower triangular inverse required by the method is performed using the built-in function LinearSolve (we could just do the back-substitution as in Chapter 9) with sparse methods, that speeds up the iteration without changing its basic form.

SORInverse[inA_, inb_, w_, Nsteps_] := Module[{matD, matL, matU,LHSmat,RHSmat,RHSvec,retx,index},
matD = DiagonalMatrix[Table[inA[[jj,jj]],{jj,1,Length[inA]}]];
matL =- Table[If[j >  k,inA[[j,k]],0],{j,1,Length[inA]},{k,1,Length[inA]}];
matU =  -Table[If[j < k,inA[[j,k]],0],{j,1,Length[inA]},{k,1,Length[inA]}];
LHSmat = matD - w matL;
RHSmat = ( (1-w) matD + w matU);
RHSvec = w inb;
LHSmat = SparseArray[LHSmat];
RHSmat = SparseArray[RHSmat];
retx = inb;
For[index = 1, index ≤ Nsteps, index=index+1,
retx = LinearSolve[LHSmat,RHSmat.retx + RHSvec];
];
Return[retx];
]

Here’s our test matrix from above

SeedRandom[11111];
Ns = 10;
A = N[Table[RandomReal[{-1/2,1/2}],{i,1,Ns},{j,1,Ns}]];
A = A + DiagonalMatrix[Table[RandomReal[{-10,10}],{j,1,Ns}]];
b = Table[1.,{i,1,Ns}];

The residual, using the inverse approximant, at ten iterations with w = 1.5

tstx = SORInverse[A,b,1.5,10];
Print[A.tstx - b];

If we make w = 1.1, the residual goes down, even at ten iterations:

tstx = SORInverse[A,b,1.1,10];
Print[A.tstx - b];

and if we run 20 steps at this w = 1.1, the residual is acceptable

tstx = SORInverse[A,b,1.1,20];
Print[A.tstx - b];

Laplace Problem - Figure 11.2

One of the examples from Section 11.2.2 - we have an empty region of charge, with prescribed boundary conditions.

RHS[x_,y_] := 0

top[x_] := Sin[3 Pi x]

bottom[x_] := Sin[Pi x]

left[y_] := Sin[2 Pi y]

right[y_] := Sin[5 Pi y]

Set up the grid,

Nx= 40;
Ny= 40;
dx = 1/(Nx + 1);
dy = 1/(Ny+1);

and construct the finite difference matrix, and source/boundary vectors:

{DX, bx} = DXmat[Nx, Ny, dx, dy, left, right] ;
{DY, by} = DYmat[Nx, Ny, dx, dy, top, bottom];
bs = MakeSource[ Nx, Ny, dx, dy,RHS] ;
rhsvec = bx + by + bs;
D2mat = DX + DY;

here’s the solution using the built-in LinearSolve function:

fxy = LinearSolve[D2mat, rhsvec];

and calculated using the SORinverse, at various levels of iteration.  Here, we are employing Gauss-Seidel, setting w = 1.

SORtab =Table[ SORInverse[D2mat,rhsvec,1.0,jj],{jj,200,1000,200}];

Gtab = Table[PlotPoisson[SORtab[[jj]],Nx,dx,Ny,dy],{jj,1,Length[SORtab]}];

Ga =PlotPoisson[fxy,Nx,dx,Ny,dy]

Grout = GraphicsGrid[{{Gtab[[1]],Gtab[[2]],Gtab[[3]]},{Gtab[[4]],Gtab[[5]],Ga}}]

It is interesting to see what happens if we increase the parameter w to, say, 1.5:

SORtab =Table[ SORInverse[D2mat,rhsvec,1.5,jj],{jj,200,1000,200}];
Gtab = Table[PlotPoisson[SORtab[[jj]],Nx,dx,Ny,dy],{jj,1,Length[SORtab]}];
Ga =PlotPoisson[fxy,Nx,dx,Ny,dy];
Grout = GraphicsGrid[{{Gtab[[1]],Gtab[[2]],Gtab[[3]]},{Gtab[[4]],Gtab[[5]],Ga}}]

Just visually (which is not enough), it appears that the convergence happens much faster using a larger SOR parameter (you should try using different values, in a more quantitative setting - one can set w too big, of course).

Lanczos Iteration

Lanczos iteration constructs the basis vectors forming the columns of Q, and the symmetric tridiagonal matrix T in (11.32) - this is an implementation of Algorithm 11.1:

Lanczos[inA_, inb_,m_] := Module[{retQ,retT,alpha,beta,n,w,v,index},
n = Length[inA];
retQ = Table[Table[0.0,{i,1,n}],{j,1,m+2}];
alpha = Table[0.0,{i,1,m+1}];
beta = Table[0.0,{j,1,m+1}];
retQ[[2]] = inb/Sqrt[inb.inb];
For[index = 2, index ≤  m+1, index=index+1,
w = inA.retQ[[index]];
alpha[[index]] = retQ[[index]].w;
v = w - beta[[index-1]] retQ[[index-1]] - alpha[[index]] retQ[[index]];
retQ[[index+1]] = v/Sqrt[v.v];
beta[[index]] = retQ[[index+1]].v;
If[Abs[beta[[index]]] ≤ 10^(-10), Print["Warning - try another vector b"];];
];
retQ = retQ[[2;;m+1]];
beta = beta[[2;;m]];
alpha = alpha[[2;;m+1]];
retT = DiagonalMatrix[alpha] + DiagonalMatrix[beta,1] + DiagonalMatrix[beta,-1];
 Return[{Transpose[retQ],retT}]; 
]

We’ll test on a symmetric, positive definite matrix (generated, by hand):

A = {{1,2,3},{4,5,6},{7,8,10}};
A = Transpose[A].A;
b = {1,1,1};

Here we get back Q and T for m = 3:

{Qm,Tm} = Lanczos[A,b,3];

The decomposition says that A = Q T Q^T:

Max[Abs[Qm.Tm.Transpose[Qm] - A]]

so that checks.  Since A is 3x3, the m=3 Krylov subspace should cover the whole space, and the eigenvalues of T should be identical to those of A (as described in Section 11.3.3)

Eigenvalues[Tm]

N[Eigenvalues[A]]

If, instead, we use a Krylov subspace smaller than the full space, then the eigenvalues of T (which are then the Ritz values of A) will only approximate the eigenvalues of A:

{Qm,Tm} = Lanczos[A,b,2];

Eigenvalues[Tm]

N[Eigenvalues[A]]

Example - Figure 11.3

This is the 10x10 example set up in Section 11.3.3.  We make a symmetric, positive definite matrix A, and then calculate the eigenvalues of T, using a Krylov subspace of dimension m = 2 -> 10.

SeedRandom[111111];
A = Table[RandomReal[{-1,1}],{i,1,10},{j,1,10}];
A = Transpose[A].A;
evals = Eigenvalues[A]

b = {1,1,1,1,1,1,1,1,1,1};

The function below will actually perform the calculations:

Evals[inA_, inb_, max_] := Module[{Qm,Tm,index,outes},
outes = {};
For[index = 2, index ≤ max, index=index+1,
{Qm,Tm} = Lanczos[A,b,index];
outes = Join[outes,{Eigenvalues[Tm]}];
];
outes = Table[Table[{j+1,outes[[j,i]]},{i,1,Length[outes[[j]]]}],{j,1,Length[outes]}];
Return[outes];
]

Run the test, and plot the result:

plotL = Evals[A,b,10];
exes = Table[{j+1,evals[[k]]},{k,1,Length[evals]},{j,1,Length[evals]-1}];

G0 = ListPlot[plotL,PlotMarkers-> {Automatic},PlotStyle->GrayLevel[.5]];
G1 = ListLinePlot[exes,PlotStyle->{GrayLevel[0.],GrayLevel[0.1],GrayLevel[0.2], GrayLevel[0.3],GrayLevel[0.4],GrayLevel[0.5],GrayLevel[0.6],GrayLevel[0.7],GrayLevel[0.8],GrayLevel[0.9]}];
Show[G0,G1,PlotRange-> {{1,10},{-1,12}},AxesOrigin->{1,-1}]

Example - infinite square well

We can apply the Lanczos iteration to find the approximate eigenvalues from any of our quantum mechanical problems.  Here’s the Schrodinger finite difference-maker from Chapter 10:

FDSchrodingerMat[V_, Nn_, dq_] := Module[{retmat,index},
retmat = Table[0.0,{i,1,Nn},{j,1,Nn}];
retmat[[1,1]] = -2.0 - dq^2 V[dq];
retmat[[1,2]] = 1.0;
For[index = 2, index ≤ Nn-1, index=index+1,
retmat[[index,index-1]] = 1.0;
retmat[[index,index]] =-2.0 - dq^2 V[index dq];
retmat[[index,index+1]] = 1.0;
];
retmat[[Nn,Nn-1]] = 1.0;
retmat[[Nn,Nn]] = -2.0 - dq^2 V[Nn dq];
retmat = -1.0/dq^2 retmat;
Return[retmat];
]

We’ll work with the infinite square well, for testing

V[x_] := 0.

Here’s the matrix, and an arbitrary b vector for forming the subspace

A = FDSchrodingerMat[V,100,1./(101)];
b= Table[j,{j,1,100}];

Now we’ll make T for m = 80, just to see how the eigenvalues line up with the known solution:

{Qm,Tm} = Lanczos[A,b,80];

Sort[Eigenvalues[Tm]]

Sort[Eigenvalues[A]]

the rule of thumb is that the largest and smallest eigenvalues of T will converge to the smallest and largest eigenvalues of A “faster”, so we expect the extremes of the two lists above to match more than the intermediate values.

Conjugate Gradient - Figure 11.4

Here is the conjugate gradient implementation from Algorithm 11.2.  Instead of the epsilon cut-off shown there, we use the dimension of the subspace to define the stopping point, this is just for testing, so we can track the accuracy improvement that comes from increasing the dimension of the subspace.

CG[inA_,inb_, maxm_] := Module[{alpha, beta, pm, xm, rm,rmm1,n,m},
n = Length[inA];
xm = Table[0.0,{i,1,n}];
rm = inb;
pm = inb;
For[m = 1, m ≤ maxm, m = m+1,
alpha = rm.rm/(pm.inA.pm);
xm = xm + alpha pm;
rmm1 = rm;
rm = rm - alpha inA.pm;
beta = rm.rm/(rmm1.rmm1);
pm = rm + beta pm;
];
Return[xm];
]

The next block makes the strongly diagonal test matrix described in the example from Section 11.3.5:

NN = 500;
SeedRandom[11111];
B = Table[RandomReal[{-1,1}],{i,1,NN},{j,1,NN}];
A = Transpose[B].B;
tau = 0.2;
A= Table[If[Abs[A[[i,j]]] > tau, 0.0,A[[i,j]]],{i,1,NN},{j,1,NN}] + DiagonalMatrix[Table[1.0,{i,1,NN}]];
b = Table[RandomReal[{-1,1}],{i,1,NN}];
MatrixPlot[A,ColorFunction->"GrayTones"]

We make a table of the log of the residual of the approximate solution for Krylov subspace dimension m = 1 -> 40

conv = Table[{mm,Log[10,Norm[A.CG[A,b,mm]-b]]},{mm,1,40}];
ListPlot[conv,PlotRange->All, PlotStyle->GrayLevel[0.5]]

so that it is clear, even at m = 40 (out of a 500 dimensional complete space), the residual error has norm on the order of 10^(-10).