MCS 571 Computer Project 2 -- Numerical EPDE Methods -- Spring 2004

MCS 571 Computer Project 2 -- Numerical EPDE -- Spring 2004

Comparison of Numerical Methods for Linear Elliptic PDEs:
MATLAB or C/Fortran Programming Computation

MATLAB or C/Fortran Programming Output are Due Wednesday 28 April 2004 in Class

General Method Specifications:

Write a single, efficient, general program, in MATLAB, Fortran, C that compares the performance of the three numerical methods:

Successive Overrelaxation (SOR),
Alternating Directions Implicit (ADI),
Conjugate Gradient Method (CGM)

for Elliptic PDEs. (See Prof. Hanson for any instruction variations.)

Apply these methods to the inhomogeneous, linear Elliptic PDE (LEPDE):

0 = F(x,y,u,u_x,u_y,u_xx,u_yy) = A₁*u_xx + A₂*u_yy + B₁(x,y)*u_x + B₂(x,y)*u_y +C(x,y)*u + D(x,y);

where a₁ < x < b₁; a₂ < y < b₂.

Caveat Usor: You are responsible for checking the correctness of the formulas, given here and for any modifications needed to transform them to useful code.

Detailed Method Specifications:

The SOR and CGM methods are all derived from the Euler's explicit method which in the linear algebra formulation (with arrays marked with an M for matrix to avoid confusion with the model) is

AM*XM = BM; XM = {U_i,j | interior, row order: for j = 1 : m₂, i = 1 : m₁};

XM_k = XM^(k); XM^(k)_i+(j-1)*m₁ = U_i,j^(k); AM = LM+UM-DM; BM = {RHS vector, including boundary conditions};

h_q = (b_q - a_q)/(m_q+1); for q = 1 : 2,

