Multigrid method

1D elliptic problem

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

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]

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;ua=1;ub=3;c=1/2;
ϕ[x_]=20π (x^3);
f[x_]=-20+c(∂_(x,x))ϕ[x]Cos[ϕ[x]]-c((((∂_x)ϕ[x]))^2)Sin[ϕ[x]];
sol=Delete[DSolve[{(∂_(x,x))U[x]==f[x],U[a]==ua,U[b] ==ub},U[x],x],0]; (* Delete[...,0] deletes outer {} *)
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])=a+j h  where h=(b-a)/(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])=a)=(Subscript[u, a]),u((Subscript[x, n+1])=b)=(Subscript[u, b])

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

Evaluating the function f(x) for many (Subscript[x, j]) takes a long time in Mathematica therefore we calculate it in advance and use the vector fx in all calculations later. We also calculate the exact solution on the grid as ux for later use.

ig=10; (* the highest level used in multigrid below *)
n=(2^ig)-1;h=(b-a)/(n+1);
Print["Highest number of points: n = ",n]
X=Range[a,b,h];
Xin=X[[2;;n+1]];
(* fx will be a packed array (machine precission) thanks to N[] *)
fx=N[f[Xin]]; 
ux=N[u[X]];

The right-hand-side will be -(h^2)f((Subscript[x, j])) but the first and last element must be modified using Dirichlet's boundary conditions u((Subscript[x, 0])) for

-(((d^2)u((Subscript[x, 1])))/(d(x^2)))≃((-u((Subscript[x, 0]))+2u((Subscript[x, 1]))-u((Subscript[x, 2])))/(h^2))=-f((Subscript[x, 1]))  ⟹ 2u((Subscript[x, 1]))-u((Subscript[x, 2])) = -(h^2)f((Subscript[x, 1]))+(Subscript[u, a])

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

rhs=-(h^2)fx;
rhs[[1]]=rhs[[1]]+ua;
rhs[[n]]=rhs[[n]]+ub;

Now we are ready to solve the problem on the grid:

T=SparseArray[{{i_,i_}->2.0,{i_,j_}/;Abs[i-j]==1->-1.0},{n,n}];
xDirect=LinearSolve[T,rhs];

Print["Maximum error = ",Max[Abs[ux[[2;;n+1]]-xDirect]]]
exactSol=Transpose[{X,ux}]; (* for later use in plots *)
ListPlot[{exactSol,Transpose[{Xin,xDirect}]},Joined->{True,False},AxesLabel->{"x","f(x)"}]

Multigrid solution - using only the same V-cycle

Multigrid V-cycle algorithm (MGV) starting on the level (i) consists of these steps:

improving the solution (here using the weighted Jacobi method with ω=2/3)

(x^((i)))=Smoothing((b^((i))),(x^((i))))

computing the residual

(r^((i)))=(T^((i)))(x^((i)))-(b^((i)))

solving the problem

(T^((i)))(d^((i)))=(r^((i)))

by recursive restriction on the coarser grids (4 is there because for the coarser grid we get (((2h))^2) on the right-hand-side)

(d^((i)))=Interpolation(MGV(4Restriction((r^((i)))),0))

correcting the fine grid solution

(x^((i)))=(x^((i)))-(d^((i)))

improving the solution again (usually using the weighted Jacobi method with ω=2/3)

(x^((i)))=Smoothing((b^((i))),(x^((i))))

If (i) is the lowest level (1) then we solve the problem directly.

The following functions do the particular jobs described above:

smoothing[b_,x0_,ω_,niter_]:=
Module[{n,xNew,R,c,ω2},
n=Length[x0];
xNew=x0;
ω2=0.5ω;
R=SparseArray[{{i_,i_}->(1.0-ω),{i_,j_}/;Abs[i-j]==1->ω2},{n,n}];
c=ω2 b;
Do[
xNew=R.xNew+c,
{i,1,niter}
];
Return[xNew]
]
restriction[r_]:=
Module[{n,P},
n=Length[r];
P=SparseArray[{{i_,j_}/;2i-1==j->0.25,{i_,j_}/;2i==j->0.5,{i_,j_}/;2i+1==j->0.25},{(n-1)/2,n}];
Return[P.r]
]
interpolation[r_]:=
Module[{n,P},
n=Length[r];
P=SparseArray[{{i_,j_}/;2j-1==i->0.5,{i_,j_}/;2j==i->1.0,{i_,j_}/;2j+1==i->0.5},{2n+1,n}];
Return[P.r]
]
MGV[b_,u_]:=
Module[{n,uNew,T,r,d},
n=Length[b];
If[n==1,
(* direct solution *)
uNew=0.5b,
(* iterations *)
uNew=smoothing[b,u,2.0/3.0,1];
T=SparseArray[{{i_,i_}->2.0,{i_,j_}/;Abs[i-j]==1->-1.0},{n,n}];
r=T.uNew-b;
d=interpolation[MGV[4.0 restriction[r],ConstantArray[0.0,(n-1)/2]]];
uNew=uNew-d;
uNew=smoothing[b,uNew,2.0/3.0,2];
];
Return[uNew];
]

