General Project Objectives:

Many industrial problems, such as chemical or biological rate problems, can be cast as a systems of Linear Ordinary Differential Equations (LODEs), i.e., Vector-Matrix Differential Equations.

dy(t)/dt = a*y(t),    y(t₀) = y₀,    t₀ < t <= t_f ,

Since the problem in linear in y(t), it be solved by eigenvalue problem (EVP) methods, by assuming an exponential solution, as usual for LODEs, of the form

y(t) = sum_k=1:n[c(k)*x(k)*exp(lambda(k)*(t-t₀))],

However, a is a real general matrix since the problem is physically real. If a is a general matrix then the EVP may have complex (real and imaginary parts) eigenvalues and eigenvectors, while there will be both left and right eigenvectors if a in not symmetric, which must be anticipated in the general case, especially since you will be asked to construct a general test matrix using random element via MATLAB's normal RNG randn(n).

The objective is to compute the solution of the LODEs by EVP methods using convenient MATLAB functions like the eigenset finder eig, the normal random number generatorrandn, the general real/complex-conjugate transpose symbol ' and complex conjugate conj.

Project Statement and Background:

Write your MATLAB code in a MATLAB m-file using the MATLAB command window menu sequence File/New/M-File and enter your MATLAB program in the MATLAB Editor/Debugger (m-file window). Save as a convenient MATLAB function in some MATLAB recognizable directory path (or add the current directory to the MATLAB path on the prompt), say as cp3me.m corresponding to function cp3me (input and output arguments are optional, while 'me' is your nick-name), then use the save and run button on top of the MATLAB m-file window or execute it in the MATLAB Command Window by function name with arguments if any. You can go back and forth between m-file editing and command line executing until output is satisfactory. Print output for each boldface heading in the code list below, appropriately labeled.

Using some pseudo-MATLAB notation:

1. Printout eps % MATLAB machine epsilon constant, eps, is min[q|q>0;float(1+q)>1].
   % (2^(-52) ~ 2.22e-16, assuming standard double precision IEEE floating point arithmetic
   % with 53 bit precision, epsilon for chopping since MATAB has turned off rounding.)

2. Printout [System,MaxArraySize]=computer, MatlabVersion=version('-release');
   % MATLAB computer system and Release Number.

3. Initialize n = 6; t0 = t₀ = 0; t_f = 4.25;
   % Project size and time horizon; "%" denotes a MATLAB style comment.

4. Let y0 = y₀ = 100*randn(n,1); a = -1.25*ones(n)+2.25*randn(n)  ;
   % Random matrices make good test examples.
   % Caution: randn(n,1) is a column n-vector, while randn(n) is a n×n matrix.

5. Let [x,dx] = eig(a)
   % x = [x(i,j)]_n×n array of right eigenvectors,
   % representing the set of all eigenvectors x(k) for k = 1:n  ;
   % dx = [lambda(i)delta(i,j)]_n×n diagonal array of right eigenvalues, satisfying a*x=x*dx  .