where m₁ and m₂ are the total number of interior points for x and y, respectively. Let k_max=75 be the maximum number of iterations and tol = 0.5e-5. Caution: UM here means upper triangular part of AM, while U_i,j^(k) denotes the kth iteration of a grid algebra component at node (i,j). (Remark: If you use MATLAB, then there may be three sets of indices to keep track of:

original Grid Algebra (i,j) = (0:m₁+1,(0:m₂+1;
base one MATLAB (im,jm) = (i+1,j+1) ;
Linear Algebra il = i+m1*(j-1) and similarly jl;

1. Successive Over Relaxation (SOR) Method:
Let XM₀ = -w*DM^-1*BM to start.
For k = 0 : k_max While (||RM_k||₂ > tol*||BM||₂) Do:

XM_k+1 = XM_k - w*DM^-1*RM_k;            %%SOR
where the SOR residual is
RM_k = BM-LM*XM_k+1-(UM-DM)*XM_k ;
which is not the same as the standard residual, and w=1.25 is the starting SOR acceleration parameter. Hence,
(DM-w*LM)*XM_k+1 = ((1-w)*DM + w*UM)*XM_k-w*BM ;
in linear algebra iteration form. The grid algebra form can also be used for SOR,
U_i,j^(k+1) = c_0,0,i,j*U_i,j^(k) + c_-1,0,i,j*U_i-1,j^(k+1) + c_0,-1,i,j*U_i,j-1^(k+1) + c_1,0,i,j*U_i+1,j^(k) + c_0,1,i,j*U_i,j+1^(k) + c_0,i,j ,
where the coefficients c_p,q,i,j for (p, q) = (-1:+1, -1:+1) and c_0,i,j depend on the PDE model and the FD scheme.

2. Alternating Directions Implicit (ADI) Method:
The ADI method is based on operator splitting by dimension for a higher order approximation to replace the Crank-Nicolson Implicit Method CNIM in higher dimension to preserve the tridiagonal advantage in using the Thomas algorithm (See your class notes or Morton and Mayer, the second opinion text). The method specified here is a modification of the PRADI variant of ADI for the variable coefficients and extra drift-source terms here. Use the same initial condition X₀ as for the SOR method. Let the ADI acceleration parameter be symbolically given by
p = padi = 0.5*Delta_t*Adh2max; where Adh2max = max_i[A_i/h_i²] ;
where Delta_t is an artificial time step from the CNI PPDE derivation, the relative diffusion coefficient is
R_i = (A_i/h_i²)/Adh2max;
delta_i² is the CFD operator for a 2nd order derivative and Delta_i is the CFD operator for a first order derivative, for i = 1:2. Let p = padi = w/Adh2max; to starting. The PRADI variant algorithm for our LEPDE here would be,
ka_max = 0.5*k_max;
for k=0 : ka_max While (||U^(k+1)-U^(k)||₂ > tol*||U^(k)||₂) Do:
%%% X-ADI:
(1-p*R₁*(delta₁² + 0.5*(B_1,i,j/A_1,i,j)*h₁*Delta₁ + 0.5*(C_i,j/A_1,i,j)*h₁²*I)) [U_i,j^(k+0.5)] =(1+p*R₂*(delta₂² + 0.5*(B_2,i,j/A_2,i,j)*h₂*Delta₂ + 0.5*(C_i,j/A_2,i,j)*h₂²*I)) [U_i,j^(k)] + 2*p*D_i,j/Adh2max ;

%%% Y-ADI:
(1-p*R₂*(delta₂² + 0.5*(B_2,i,j/A_2,i,j)*h₂*Delta₂ + 0.5*(C_i,j/A_2,i,j)*h₂²*I)) [U_i,j^(k+1)] =(1+p*R₁*(delta₁² + 0.5*(B_1,i,j/A_1,i,j)*h₁*Delta₁ + 0.5*(C_i,j/A_1,i,j)*h₁²*I)) [U_i,j^(k)] + 2*p*D_i,j/Adh2max ;
These must be converted to a suitable form for the Thomas algorithm.

3. Conjugate Gradient Method (CGM):
The CGM method is best modified to account for non-symmetric matrices A and for program efficiency.
For this method, use one of the following professional codes given in MATLAB from the NeTLib software archive (the class pseudo-code only gives the major ideas):

Conjugate Gradiant Squared Method (CGSM), see "help cgs" in MATLAB or use a text editor for the lead comments specifying the meaning of the input and output arguments,
function [x, error, iter, flag] = cgs(A, x, b, M, max_it, tol)
Bi-Conjugate Gradiant Method (BCGM), see "help cgs" in MATLAB or use a text editor for the lead comments specifying the meaning of the input and output arguments,
function [x, error, iter, flag] = bicg(A, x, b, M, max_it, tol)
For each of these method M is the preconditioner or starting approximate for A = AM which can be taken as M = -DM.

Test Case Specifications:
Your program should be fairly general and not made just for the particular values or functions here, but you can use the following values as a test case as long as your program is not specifically tailored to it. where a₁ < x < b₁, a₂ < y < b₂. On
(a₁,b₁) ×(a₂,b₂) = (exp(-1),exp(1))×(-PI,PI) ;
where PI = 4*atan(1), and the spatially variable coefficients are
A₁ = exp(1.5);    A₂ = exp(1.25);
B₁(x,y) = exp(-1.3)*(x-a₁)/sqrt(1+x²+y²);    B₂(x,y)=exp(-1.1)*(y-a₂)/sqrt(1+x²+y²);
C(x,y) = -exp(1.5)/sqrt(1+x²+y²);    D(x,y) = -exp(-2.5)/sqrt(1+x²+y²);
with boundary conditions,
u(x,a₂) = 6*exp(0.5)*(x-a₁)/(b₁-a₁) ;
u(a₁,y) = 3*exp(0.5)*(y-a₂)/(b₂-a₂) ;
u(x,b₂) = exp(0.5)*(3+(x-a₁)/(b₁-a₁)) ;
u(b₁,y) = 2*exp(0.5)*(3-(y-a₂)/(b₂-a₂)) ;
The mesh requirements are that m₁ = 6 and m₂ = 5. Let tol = 0.5e-6 and kmax = 10 iterates.

Output Specifications:

In a well-labeled rectangular table, tabulate the SOR, ADI and CGM
U_i,j^(k)
for (i,j) = (0,0) : (m₁+1,m₂+1) at the starting iterate, for the first two new iterates (for SOR and ADI only) and the last full iterate, to 5 significant digits in U_i,j^(k). Also make a plot of the final iterate. Print out the final value (successful or not) of the stopping criterion and the total number of iterations. For CGM methods also print out the estimated error. Note that you have to convert XM to U, including boundary conditions for CGM and other methods where the linear algebra form was used.
Report the timings for each method of solution. You can use clock with etime or btter timer in MATLAB or C or Fortran. Specify what you are actually timing and make sure it is a fair comparison, being sure to mention if you have excluded a common setup of the problem prior to solving.
Comment on the nature of the difference among SOR, ADI and CGM and say why. State what computer and what compiler were used.

Efficiency Hints:

Compute once and store the spatially dependent coefficients in an array to use in all iteration steps.
Avoid calculating items like constant expressions in loops.

Web Source: http://www.math.uic.edu/~hanson/mcs571/ecp571s04:.html
MCS 571 HomePage: http://www.math.uic.edu/~hanson/mcs571/
Email Comments or Questions to hanson A T uic edu

MCS 571 Computer Project 2 -- Numerical EPDE -- Spring 2004

Comparison of Numerical Methods for Linear Elliptic PDEs: MATLAB or C/Fortran Programming Computation

Comparison of Numerical Methods for Linear Elliptic PDEs:
MATLAB or C/Fortran Programming Computation