To test one V-cycle we use as an initial guess a straight line satisfying the boundary conditions:

u0[x_]=InterpolatingPolynomial[{{a,ua},{b,ub}},x]
x0=N[u0[Xin]];

Let us apply the V-cycle recursively several times and see how the error is decreased. We compare the iterative solution with the direct one obtained above:

niter=10;
xMGV=x0;
errorMGV=ConstantArray[1.0 (10^(-20)),niter+1];
errorMGV[[1]]=Max[Abs[xMGV-xDirect]];
Do[
xMGV=MGV[rhs,xMGV];
errorMGV[[i+1]]=Max[Abs[xMGV-xDirect]],
{i,1,niter}
]
eMGV=xMGV-xDirect;

ListLogPlot[{errorMGV},PlotRange->All,AxesLabel->{"niter","error"}]
nit1=6;nit2=10;
ρMGV=Exp[(Log[errorMGV[[nit2]]]-Log[errorMGV[[nit1]]])/(nit2-nit1)];
Print["estimated ρ((Subscript[R, MGV])) = ",ρMGV]
ListPlot[{exactSol,Transpose[{Xin,xMGV}]},Joined->{True,False},AxesLabel->{"x","f(x)"}]
ListPlot[{Transpose[{Xin,eMGV}]},Joined->{False},AxesLabel->{"x","error"}]

Full multigrid solution

The full multigrid starts at the level (1) by solving the problem exactly, then we interpolate the result to the next level and run one V-cycle algorithm (MGV), then we again interpolate to the next level and run another V-cycle etc.

To accomplish this we need a slightly different function for interpolation:

interpolationWithBoundaries[u_,ua_,ub_]:=
Module[{n,P},
n=Length[u];
P=SparseArray[{{i_,j_}/;2j-1==i->0.5,{i_,j_}/;2j-2==i->1.0,{i_,j_}/;2j-3==i->0.5},{2n+1,n+2}];
Return[P.Join[{ua},u,{ub}]]
]
FMG[k_,fx_,a_,b_,ua_,ub_]:=
Module[{n,maxk,i2,h,rhs,u},
n=Length[fx];maxk=Log[2,n+1];
If[k>maxk,Print["Full multigrid can be run only up to the level (",maxk,")"];];
h=(b-a)/(2^maxk);
u={(-fx[[(n+1)/2]]+ua+ub)/2}; (* exact solution at the level (1) *)
Do[
i2=(2^i);
h=(b-a)/i2;
rhs=-(h^2)fx[[Range[(n+1)/i2,n,(n+1)/i2]]];
rhs[[1]]=rhs[[1]]+ua;
rhs[[i2-1]]=rhs[[i2-1]]+ub;
u=MGV[rhs,interpolationWithBoundaries[u,ua,ub]],
{i,2,Min[k,maxk]}
];
Return[u];
]

errorFMG=ConstantArray[1.0 (10^(-20)),ig];
timeFMG=ConstantArray[0,ig];
Do[
t=Timing[xFMG=FMG[i,fx,a,b,ua,ub];];
timeFMG[[i]]=t[[1]];
h=(b-a)/(2^i);
Xin=Range[a+h,b-h,h];
errorFMG[[i]]=Max[Abs[u[Xin]-xFMG]],
{i,1,ig}
];
eFMG=xFMG-xDirect;

ListLogPlot[{timeFMG[[2;;ig]]},PlotRange->All,AxesLabel->{"level","time"}]
ListLogPlot[{errorFMG},PlotRange->All,AxesLabel->{"level","error"}]
ListPlot[{exactSol,Transpose[{Xin,xFMG}]},Joined->{True,False},AxesLabel->{"x","f(x)"}]
ListPlot[{Transpose[{Xin,eFMG}]},Joined->{False},AxesLabel->{"x","error"}]