Path: news.mathworks.com!not-for-mail
From: "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
Newsgroups: comp.soft-sys.matlab
Subject: Re: calculating a definite integral
Date: Sat, 19 Apr 2008 00:28:02 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 95
Message-ID: <fubeai$qaa$1@fred.mathworks.com>
References: <fu2f72$s9b$1@fred.mathworks.com> <fuas3d$du7$1@fred.mathworks.com>
Reply-To: "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1208564882 26954 172.30.248.35 (19 Apr 2008 00:28:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 19 Apr 2008 00:28:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:463999


"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in 
message <fuas3d$du7$1@fred.mathworks.com>...
>   Aviral, if you substitute z = exp(j*w) in your definite integrals, they 
become 
> contour integrals taken with respect to z and moving counterclockwise 
> around the unit circle in the complex plane for the matrix-valued 
integrands
> 
>  -j*inv(z*I-A)*B*C*inv(1/z*I-A)/z.
> 
> Cramer's rule shows us that both these two inverses are rational functions 
of 
> z for each of the elements of the matrices.  Thus each of the elements of 
the 
> above integrands is some rational and therefore meromorphic function of 
z.  
> We can find their integrals around the closed circle by simply evaluating all 
> their single poles which lie within the unit circle, by virtue of the theory of 
> residues.
> 
>   In principle, matlab's symbolic toolbox ought to be able to determine 
these 
> residues and therefore evaluate the integrals without resorting to numerical 
> quadrature.  I am not sure how well the toolbox can actually perform at this 
> for large size matrices because the complexity of expressions increases 
> rapidly as their size increases, but you can try putting it through its paces 
> along these lines and see how well it performs.
> 
> Roger Stafford
------------
  It has occurred to me on thinking of the matter at greater length, that the 
poles of the meromorphic function I referred to earlier are precisely the 
eigenvalues of the matrix A and their reciprocals, a fact which should 
considerably ease the burden of symbolic toolbox manipulation.

  Again referring to Cramer's rule, the inverse of z*I-A can be seen to be the 
ratio of two polynomials in z.  The denominator is just det(z*I-A), which in 
turn equals (z-L1)*(z-L2)*...*(z-Ln) where L1, L2, ..., Ln are the eigenvalues of 
A.  The numerator will be some polynomial of no more than the n-1 degree.

  A similar situation exists for the inverse of 1/z*I-A in respect to the 
reciprocals of A's eigenvalues.  If the numerator and denominator from 
Cramer's rule are multiplied by (-z)^n, and divided by the product of the 
eigenvalues, then the denominator becomes (z-1/L1)*(z-1/L2)*...*(z-1/Ln) 
and the numerator will become z times a polynomial of degree no more than 
n-1.  This first z factor will cancel out the z in the denominator of the original 
integrand, namely at the right end of 

 -j*inv(z*I-A)*B*C*inv(1/z*I-A)/z.

  Provided no eigenvalue lies right on the unit circle itself, (which would create 
a singularity smack on the path of integration,) for each eigenvalue, either it 
lies in the interior of the unit circle or its reciprocal does, but not both.  This 
gives us precisely n poles to be found in the interior of our meromorphic 
integrand.  The definite integral you seek is then equal to 2*pi*j times the 
sum of the residues at these n poles.

  One way of determining these residues is this.  Suppose eigenvalue L lies 
inside the unit circle.  If we take the symbolic derivative of the above 
denominator, det(z*I-A), with respect to z, and then assign z in it to the 
numerical value z = L, that is the equivalent of having divided det(z*I-A) by 
z-L symbolically and then set z numerically equal to L in the remaining 
factor.  However, doing it this way avoids the problem of not possessing a 
symbolic expression for L for sizes of A greater than 4x4.  The entire 
integrand, but corrected by having the z-L factor in the denominator 
removed in this fashion, can then be evaluated numerically at z = L and the 
result would be just the residue at this pole.  An analogous method would be 
used to adjust det(1/z*I-A) in a similar way for an eigenvalue occurring 
outside the circle.

  This is admittedly something I haven't tried out on my computer as yet, but 
the theory appears valid to me.  The main practical difficulty I foresee is that 
the toolbox may have trouble obtaining the above symbolic inverses for large 
size matrices, A.  Some of the expressions may not bear looking at by us 
mere mortals for fear we will turn to stone, though we probably don't really 
need to look at them to obtain numerical results from them.  (On my system 
no expression string is permitted to exceed 8192 characters, so I am 
protected to a degree from excessive shock from such mind-boggling 
visions.  On the other hand, I can't get the symbolic inverse of z*I-A for 
anything past a 3 x 3 matrix, A.)

  As to whether this method is superior to doing the integration by numerical 
quadrature I cannot be sure.  Numerical quadrature could encounter severe 
accuracy problems if some of the eigenvalues happen to lie very near the unit 
circle, as might be anticipated.  Using the residue method would avoid such 
problems for the most part, except that the residue obtained for the pole just 
inside the circle, even though accurate, would presumably be very large.  On 
the other hand, it is probably much less laborious to prepare the code for 
doing this problem numerically.  In any case, you could regard this little 
dissertation as strengthening your theoretical understanding of the nature of 
the task facing you.

Roger Stafford