General Project Objectives: Often in Industry, solutions are developed for small scale test problems using some high level computational system like MATLAB and later solutions for larger production problems are developed using an intermediate level programming language like C or perhaps Fortran so that the final product is more efficient for large scale production problems, i.e., problems that are very large with many computing runs.

Problem Options:

1. MATLAB Small Scale & Programming Large Scale: Hence, develop a procedure for the following project using MATLAB and then use the MATLAB procedure to develop a larger scale programming language (like in C or other scientific language such as FORTRAN, but MATLAB may also be used for the large scale problem too) procedure.

2. Entirely in MATLAB: Alternatively, you can do the problem entirely in MATLAB by finding the largest size problem you can run on the system of your choice based upon the initial test problem, but give justification for using the largest size.

3. Entirely in Ch: Another alternative is to use the Ch Language System which combines the benefits of Unix, C and MATLAB in one self-contained package for MSwindows and other operating systems and is free for academic use.
   Download Ch
   (Developed by Prof. Harry Cheng of UC/Davis; he is a UIC Ph.D. in ME, Prof. K. Gupta was his primary advisor and Hanson was his secondary advisor).

Project Statement: An industrial chemical diffusion reactor has been modeled as a partial differential equations, reduced to a linear algebraic problem, a*Z = b, bZ finite differences and requires a numerical solution in terms of a LU decomposition solution technique since the solution procedure may be performed multiple times for different parameter values that do not affect the

• Coefficient Matrix:   a = [a(i,j)]_n×n.

• Discretized Physical Coordinates (meters): (x1,x2) = (i1*h1, i2*h2) for i1=1 to n1 and i2=1 to n2.
  For MATLAB test use n1=4 and n2=3 interior points per dimension, respectively; but for programming procedure use n1=40 and n2=30 or for pure MATLAB project use the largest since that roughly scales the same way.

• Side Physical Lengths (meters): L1 = 7.891; L2 = 10.156;

• Step Sizes (meters): h1 = L1/(n1+1); h2 = L2/(n2+1);

• Boundary Conditions: u = 0 on boundaries at i1 = 0 or i1 = n1+1 or i2 = 0 or i2 = n2+1, but also at i1 = n1+1 and i2 = n2+1 for upwinding.

• Constant Coefficients: c0 = -0.001758 (Note c0 < 0); c1 = +0.01189; c2 = +1.392; c3 = +1.284;

• Side Indices for 2D reactor geometry u(i1,i2): i1=1:n1; i2=1:n2;     

• Finite Difference: Central finite differences are generally used for linear algebraic vector form Z(i) = u(i1,i2), where the elemental form of the equation is

  sum_j=1:na(i,j)*Z(j) = b(j), for i=1:n;

  but forward upwinding in drift direction for stability.

• Right-Hand-Side Vector: B = [b(i)] with n = n1*n2 components

  for i1 = 1:n1, for i2 = 1:n2, i = i1+n1*(i2-1);
  b(i) = b(i1+n1*(i2-1)) = c0*(i1+i2)/(n1+n2); end, end