6. Let [z,dz] = eig(a')
   % z = [z(i,j)]_n×n array of left eigenvectors, the set z(p(k));
   % dz = [mu(i)delta(i,j)]_n×n diagonal array of corresponding left eigenvalues,
   % satisfying a'*z=z*dz or z'*a=dz'*z'  .
   % Likely, mu = conj(lambda(:,p(k))) rather than conj(lambda(:,k)) for some pivot vector p
   % as needed for biorthogonality of z and x .
   % The a' is complex transpose or a, since MATLAB recognizes the right type
   % and does the right thing.
   %
   % Caution 1: For MATLAB, the two calls to function "eig" are independent
   % according to MATLAB, so ordering of eigenvalues in arrays dz and dx
   % will likelyNOT be in the same order as expected from the theory.
   % Hence, the user must sort the complex conjugate transposed array dz',
   % so that it has the same order of eigenvalues dx within the accuracy of finite precision.
   % Also while sorting the elements of dz', the eigenvectors of z must be interchanged
   % in the same order as the eigenvalues in dz'. This can be done using a for/if loop to order dz'
   % by comparison to the order in dx, modulo floating point precision.
   %
   % Caution 2: With floating point precision,
   % NEVER ask if two real or floating point numbers A and B are equal,
   % as if the precision were infinite precision, it is NOT.
   % Always ask if |A-B| < tol, for some small tolerance tol,
   % where tol several multiples of machine epsilon,
   % depending on the application, see below.
   %
   % Caution 3: You can assume that the eigenvalues are distinct, unless you find otherwise,
   % but not their magnitudes since complex eigenvalues of real matrices
   % must come in complex conjugate pairs and complex conjugate pairs have the same magnitude.
   %
   % Caution 4: z' = ctranspose(z) is, in general, complex conjugate transpose for MATLAB,
   % i.e., if z = [real(z(j,k)) + I*imag(z(j,k))],
   % where here "I" or "J" are the MATLAB imaginary unit functions
   % (but are over-ridden if used as a variable as in 'for i=1:n');
   % then z'=[conj(z(k,j))]=[real(z(k,j)) - I*imag(z(k,j))], complex conjugate transpose.
   %
   % Caution 5: z'x should be diagonal except for floating point errors,
   % but not necessarily real, if z and x have the same ordering,
   % else it will be the permutation of a diagonal matrix.

7. Printout ordered dzr' and zr,
   % where zr is the ordered eigenvector matrix from the original z .

8. Compute 2-norm of difference (dzr'-dx) relative to eps;
   % If huge then dzr' is not ordered the same as dx, see above comments.

9. Compute 2-norm of residual (a*x-x*dx) relative to eps*norm(x);
   % If norms are large, check code your code for errors.

10. Compute 2-norm of residual (a'*zr-zr*dzr) relative to eps*norm(zr);
    % If norms are large, check code your code for errors and that you have not used dzr' here.

11. Compute diagzrtx = diag(zr'*x)
    % If zr'*x is not diagonal then the ordering of (dz,z) to (dzr,zr) is not the same as (dx,x)
    % or there is some linear algebraic error.

12. Compute c=(zr'*y0)./diagzrtx
    % See class notes for a derivation of this form
    % and note the element division operator "./" is not "/".

13. Let np = 41; dt = (tf-t0)/(np-1); tv = to:dt:tf;
    % Discretize time to keep problem algebraic, using MATLAB colon constructor.

14. Compute eigen-solution for (n×np) Y = y(tv),
    % using the for/end loop or else in full MATLAB array form (see class notes).

15. Compute 2-norm of the residual dY-a*Y relative to eps and relative (eps*norm(dY)),
    % where dY is the time rate of change and a*Y is the matrix vector product.
    % If the latter is huge (say > 1000, i.e., off by 3 significant decimal digits
    % relative to eps, but is scale dependent), then you MAY have some ERROR.

16. Professional Plot, using MATLAB plot function, the constructor solution vectors Y(i,:) as a function of time.
    % {":" means over all "j"} for i=1:n in a single professionally labeled figure,
    % combining all n components in one plot for comparison, with
    % different markings for each solution vector component.
    % Fully label plot (table, xlabel, ylabel and legend),
    % e.g., legend('y(1,:)','y(2,:)','y(3,:)','y(4,:)','y(5,:)','y(6,:)',0)

Professional Report and Output Requirements:

1. Cover Page (Title, Course, Name, Affiliation, Date, etc.).

2. Table of Contents.

3. Hardware and Compiler Specifications.

4. Executive Summary: Itemized list briefly summarizing in complete sentences the project report, describing problem, results and conclusions for the busy Boss, possibly including a significant graphical illustration of the results. Also give your recommendations to the boss for doing the problem.

5. Introduction: Motivation for the project.

6. Project: Well-documented with descriptions of the project problem.

7. Methods: Well-documented with descriptions of the variables and major steps.

8. Results: Description of output in words to explain well-labeled output tables, including a table of all "per eps" error checks, and figures. Note your MATLAB figures should be professionally labeled, i.e, the output corresponding to all bold headings 1-16.

9. Discussion: Give a good discussion of your results: with regard to prominent features and error checks.

10. Conclusions: What was primarily accomplished and what methods are recommended?

11. Acknowledgements: Acknowledge those who helped you and what nonpersonal computing resources that you used. Otherwise, it is a form of cheating and subject to servere penalties.

12. References: List the resources like books, journal papers, online information, etc., the you used. Otherwise, it is called plagerism and subject to servere penalties.

13. Appendices: Codes and other supporting material.

Project Resources:

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

2. MATLAB Help Page (Hanson).