Iterative methods - 1D elliptic problem

Taken from R.J.LeVeque: Finite Difference Methods for Ordinary and Partial Differential Equations - Steady-State and Time-Dependent Problems, SIAM, Philadephia 2007, chapter 4.6

Preliminaries

Clear all symbols from previous evaluations to avoid conflicts

Clear["Global`*"]

Problem

Differential equation

We would like to solve numerically the differential equation

(((d^2)u(x))/(d(x^2)))=f(x)

with Dirichlet  boundary conditions

u(a)=(Subscript[u, a])
u(b)=(Subscript[u, b])

Particular problem

As the right-hand-side we will take

f(x)=-20+c ϕ''(x)cos(ϕ(x))-((c(ϕ'(x)))^2)sin(ϕ(x))

where

c=1/2
ϕ(x)=20π (x^3)

and boundary conditions are

u(0)=1
u(1)=3

The analytic solution

This problem can be solved in the closed form

a=0;b=1;(Subscript[u, a])=1;(Subscript[u, b])=3;c=1/2;
ϕ[x_]=20π (x^3);
f[x_]=-20+c(∂_(x,x))ϕ[x]Cos[ϕ[x]]-c((((∂_x)ϕ[x]))^2)Sin[ϕ[x]];
(* Delete[...,0] deletes outer {} of a list *)
sol=Delete[DSolve[{(∂_(x,x))U[x]==f[x],U[a]==(Subscript[u, a]),U[b] ==(Subscript[u, b])},U[x],x],0]; 
u[x_]=Expand[U[x]/.sol]
Plot[u[x],{x,a,b},AxesLabel->{"x","f(x)"}]

Direct numerical solution

We discretize the equation () using a standard finite difference formula on the equidistant grid

(Subscript[x, j])=j h  where h=1/(n+1), j=0,...,n+1

i.e. we have to solve the system of n equations

(((d^2)u((Subscript[x, j])))/(d(x^2)))≃((u((Subscript[x, j+1]))-2u((Subscript[x, j]))+u((Subscript[x, j-1])))/(h^2))=f((Subscript[x, j]))   for  j=1,...,n

with boundary conditions

u((Subscript[x, 0])=0)=1,u((Subscript[x, N+1])=1)=3

We can use Mathematica to solve this problem directly. For later convenience in gradient methods we actually solve

-(((d^2)u((Subscript[x, j])))/(d(x^2)))≃((-u((Subscript[x, j+1]))+2u((Subscript[x, j]))-u((Subscript[x, j-1])))/(h^2))=-f((Subscript[x, j]))   for  j=1,...,n

n=100;h=(b-a)/(n+1);
X=Range[a,b,h];Xin=X[[2;;n+1]];
T=SparseArray[{{i_,i_}->2.0,{i_,j_}/;Abs[i-j]==1->-1.0},{n,n}];
rhs=-N[h^2f[X[[2;;n+1]]],16];
rhs[[1]]=rhs[[1]]+(Subscript[u, a]); (* boundary condtions *)
rhs[[n]]=rhs[[n]]+(Subscript[u, b]);
xDirect=LinearSolve[T,rhs];
Print["MaxError = ",Max[Abs[u[X[[2;;n+1]]]-xDirect]]]
exactSol=Transpose[{X,u[X]}]; (* for later use in plots *)
ListPlot[{exactSol,Transpose[{Xin,xDirect}]},Joined->{True,False},PlotLegends->{"Exact","Numerical"},AxesLabel->{"x","f(x)"}]

The matrix is symmetric positive definite, all eigenvalues are positive, all methods should converge

eval=Sort[Eigenvalues[T]];
ListPlot[eval,AxesLabel->{"n","eigenvalues"}]
ListPlot[{eval-(2-2Cos[π Range[n]/(n+1)])},AxesLabel->{"n","eigenvalue error"}]

Basic iterative methods

In general, we solve a system of linear equations

A x=b

and use a standard decomposition

A=D-(OverTilde[L])-(OverTilde[U])=D(I-L-U)

where D  is the diagonal of A,-(OverTilde[L]) and -(OverTilde[U]) are the strictly lower and upper triangular parts respectively and  L=(D^(-1))(OverTilde[L]), U=(D^(-1))(OverTilde[U]),
More generally will write A=M-K.

One step of the general iterative method can be written as

(Subscript[x, m+1])=R (Subscript[x, m])+c  where  R=(M^(-1))K, c=(M^(-1))b

and the method is convergent if and only if the spectral radius of the matrix R

ρ(R) =max|(Subscript[λ, i])|   where  (Subscript[λ, i]) are eigenvalues of R

satisfies

ρ(R) <1.

For our particular problem ()-() we use the same discretization () on the equidistant grid as in the direct method and compare iterative solutions with the solution obtained by the direct method. The standard decomposition is (note that we have A=T and b= rhs and we set D=Diag):

Diag=DiagonalMatrix[Diagonal[T]];
D1=Inverse[Diag];
L=-D1.LowerTriangularize[T,-1];
U=-D1.UpperTriangularize[T,1];
Id=IdentityMatrix[n];
(* MatrixForm[U] *)

An initial guess and the maximum number of iterations will be the same for all methods. As the initial guess we will take the linear function satisfying the boundary conditions ().

u0[x_]=InterpolatingPolynomial[{{a,(Subscript[u, a])},{b,(Subscript[u, b])}},x]
x0=N[u0[X[[2;;n+1]]],16];

And here is a common part of all basic iterative algorithms. It is supposed that (Subscript[n, iter]), (Subscript[x, 0]) and (Subscript[x, direct]) are already assigned.

basicIterativeMethod[matrixR_,vectorC_,niter_]:=
Module[{x,e},
x=x0;
e=ConstantArray[1.0 (10^(-20)),niter+1];
e[[1]]=Max[Abs[x-xDirect]];
Do[
x=matrixR.x+vectorC;
e[[i+1]]=Max[Abs[x-xDirect]],
{i,1,niter}
];
{x,e}
]

Jacobi method

One step of the Jacobi method is

(Subscript[x, m+1,j])=(1/(Subscript[a, jj]))((Subscript[b, j])-(∑_(k≠j))(Subscript[a, jk])(Subscript[x, m,k])) ,j=1,...,n

or in a matrix form

(Subscript[x, m+1])=(Subscript[R, Jac])(Subscript[x, m])+(Subscript[c, Jac])   where (Subscript[R, Jac])=(D^(-1))((OverTilde[L])+(OverTilde[U]))=L+U  and  (Subscript[c, Jac])=(D^(-1))b

Properties of (Subscript[R, Jac]) for the problem ()-():

RJac=L+U;
Print["Part 5x5 of (Subscript[R, Jac]):"]
MatrixForm[RJac[[1;;5,1;;5]]]
cJac=D1.rhs;
(* MatrixForm[R]*)
eigenJac=Eigenvalues[RJac];
ρJac=Max[Abs[eigenJac]];
Print["Spectral radius ρ((Subscript[R, Jac])) = ",ρJac]
ListPlot[Sort[eigenJac],PlotLegends->{"Eigenvalues of (Subscript[R, Jac])"},AxesLabel->{"n","eigenvalues"}]

The main algorithm of the Jacobi method for our problem ()-(). In each step we save the maximum error of the solution.

niter=1000;
{xIter,errorJac}=basicIterativeMethod[RJac,cJac,niter];
xJac=Transpose[{Xin,xIter}]; (* for later use in plots *)
(* Plot results *)
ListLogPlot[{errorJac},PlotRange->All,PlotStyle->{Blue},PlotLegends->{"Jacobi"},AxesLabel->{"niter","max error"}]
ListPlot[{exactSol,xJac},Joined->{True,False},PlotStyle->{Black,Blue},PlotLegends->{"Exact","Jacobi"},AxesLabel->{"x","f(x)"}]

Gauss-Seidel method

One step of the Gauss-Seidel method is

(Subscript[x, m+1,j])=(1/(Subscript[a, jj]))((Subscript[b, j])-(∑_(k=1)^(j-1))(Subscript[a, jk])(Subscript[x, m+1,k])-(∑_(k=j+1)^n)(Subscript[a, jk])(Subscript[x, m,k])) ,j=1,...,n

or in a matrix form

(Subscript[x, m+1])=(Subscript[R, GS])(Subscript[x, m])+(Subscript[c, GS])   where (Subscript[R, GS])=(((D-(OverTilde[L])))^(-1))(OverTilde[U])=(((I-L))^(-1))U  and  (Subscript[c, GS])=(((D-(OverTilde[L])))^(-1))b=(((I-L))^(-1))(D^(-1))b

Properties of (Subscript[R, GS]) for the problem ()-():

ImL1=Inverse[Id-L];
RGS=ImL1.U;
Print["Part 5x5 of (Subscript[R, GS]):"]
MatrixForm[RGS[[1;;5,1;;5]]]
cGS=ImL1.D1.rhs;
(*MatrixForm[R]*)
eigenGS=Eigenvalues[RGS];
ρGS=Max[Abs[eigenGS]];
Print["Spectral radius ρ((Subscript[R, GS])) = ",ρGS]
ListPlot[Sort[Re[eigenGS]],PlotLegends->{"Eigenvalues of (Subscript[R, GS])"},AxesLabel->{"n","eigenvalues"}]

The main algorithm of the Jacobi method for our problem ()-(). In each step we save the maximum error of the solution.

{xIter,errorGS}=basicIterativeMethod[RGS,cGS,niter];
xGS=Transpose[{Xin,xIter}]; (* for later use in plots *)
(* Plot results *)
ListLogPlot[{errorJac,errorGS},PlotRange->All,PlotStyle->{Blue,Red},PlotLegends->{"Jacobi","Gauss-Seidel"},AxesLabel->{"niter","max error"}]
ListPlot[{exactSol,xJac,xGS},Joined->{True,False,False},PlotStyle->{Black,Blue, Red},PlotLegends->{"Exact","Jacobi","Gauss-Seidel"},AxesLabel->{"x","f(x)"}]

Successive overrelaxation method

One step of the SOR(ω) method is

(Subscript[x, m+1,j])=(1-ω)(Subscript[x, m,j])+(ω/(Subscript[a, jj]))((Subscript[b, j])-(∑_(k=1)^(j-1))(Subscript[a, jk])(Subscript[x, m+1,k])-(∑_(k=j+1)^n)(Subscript[a, jk])(Subscript[x, m,k])) ,j=1,...,n

or in a matrix form

(Subscript[x, m+1])=(Subscript[R, SOR])(Subscript[x, m])+(Subscript[c, SOR])   where (Subscript[R, SOR])=(((D-ω(OverTilde[L])))^(-1))[(1-ω)D+ω(OverTilde[U])]=(((I-ω L))^(-1))[(1-ω)I+ω U] 
and  (Subscript[c, GS])=((ω(D-ω(OverTilde[L])))^(-1))b=((ω(I-ω L))^(-1))(D^(-1))b

Properties of (Subscript[R, SOR]) for the problem ()-() for an optimal certain ω:

ω=2/(1+Sqrt[1-((ρJac)^2)]);
Print["ω = ",ω]
ImL1=Inverse[Id-ω L];
RSOR=ImL1.((1-ω)Id+ω U);
Print["Part 5x5 of (Subscript[R, SOR]):"]
MatrixForm[RSOR[[1;;5,1;;5]]]
cSOR=ω ImL1.D1.rhs;
(*MatrixForm[R]*)
eigenSOR=Eigenvalues[RSOR];
ρSOR=Max[Abs[eigenSOR]];
Print["Spectral radius ρ((Subscript[R, SOR])) = ",ρSOR]
ListPlot[Sort[Re[eigenSOR]],PlotLegends->{"Eigenvalues of SOR"}]

The main algorithm of the Jacobi method for our problem ()-(). In each step we save the maximum error of the solution.

{xIter,errorSORopt}=basicIterativeMethod[RSOR,cSOR,niter];
xSORopt=Transpose[{Xin,xIter}]; (* for later use in plots *)
(* Plot results *)
ListLogPlot[{errorJac,errorGS,errorSORopt},PlotRange->All,PlotStyle->{Blue,Red,Green},PlotLegends->{"Jacobi","Gauss-Seidel","optimal SOR"},AxesLabel->{"niter","max error"}]
ListPlot[{exactSol,xJac,xGS,xSORopt},Joined->{True,False,False,False},PlotStyle->{Black,Blue,Red,Green},PlotLegends->{"Exact","Jacobi","Gauss-Seidel","optimal SOR"},AxesLabel->{"x","f(x)"}]

Properties of (Subscript[R, SOR]) for the problem ()-() for an estimated ω:

Print["ρ((Subscript[R, Jac])) = ",ρJac]
(* Determination of ρ((Subscript[R, Jac])) from the error *)
nit1=50;nit2=60;
ρJacEst=Exp[(Log[errorJac[[nit2]]]-Log[errorJac[[nit1]]])/(nit2-nit1)];
Print["estimated ρ((Subscript[R, Jac])) = ",ρJacEst]
(* Use this value in SOR algorithm *)
ω=2/(1+Sqrt[1-((ρJacEst)^2)]);
Print["ω = ",ω]
ImL1=Inverse[Id-ω L];
RSORest=ImL1.((1-ω)Id+ω U);
cSORest=ω ImL1.D1.rhs;
(*MatrixForm[R]*)
eigenSORest=Eigenvalues[RSORest];
ρSORest=Max[Abs[eigenSORest]];
Print["estimated ρ((Subscript[R, SOR])) = ",ρSORest, " (Optimal ρ((Subscript[R, SOR])) = ",ρSOR,")"]
ListPlot[Sort[Re[eigenSORest]],PlotLegends->{"Eigenvalues of estimated SOR"}]

{xIter,errorSORest}=basicIterativeMethod[RSORest,cSORest,niter];
xSORest=Transpose[{Xin,xIter}]; (* for later use in plots *)
(* Plot results *)
ListLogPlot[{errorJac,errorGS,errorSORopt,errorSORest},PlotRange->All,PlotStyle->{Blue,Red,Green,Orange},PlotLegends->{"Jacobi","Gauss-Seidel","optimal SOR","estimated SOR"},AxesLabel->{"niter","max error"}]
ListPlot[{exactSol,xJac,xGS,xSORopt,xSORest},Joined->{True,False,False,False,False},PlotStyle->{Black,Blue,Red,Green,Orange},PlotLegends->{"Exact","Jacobi","Gauss-Seidel","optimal SOR","estimated SOR"},AxesLabel->{"x","f(x)"}]

Gradient iterative methods

Let us solve a system of linear equations

A x=b

with a symmetric positive definite matrix A by searching a minimum of the function

ϕ(x)=(1/2)(x^T)A x-(x^T) b

An initial guess and the maximum number of iterations will be the same as in basic iterative method.

Method of the steepest descent

xIter=x0;
r=rhs-T.x0;
errorMSD=ConstantArray[1.0 (10^(-50)),niter+1];
errorMSD[[1]]=Max[Abs[xIter-xDirect]];
Do[
w=T.r;
α=r.r/r.w;
xIter=xIter+α r;
r=r-α w;
errorMSD[[i+1]]=Max[Abs[xIter-xDirect]],
{i,1,niter}
]
xMSD=Transpose[{Xin,xIter}]; (* for later use in plots *)
(* Plot results *)
ListLogPlot[{errorJac,errorMSD},PlotRange->All,PlotStyle->{Blue,Orange},PlotLegends->{"Jacobi","steepest descent"},AxesLabel->{"niter","max error"}]
ListPlot[{exactSol,xMSD},PlotStyle->{Black,Orange},Joined->{True,False},PlotLegends->{"Exact","steepest descent"},AxesLabel->{"x","f(x)"}]

Conjugate Gradient Method

xIter=x0;
r=rhs-T.x0;
p=r;
γ=r.r;
errorCG=ConstantArray[1.0 (10^(-50)),niter+1];
errorCG[[1]]=Max[Abs[xIter-xDirect]];
Do[
w=T.p;
α=γ/p.w;
xIter=xIter+α p;
r=r-α w;
β=1/γ;
γ=r.r;
β=β γ;
p=r+β p;
errorCG[[i+1]]=Max[Abs[xIter-xDirect]],
{i,1,niter}
]
xCG=Transpose[{Xin,xIter}]; (* for later use in plots *)
(* Plot results *)
ListLogPlot[{errorSORopt,errorMSD,errorCG},PlotRange->All,PlotStyle->{Green,Orange,Blue},PlotLegends->{"optimal SOR","steepest descent","CG"},AxesLabel->{"niter","max error"}]
ListLogPlot[{errorSORopt,errorMSD,errorCG},PlotRange->{{0,200},{0.01,4}},PlotStyle->{Green,Orange,Blue},PlotLegends->{"optimal SOR","steepest descent","CG"},AxesLabel->{"niter","max error"}]
ListPlot[{exactSol,xSORopt,xMSD,xCG},PlotStyle->{Black,Green,Orange,Blue},PlotLegends->{"Exact","optimal SOR","steepest descent","CG"},Joined->{True,False,False,False},AxesLabel->{"x","f(x)"}]

Conjugate Gradient Method with Preconditioning

In this simple case of tridiagonal matrix, it is difficult to choose some reasonable preconditioner. Therefore I only illustrate here the effect by choosing as a preconditioner the original matrix with smaller values of off-diagonal elements.

M=SparseArray[{{i_,i_}->2.0,{i_,j_}/;Abs[i-j]==1->-0.8},{n,n}];
niter=200;
xIter=x0;
r=rhs-T.x0;
errorCGP=ConstantArray[1.0 (10^(-50)),niter+1];
errorCGP[[1]]=Max[Abs[xIter-xDirect]];
Do[
z=LinearSolve[M,r];
If[i==1,
p=z;γ=r.z,
β=1/γ;γ=r.z;β=β γ;p=z+β p
];
w=T.p;
α=γ/p.w;
xIter=xIter+α p;
r=r-α w;
errorCGP[[i+1]]=Max[Abs[xIter-xDirect]],
{i,1,niter}
]
xCGP=Transpose[{Xin,xIter}]; (* for later use in plots *)
(* Plot results *)
ListLogPlot[{errorSORopt,errorMSD,errorCG,errorCGP},PlotRange->All,PlotStyle->{Green,Orange,Blue,Red},PlotRange->All,PlotLegends->{"optimal SOR","steepest descent","CG","CG with preconditioning"},AxesLabel->{"niter","max error"}]
ListLogPlot[{errorSORopt,errorMSD,errorCG,errorCGP},PlotRange->{{0,110},{10^-15,4}},PlotStyle->{Green,Orange,Blue,Red},PlotLegends->{"optimal SOR","steepest descent","CG","CG with preconditioning"},AxesLabel->{"niter","max error"}]
ListPlot[{exactSol,xMSD,xCG,xCGP},PlotStyle->{Black,Orange,Blue,Red},PlotLegends->{"Exact","steepest descent","CG","CG with preconditioning"},Joined->{True,False,False,False},AxesLabel->{"x","f(x)"}]