• Coefficient Matrix: a = [a(i,j)]_n×n, where n = n1*n2
  Elements in Composite Form:

  • Initialization:           a = zeros(n,n); % Handle sparseness.

  1. For Z(i) = u(i2,i2), diagonal (CFD) terms:
     for i = 1:n, a(i,i) = -2*c1*(1/(h1*h1)+1/(h2*h2)); end

     Corrections for forward (FFD) upwinding:
     for i1=1:n1-1, for i2=1:n2, i=i1+n1*(i2-1); a(i,i) = a(i,i) - 1.5*c2/h1; end, end
     for i1=1:n1, for i2=1:n2-1, i=i1+n1*(i2-1); a(i,i) = a(i,i) - 1.5*c3/h2; end, end

  2. For Z(i-1) = u(i1-1,i2), super-diagonal (CFD) terms :
     for i1=2:n1, for i2=1:n2, i=i1+(i2-1)*n1; j=i-1; a(i,j) = c1/(h1*h1); end, end

     Corrections for near boundary (CFD) (BC: u(n1+1,i2)=0;):
     i1 = n1; for i2=1:n2, i=i1+(i2-1)*n1; j=i-1; a(i,j) = a(i,j) - 0.5*c2/h1; end

  3. For Z(i+1) = u(i1+1,i2), sub-diagonal terms with (FFD) upwinding (u(n1,i2)=0;):
     for i1=1:n1-1, for i2=1:n2, i=i1+(i2-1)*n1; j=i+1; a(i,j) = c1/(h1*h1) + 2*c2/h1; end, end

  4. For Z(i+2) = u(i1+2,i2), 2nd sub-diagonal terms with (FFD) upwinding (u(n1+1:n1+2,i2)=0;):
     for i1=1:n1-2, for i2=1:n2, i=i1+n1*(i2-1); j=i+2; i1=1:n1-2, a(i,j) = - 0.5*c2/h1; end, end

  5. For Z(i-n1) = u(i1,i2-1), n1th super-diagonal terms (u(i1,0)=0;):
     for i1=1:n1, for i2=2:n2, i=i1+(i2-1)*n1; j=i-n1; a(i,j) = c1/(h2*h2); end, end

     Corrections for near boundary CFD:
     i2 = n2; for i1=1:n1, i=i1+(i2-1)*n1; j=i-n1; a(i,j) = a(i,j) - 0.5*c3/h2; end

  6. For Z(i+n1) = u(i1,i2+1), n1th super-diagonal terms with FFD upwinding (u(i1,n2+1)=0;):
     for i1=1:n1, for i2=1:n2-1, i=i1+(i2-1)*n1; j=i+n1; a(i,j)=c1/(h2*h2)+2*c3/h2; end, end

  7. For Z(i+2*n1)=u(i1,i2+2), 2*n1th super-diagonal terms with FFD upwinding (u(i1,n2+1:n2+2)=0;):
     for i1=1:n1, for i2=1:n2-2, i=i1+n1*(i2-1); j=i+2*n1; a(i,j) = -0.5*c3/h2; end, end

General Instructions:

