MCS 507 --- Individual MATLAB Project 3 --- Fall 2004
Eigen-Solutions of Linear ODEs for Industrial Rate Computations:
MATLAB Computation: m-files, eig, randn and complex conj
MATLAB Output is
Due Monday 08 November 2004 in Class
.
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(t0) = y0,
t0 < t <= tf ,
where y(t) = [yi(t)]n×1
is a time-dependent vector,
a = [ai,j]n×n
is a given constant array, and
y0 = [y0,i]n×1
is a specified initial vector.
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) = sumk=1:n[c(k)*x(k)*exp(lambda(k)*(t-t0))],
after linear superposition,
where lambda(k) are set n scalar eigenvalues, likely distinct,
x(k) = [x(i,k)]n×1
are the corresponding eigenvectors,
and c = [c(i)]n×1 are multiplicative constants of integration.
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:
- 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.)
- Printout [System,MaxArraySize]=computer,
MatlabVersion=version('-release');
% MATLAB computer system and Release Number.
- Initialize n = 6; t0 = t0 = 0;
tf = 4.25;
% Project size and time horizon; "%" denotes a MATLAB style comment.
- Let y0 = y0 = 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.
-
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 .
-
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.
- Printout ordered dzr' and zr,
% where zr is the ordered eigenvector matrix from the original z .
- Compute 2-norm of difference (dzr'-dx) relative to eps;
% If huge then dzr' is not ordered the same as dx, see above comments.
- Compute 2-norm of residual (a*x-x*dx) relative to eps*norm(x);
% If norms are large, check code your code for errors.
- 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.
- 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.
- Compute c=(zr'*y0)./diagzrtx
% See class notes for a derivation of this form
% and note the element division operator "./" is not "/".
- Let np = 41; dt = (tf-t0)/(np-1); tv = to:dt:tf;
% Discretize time to keep problem algebraic,
using MATLAB colon constructor.
- Compute eigen-solution for (n×np) Y = y(tv),
% using the for/end loop or else in full MATLAB array form (see class notes).
- 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.
- 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:
- Cover Page (Title, Course, Name, Affiliation, Date, etc.).
- Table of Contents.
- Hardware and Compiler Specifications.
- 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.
- Introduction: Motivation for the project.
- Project:
Well-documented with descriptions of the project problem.
- Methods:
Well-documented with descriptions of the variables and major steps.
- 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.
- Discussion: Give a good discussion of your results:
with regard to prominent features and error checks.
- Conclusions: What was primarily accomplished and what methods are
recommended?
- 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.
- References: List the resources like books, journal papers, online
information, etc., the you used. Otherwise, it is called plagerism and subject
to servere penalties.
- Appendices: Codes and other supporting material.
Project Resources:
- UIC PCLabs MATLAB available software. Also check out your
departmental computers such as in EECS and MSCS.
-
MATLAB Help Page (Hanson).
Web Source: http://www.math.uic.edu/~hanson/mcs507/cp3f03.html
MCS 507 HomePage: http://www.math.uic.edu/~hanson/mcs507/
Email Comments or Questions to Professor Hanson,
hanson A T math uic edu .