1. Small Scale MATLAB Test Case Part: In MATLAB (see help resources below), use LU decomposition function "[L,U,P] = lu(a)" where "P" is the pivot matrix such that "P*a = L*U" approximately, and the forward with backward solve steps like "Y₁ = L\P*b" and "Z₁ = U\Y₁" (backslash alert, see "help \"), to get an approximate solution "Z₁" for "Z" to the smaller test problem "a*Z = b". Then, use one iterative refinement with the residual "r₁ = b - a*Z₁", and another solve for the approximate error solution "e₁" for "e" in "L*U*e=P*r₁" using the same LU-decomposition. Get the second approximation to the exact solution "Z", "Z₂ = Z₁ + e₁", the infinity-norm of the difference "||Z₂ - Z₁||", the infinity-norm of the first relative residual "||r₁||/||b||", and the infinity-norm of the second relative residual "||r₂||/||b||",

2. Small Scale Output: Hand in your printed MATLAB (*.m file) code and output with explanatory comments for documentation. (Code documentation is very important in Industry.) Report your output in 5 significant digits by print formatting using exponential notation, unless you are tabulating a column of numbers that do not contain very small or very large numbers. Errors and changes should always be printed in e-notation for professional industrial output. In this project the approximate solutions, "Z₁" and "Z₂", must be printed in e-notation (use fprintf in both MATLAB (or variants in C and Fortran languages).

3. Large Scale Part:

   • Programming Language: If you are going to use a programming procedure (C, F or equivalent), solve the large scale version of "a*Z = b" using the Numerical Recipes, or similar, LU decomposition program procedure "ludcmp" and the backsolve companion "lubksb" (or equivalent procedures: it is recommended that you modify these codes to remove special NRC conversion header utilities that are only simulations Fortran style for the conversion to C, C++ or F.

   • OR MATLAB at Maximal Array Size: If you are doing a MATLAB large scale version, doing the largest matrix "a" possible on your system, then make sure your code is in maximal array-vector format, since MATLAB is more efficient in the format due to background use of LAPACK.

   • Either Case: Use your LU decomposition for the initial solution "Z₁" as well as the single iterative refinement to compute second approximation to the exact solution "Z", "Z₂ = Z₁ + e₁", and compute the infinity-norm of the difference "||Z₂ = Z₁||₁", the infinity-norm of the difference "||Z₂ - Z₁||", the infinity-norm of the first relative residual "||r₁||/||b||", and the infinity-norm of the second relative residual "||r₂||/||b||",

4. Plot test, and if possible large scale, results for "Z₁", say using MATLAB's surf or mesh 3D plotting functions, on the original rectangular model, i.e., convert the LU solution "Z₁ = [Z₁(i)]_n×1" to rectangular form
   "Urect =[urect(i1,i2)]_n1×n2" as a function of "x1 = [i1*h1]_n1×1" and "x2 = [i2*h2]_n2×1",
   adding the boundary condtions "u = 0" on the rectangular boundaries, "x1=0", x2=0" "x1=L1" and "x2=L2".

5. Hand In: your printed programming code and truncated output with explanatory comments for documentation together with your MATLAB test.

6. Printing Output: In order to conserve printing and paper resources, you need only print out the coefficient matrix "a", right hand side vector "b" and solution "Z" for i = 1 to n in steps of "K = (n-1)/(nml-1)", designed to make large scale values comparable to small scale values, and the same for the column index "j" of "a", where "nml = 3*4" is the MATLAB (ml) array size. You will have to fold the printing of the rows of a in half, at least.

7. Computer Precision: In your code calculations, you must use at least double precision (64 bit; MATLAB is double by default), since float or single is not suitable for "non-toy" problems. Report your output in 5 significant digits by print formatting using exponential notation, unless you are tabulating a column of numbers that do not contain very small or very large numbers. Errors and changes should always be printed in e-notation for professional output. The solution "Z" must be printed in e-notation.

Project Report: Your report needs the following parts:

1. Cover Page: Put a project title, your name, your affiliation, date and other identifying information on this individual computer project. What you submit must be your own work (this in NOT a group project), else points will be deducted for similar work, i.e. your grade will receive a discount.

2. Executive Summary: This is about a page briefly summarizing the project and results for a busy boss . This should be in the form of a page outline or itemized list for easy and fast reading. This summary is not the same thing as an abstract. To give the boss a taste of the results, show a graph of the most important results.

3. Project Description or Introduction: Describe the project in your own words as an introduction to your report, in sufficient depth so that a reader such as yourself would understand it.

4. Methods: Describe the methods used to solve the problem, giving both advantages and disadvantages in a fair manner. Describe the theory and model behind the computer problem.

5. Results: Describe the nature of the results and illustrate them with appropriate tables or plots. You can use MATLAB for plotting both your MATLAB and C or equivalent results, importing the latter's output date into MATLAB. Clearly label tables and plot figures. Also give the maximum and minimum of "u" in both small and large scale cases.

6. Discussion of Results: Discuss the results, including how they can be used elsewhere for different industrial applications.

7. Acknowledgements: Acknowledge what resources you used in this project, including what versions of MATLAB and C or equivalent that you used, the operating system, the computer or hardware platform, persons consulted (this is important to help avoid deductions for plagerism), and any other resources (references are listed in the next section) used.

8. Conclusions: List what you have learned from this project and explain why it is significant.

9. References: Cite all books, scientific papers, websites and other library or web resources that you used. Give author, title, journal name or book publisher or URL where appropriate, and date of publication or web access.

10. Appendices: Include MATLAB and C or equivalent documented source code and output.

Project Resources:

1. UIC PCLabs MATLAB available software. Also check out your departmental or laboratory computers such as in EECS and MSCS.

2. MATLAB Help Page (Hanson).

3. (NRC): William H. Press et al., Numerical Recipes in C, C++ or F: The Art of Scientific Computing, CambridgeU. Press, ISBN: 0521431085 (See UIC TextBookStore or www.amazon.com price: $57.95 last year; also available for F77 and/or F90; Disk of programs available in Fortran and nonstandard fortran translated C); (See www.amazon.com or similar online bookstores). Also, available at
   
   
   • UIC Libraries 430 SEO:
     
     
     • QA76.73.C15 N865 1992 Main Library for C;
       
     • QA297.N866 1992 Math Library 430 SEO for Fortran;
   
   
   • Online for NRC or NRF or NRF90:
     
     
     • NR Texts:
       Click on Numerical Recipes Home Page
       OR
       BOOKS On-Line link.

       To Get Postscript Copy of pages from "Numerical Recipes": C, F and F90. Unfortunately, C++ version is not available due to copyright restrictions.
       OR
       NR Forum, the NR news list.
       

     • NR Software: Click on NR Software Store On-Line Link.
       You have to buy a single-site license on-line if you do not buy a tex containing software to use NR software.

4. Google Numerical Recipies or MATLAB. Many students have "googled" online to find NR information and codes. For example, search for "ludcmp.c" or "ludcmp.f".

5. C/C++ Help Page (Hanson).

6. Fortran Help Page (Hanson).

7. Ch Homepage, a combined Unix, C and MATLAB like package, so no additional software, no C compiler, no Unix system or no MATLAB, are needed, even in MSwindows. The primary interface is a Unix-like terminal window in which commands can be issued, including the execution of C programs. The Ch command names are typically longer than MATLAB names.
   
   Web Based Ch Computing, a simple graphic interface Ch online tester.