Home > Archive > Fortran > May 2004 > Searching zeros of complex function
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
Searching zeros of complex function
|
|
| Miguel Ángel Solano Vérez 2004-05-14, 5:31 am |
| I have a square complex matrix. A cuadrant of this matrix has terms of the
type exp(j*beta*d). being beta the unknown. I'm looking for a rutine
(fortran) which gives the zeros of the determinant of this matrix, or simple
a rutine which gives zeros of a complex function. I have been looking in the
IMSL package but i do not find anything.
Thanks,
Miguel A. Solano
****************
| |
| Gerry Thomas 2004-05-14, 9:31 am |
|
"Miguel Ángel Solano Vérez" <solanom@unican.es> wrote in message
news:c81se1$bbe$1@ccserver13.unican.es...
> I have a square complex matrix. A cuadrant of this matrix has terms of
the
> type exp(j*beta*d). being beta the unknown. I'm looking for a rutine
> (fortran) which gives the zeros of the determinant of this matrix, or
simple
> a rutine which gives zeros of a complex function. I have been looking in
the
> IMSL package but i do not find anything.
>
Miguel,
As det(A) is a constant it doesn't have any roots. det(A-zI)=0 is a
polynomial of the order of A, n say, and its n roots are the eigenvalues of
A. IMSL has routines to determine the eigenvalues of general complex A.
--
E&OE
HTH,
Gerry T.
| |
| Rich Townsend 2004-05-14, 10:32 am |
| Gerry Thomas wrote:
> "Miguel Ángel Solano Vérez" <solanom@unican.es> wrote in message
> news:c81se1$bbe$1@ccserver13.unican.es...
>
>
> the
>
>
> simple
>
>
> the
>
>
>
> Miguel,
>
> As det(A) is a constant it doesn't have any roots. det(A-zI)=0 is a
> polynomial of the order of A, n say, and its n roots are the eigenvalues of
> A. IMSL has routines to determine the eigenvalues of general complex A.
>
The OP specifies that beta is an unknown, therefore det(A) is not a
constant, but a function of beta.
For general-purpose root finding, I use Muller's method. Although it is
usually presented as an algorithm for finding roots of polynomials, it
works perfectly well for finding the roots of generic complex functions.
A good implementation of the method is given by Traub (1964,
"Iterative Methods for the Solution of Equations", Prentice-Hall,
Englewood Cliffs).
Other options are a simple secant method (see, e.g., Numerical Recipes
in Fortran), or -- if you can somehow calculate the derivative of the
determinant wrt beta -- a Newton-Raphson algorithm.
cheers,
Rich
| |
| Gerry Thomas 2004-05-14, 10:32 am |
|
"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:c82h2k$dkd$1@scrotar.nss.udel.edu...
> Gerry Thomas wrote:
in[color=darkred]
eigenvalues of[color=darkred]
>
> The OP specifies that beta is an unknown, therefore det(A) is not a
> constant, but a function of beta.
det(A) is zero iff A is singular. As this isn't necessarily true, det(A) is
a constant.
--
E&OE
HTH,
Gerry T.
| |
| Rich Townsend 2004-05-14, 10:32 am |
| Gerry Thomas wrote:
> "Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
> news:c82h2k$dkd$1@scrotar.nss.udel.edu...
>
>
> in
>
>
> eigenvalues of
>
>
>
> det(A) is zero iff A is singular. As this isn't necessarily true, det(A) is
> a constant.
>
>
But the elements of A depend on beta. Therefore, det(A) depends on beta,
and may be zero (indicating a singular matrix) for certain values of
beta. What part of "depends on beta" are you unable to comprehend?
| |
| Gerry Thomas 2004-05-14, 11:32 am |
|
"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:c82ipj$dvi$1@scrotar.nss.udel.edu...
> Gerry Thomas wrote:
of[color=darkred]
A.[color=darkred]
det(A) is[color=darkred]
>
> But the elements of A depend on beta. Therefore, det(A) depends on beta,
> and may be zero (indicating a singular matrix) for certain values of
> beta. What part of "depends on beta" are you unable to comprehend?
You're making an unsubstantiated assumption about the structure of A.
Here's a hint:
abs(det(A)) <= prod(1<=k<=n)norm-2(A(:,k) )
There are circumstances in which det(A)=0 but the OP hasn't stated this to
be the case and neither should you.
--
Ciao,
Gerry T.
______
"...as I further could have told you as brisk as your D.B.C.
behaviouristically paillet, with a coat of homoid icing which is in reality
only a done by chance ridiculisation of the whoo-whoo and where's hairs
theorics of Winestain." -- James Joyce, in Finnegans Wake.
| |
| Rich Townsend 2004-05-14, 11:32 am |
| Gerry Thomas wrote:
> "Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
> news:c82ipj$dvi$1@scrotar.nss.udel.edu...
>
>
> of
>
>
> A.
>
>
> det(A) is
>
>
>
> You're making an unsubstantiated assumption about the structure of A.
> Here's a hint:
>
> abs(det(A)) <= prod(1<=k<=n)norm-2(A(:,k) )
>
> There are circumstances in which det(A)=0 but the OP hasn't stated this to
> be the case and neither should you.
>
I'm making no assumptions whatsoever. But you have claimed (twice) that
det(A) is a constant. How can you reconcile this with the fact that some
of its elements depend on the variable beta -- as the OP has stated
quite clearly?
| |
| Gerry Thomas 2004-05-14, 12:32 pm |
|
"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:c82mig$f2i$1@scrotar.nss.udel.edu...
> Gerry Thomas wrote:
rutine[color=darkred]
or[color=darkred]
looking[color=darkred]
beta,[color=darkred]
to[color=darkred]
>
> I'm making no assumptions whatsoever. But you have claimed (twice) that
> det(A) is a constant. How can you reconcile this with the fact that some
> of its elements depend on the variable beta -- as the OP has stated
> quite clearly?
Your logic is lacking. I fully appreciate that a quadrant of A depends on b
but that doesn't imply that det(A) is zero, just that it's a constant and
is a function of b.
I'm sure the OP will figure it out.
--
E&OE
Ciao,
Gerry T.
| |
| Richard Maine 2004-05-14, 12:32 pm |
| "Gerry Thomas" <gfthomas@sympatico.ca> writes:
> You're making an unsubstantiated assumption about the structure of A.
No assumptions seem needed. The OP seemed pretty explicit about it
to me.
> There are circumstances in which det(A)=0 but the OP hasn't stated this to
> be the case and neither should you.
The circumstances are these things called "roots". You know - what
the OP said he was looking for.
Now perhaps the OP, Rich Townsend, and/or now myself might have misused
some term in some trivial way that you are trying to be obscure about
pointing out. Could be; I don't much care if you want to belabor
whatever point that might be.
In the meantime, the OP has sucessfully communicated his question
to me, and apparently also to Rich Townsend. If you were unable to
understand the question, don't worry about it, as Rich Townsend
already offered a constructive reply.
--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain
| |
| Rich Townsend 2004-05-14, 12:32 pm |
| Gerry Thomas wrote:
<snip>
> Your logic is lacking. I fully appreciate that a quadrant of A depends on b
> but that doesn't imply that det(A) is zero, just that it's a constant and
> is a function of b.
LOL! You remind me of that great quote "Constants aren't, variables
won't". I think your last sentence gets the prize for the most
contradictory statement of the year.
| |
| Richard Maine 2004-05-14, 12:32 pm |
| "Gerry Thomas" <gfthomas@sympatico.ca> writes:
> ...it's a constant and is a function of b.
I'm clearly doing the wrong drugs. Caffeine isn't going to get
me there. :-(
--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain
| |
| Gerry Thomas 2004-05-14, 1:31 pm |
|
"Richard Maine" <nospam@see.signature> wrote in message
news:m1pt97ni18.fsf@macfortran.local...
> "Gerry Thomas" <gfthomas@sympatico.ca> writes:
>
>
> I'm clearly doing the wrong drugs. Caffeine isn't going to get
> me there. :-(
>
The entries of a Hilbert matrix are functions of its row and column indices
and the entries are constants.
--
E&OE
HTH,
Gerry T.
______
"Just say no to drugs." -- Nancy Reagan.
| |
| Rich Townsend 2004-05-14, 1:31 pm |
| Gerry Thomas wrote:
> "Richard Maine" <nospam@see.signature> wrote in message
> news:m1pt97ni18.fsf@macfortran.local...
>
>
>
> The entries of a Hilbert matrix are functions of its row and column indices
> and the entries are constants.
>
And where did the OP mention that they have a Hilbert matrix? LOL, you
just keep on digging, don't you!
| |
| Gerry Thomas 2004-05-14, 1:31 pm |
|
"Miguel Ángel Solano Vérez" <solanom@unican.es> wrote in message
news:c81se1$bbe$1@ccserver13.unican.es...
> I have a square complex matrix. A cuadrant of this matrix has terms of
the
> type exp(j*beta*d). being beta the unknown. I'm looking for a rutine
> (fortran) which gives the zeros of the determinant of this matrix, or
simple
> a rutine which gives zeros of a complex function. I have been looking in
the
> IMSL package but i do not find anything.
>
> Thanks,
>
> Miguel A. Solano
> ****************
>
>
Since the determinant of the quadrant its simply the product of its
eigenvalues the roots of the determinant are the eigenvalues of the
quardrant, a freshman exercise, :-).
--
E&OE
HTH,
Gerry T.
| |
| Gerry Thomas 2004-05-14, 1:31 pm |
|
"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:c82r28$gc9$1@scrotar.nss.udel.edu...
> Gerry Thomas wrote:
indices[color=darkred]
>
> And where did the OP mention that they have a Hilbert matrix?
Strawmaning again, eh?
>LOL, you just keep on digging, don't you!
I believe the feeble minded like to laugh at their own lack of wit. Thanks
for the demo.
--
E&OE
Ciao,
Gerry T.
______
"Ah, Klinger, my constant reminder that Darwin was right!." -- Maj Charles
Winchester III, the 4077th M*A*S*H
| |
| Steven G. Kargl 2004-05-14, 2:31 pm |
| In article <f73pc.34245$dr1.940468@news20.bellglobal.com>,
"Gerry Thomas" <gfthomas@sympatico.ca> writes:
>
> "Miguel Ángel Solano Vérez" <solanom@unican.es> wrote in message
> news:c81se1$bbe$1@ccserver13.unican.es...
>
> As det(A) is a constant it doesn't have any roots. det(A-zI)=0 is a
> polynomial of the order of A, n say, and its n roots are the eigenvalues of
> A. IMSL has routines to determine the eigenvalues of general complex A.
>
| 1 exp(j*beta*d) |
| | = det(A) = -exp(j*beta*d)
| 1 0 |
How is det(A) a constant when the OP stated beta is an unknown?
--
Steve
http://troutmask.apl.washington.edu/~kargl/
| |
| Steven G. Kargl 2004-05-14, 2:31 pm |
| In article <c81se1$bbe$1@ccserver13.unican.es>,
"Miguel Ángel Solano Vérez" <solanom@unican.es> writes:
> I have a square complex matrix. A cuadrant of this matrix has terms of the
> type exp(j*beta*d). being beta the unknown. I'm looking for a rutine
> (fortran) which gives the zeros of the determinant of this matrix, or simple
> a rutine which gives zeros of a complex function. I have been looking in the
> IMSL package but i do not find anything.
>
Miguel,
Do you have some idea where in the complex plane that values
of beta may cause the determinant to be zero? If so, you can
compute the winding number for a simply closed contour to
determine the number of zeros and poles enclosed by the contour.
If the winding number is 1 and you have reasons to believe that
no poles are within the contour, then computing a simple cauchy
integral will give you the value of beta. You can find a detailed
description of the method in Appendix B of SG Kargl and PL Marston,
J. Acoust. Soc. Am., 85, 1014-1028 (1989).
--
Steve
http://troutmask.apl.washington.edu/~kargl/
| |
| glen herrmannsfeldt 2004-05-14, 3:31 pm |
| Gerry Thomas wrote:
> "Miguel Ángel Solano Vérez" <solanom@unican.es> wrote in message
> news:c81se1$bbe$1@ccserver13.unican.es...
(snip)
[color=darkred]
> Since the determinant of the quadrant its simply the product of its
> eigenvalues the roots of the determinant are the eigenvalues of the
> quardrant, a freshman exercise, :-).
Some of us haven't been freshman for a fairly large number
of years, but that one does sound familiar.
The ones I am remembering were in the form
exp(t*A) where t is a variable and A is a matrix,
and * is the Fortran multiplication operator.
I don't have my Freshman calculus book anymore, though.
-- glen
| |
| Gerry Thomas 2004-05-14, 3:31 pm |
|
"Steven G. Kargl" <kargl@troutmask.apl.washington.edu> wrote in message
news:c82tii$ahq$1@nntp6.u.washington.edu...
> In article <f73pc.34245$dr1.940468@news20.bellglobal.com>,
> "Gerry Thomas" <gfthomas@sympatico.ca> writes:
rutine[color=darkred]
eigenvalues of[color=darkred]
>
> | 1 exp(j*beta*d) |
> | | = det(A) = -exp(j*beta*d)
> | 1 0 |
>
> How is det(A) a constant when the OP stated beta is an unknown?
The problem posed was how to find the roots of det(A) for arbitrary b and
was not how to find b. In the current post to the OP I gave him a hint and
in a followup post I gave him the answer to his homework assignment. As for
your question, specify b and you'll know det(A), a constant. Keep in mind
that this burlesque is not MATLAB/Maple or Mathematica. It's interesting to
note that on sci.math.num-analysis his post was pretty well ignored as it's
customary for respondents there to spot homework and at most to provide a
hint. Here it degenerates to vain pillory by delicate and fragile
narcissistic egos who provide a veritable anthropological feast of their
untempered ids.
The show's over, or is it Mission Accomplished, :-).
--
E&OE
Ciao,
Gerry T.
| |
| Gerry Thomas 2004-05-14, 3:31 pm |
|
"glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:YY7pc.1076$gr.63328@attbi_s52...
> Gerry Thomas wrote:
>
>
>
> (snip)
>
>
> Some of us haven't been freshman for a fairly large number
> of years, but that one does sound familiar.
>
It's about the most popular thing that freshmen remember about det's.
> The ones I am remembering were in the form
>
> exp(t*A) where t is a variable and A is a matrix,
> and * is the Fortran multiplication operator.
>
But don't use the Fortran exp to do this if you're interested in a
meaningful result (barring that you not just doing an elementwise
exponentiation).
> I don't have my Freshman calculus book anymore, though.
>
Nor do I.
--
E&OE
Ciao,
Gerry T.
| |
| Rich Townsend 2004-05-14, 3:31 pm |
| Gerry Thomas wrote:
> The problem posed was how to find the roots of det(A) for arbitrary b and
> was not how to find b.
And here is were the disagreement lies. From the title of this thread,
plus this phrase from the original post:
"A cuadrant of this matrix has terms of the
type exp(j*beta*d). being beta the unknown."
....it seems clear to me that the OP s s to find beta, beta being the
unknown (variable) quantity which allows the equation det(A)=0 to be
satisfied. Your reading of the problem appears to be different, and I
guess we'll have to wait for the OP to get back with a clarification.
In the meantime, please feel free to carry on entertaining us all.
cheers,
Rich
| |
| Gerry Thomas 2004-05-14, 3:31 pm |
|
"Gerry Thomas" <gfthomas@sympatico.ca> wrote in message
news:l68pc.80432$FH5.1838553@news20.bellglobal.com...
>
> "glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
> news:YY7pc.1076$gr.63328@attbi_s52...
>
> It's about the most popular thing that freshmen remember about det's.
>
>
> But don't use the Fortran exp to do this if you're interested in a
> meaningful result (barring that you not just doing an elementwise
> exponentiation).
PS
I just remembered that det(exp(A))=exp(Trace(A)) which is neat.
[...]
| |
| Gerry Thomas 2004-05-14, 3:31 pm |
|
"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:c833l2$iqk$1@scrotar.nss.udel.edu...
[...]
> And here is were the disagreement lies. From the title of this thread,
> plus this phrase from the original post:
>
> "A cuadrant of this matrix has terms of the
> type exp(j*beta*d). being beta the unknown."
>
> ...it seems clear to me that the OP ss to find beta, beta being the
> unknown (variable) quantity which allows the equation det(A)=0 to be
> satisfied. Your reading of the problem appears to be different, and I
> guess we'll have to wait for the OP to get back with a clarification.
>
That's one take.
> In the meantime, please feel free to carry on entertaining us all.
>
Careful, your id is out.
--
E&OE
Ciao,
Gerry T.
______
"The best lesson life has taught me is that the idiots in many cases are
right." -- Sir Winston S. Churchill.
| |
| Richard Maine 2004-05-14, 4:32 pm |
| "Gerry Thomas" <gfthomas@sympatico.ca> writes:
> The problem posed was how to find the roots of det(A) for arbitrary b and
> was not how to find b.
Ah. Finally a statement on this subject that helps me to actually
understand what it is that we disagree about. (Yes, I'm sure my former
failings say no end of uncomplimentary things about me, as can
no doubt be detailed by those interested in detailing such things.)
I still disagree, but now I understand at least one basis of the
disagreement.
How to find value of b that makes the determinant zero was *EXACTLY*
what I understood the OP to be asking. I still think that the most
likely interpretation of the request, though I can see at least the
possibility of reading it the other way.
Just as I try to carefully distinguish my personal Fortran style
preferences from the dictates of the standard, likewise I prefer
to make it clear when I am giving my personal interpretation of
what someone probably meant. Thus my paragraph above avoids
wording that implies absolute truth; I don't have the necessary
divine connections.
I'll ignore the rest of your post as being at a ...um... different
level of professionalism than I prefer to exhibit. If you are so
inclined, feel free to interpret "different" as "higher".
--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain
| |
| Miguel Ángel Solano Vérez 2004-05-14, 5:31 pm |
| Thanks very much for the answers and the "discussion". I think that my
problem is not if det(A)=constant or not. Obvioulsy, in a numerical
computing language (fortran, for example) det(A) is constant for a value of
beta. But I have to find the value of beta that makes det(A)=0. This is very
common in electromagnetics: finding zeros of "complex characeristic
equations". So, I can compute det(A) without problems (this would be the
characteristic function), but when the matrix is a bit large (20x20 or more)
I know that in some cases I have to find a real value of beta (which makes
det(A)=0) and I do not. The strategy for finding the zeros, that I are
using, is very poor, because is not based in any mathematical basis: simply
I see where the sing of det(A) changes (both real and imaginary parts) and I
go narrowing the interval. Also, I have a fortran routine based in Muller's
Method for complex functions that I used some time ago for a complex
trigonometric function, and it is not work as well as I would like: some
solutions were not found.
In summary, as you see I'm not an expert in Mathematics and perhaps this is
the problem. But I am very gratefully for yor help. And, please, don't get
 .
I only was looking for a "good" fortran routine which solves me the problem.
Thanks very much to everybody.
Cheers,
Miguel Angel
*************
"Steven G. Kargl" <kargl@troutmask.apl.washington.edu> escribió en el
mensaje news:c82tc3$aho$1@nntp6.u.washington.edu...
> In article <c81se1$bbe$1@ccserver13.unican.es>,
> "Miguel Ángel Solano Vérez" <solanom@unican.es> writes:
the[color=darkred]
simple[color=darkred]
the[color=darkred]
>
> Miguel,
>
> Do you have some idea where in the complex plane that values
> of beta may cause the determinant to be zero? If so, you can
> compute the winding number for a simply closed contour to
> determine the number of zeros and poles enclosed by the contour.
> If the winding number is 1 and you have reasons to believe that
> no poles are within the contour, then computing a simple cauchy
> integral will give you the value of beta. You can find a detailed
> description of the method in Appendix B of SG Kargl and PL Marston,
> J. Acoust. Soc. Am., 85, 1014-1028 (1989).
>
> --
> Steve
> http://troutmask.apl.washington.edu/~kargl/
| |
| Gerry Thomas 2004-05-14, 5:32 pm |
|
"Miguel Ángel Solano Vérez" <solanom@unican.es> wrote in message
news:c837pp$jbv$1@ccserver13.unican.es...
> Thanks very much for the answers and the "discussion". I think that my
> problem is not if det(A)=constant or not. Obvioulsy, in a numerical
> computing language (fortran, for example) det(A) is constant for a value
of
> beta. But I have to find the value of beta that makes det(A)=0. This is
very
> common in electromagnetics: finding zeros of "complex characeristic
> equations". So, I can compute det(A) without problems (this would be the
> characteristic function), but when the matrix is a bit large (20x20 or
more)
> I know that in some cases I have to find a real value of beta (which
makes
> det(A)=0) and I do not. The strategy for finding the zeros, that I are
> using, is very poor, because is not based in any mathematical basis:
simply
> I see where the sing of det(A) changes (both real and imaginary parts)
and I
> go narrowing the interval. Also, I have a fortran routine based in
Muller's
> Method for complex functions that I used some time ago for a complex
> trigonometric function, and it is not work as well as I would like: some
> solutions were not found.
>
> In summary, as you see I'm not an expert in Mathematics and perhaps this
is
> the problem. But I am very gratefully for yor help. And, please, don't
get
> .
> I only was looking for a "good" fortran routine which solves me the
problem.
>
> Thanks very much to everybody.
>
> Cheers,
>
> Miguel Angel
> *************
>
Thank you for the clarification. If you have access to MATLAB I suggest
looking at
http://web.comlab.ox.ac.uk/projects...pectra/eigtool/
which I have found to be an insightful tool for the kind of problem you are
dealing with.
--
E&OE
Good Luck,
Gerry T.
| |
| glen herrmannsfeldt 2004-05-14, 5:32 pm |
| Gerry Thomas wrote:
(I wrote)
[color=darkred]
[color=darkred]
> But don't use the Fortran exp to do this if you're interested in a
> meaningful result (barring that you not just doing an elementwise
> exponentiation).
There is always the problem of how to write mathematical
operations in ASCII (or EBCDIC) characters. I do like the
exp() notation better than fake superscripts.
-- glen
| |
| Rich Townsend 2004-05-14, 5:32 pm |
| Miguel Ángel Solano Vérez wrote:
> Thanks very much for the answers and the "discussion". I think that my
> problem is not if det(A)=constant or not. Obvioulsy, in a numerical
> computing language (fortran, for example) det(A) is constant for a value of
> beta. But I have to find the value of beta that makes det(A)=0.
Which now resolves the dispute!
> This is very
> common in electromagnetics: finding zeros of "complex characeristic
> equations". So, I can compute det(A) without problems (this would be the
> characteristic function), but when the matrix is a bit large (20x20 or more)
> I know that in some cases I have to find a real value of beta (which makes
> det(A)=0) and I do not. The strategy for finding the zeros, that I are
> using, is very poor, because is not based in any mathematical basis: simply
> I see where the sing of det(A) changes (both real and imaginary parts) and I
> go narrowing the interval. Also, I have a fortran routine based in Muller's
> Method for complex functions that I used some time ago for a complex
> trigonometric function, and it is not work as well as I would like: some
> solutions were not found.
This is a general problem with root-finding algorithm's -- they don't
always find all the roots one wants. It sounds like the approach you are
using at the moment is a form of bisection, generalized to complex
variables. A related approach to this bisection is Lehmer's method,
which you can read about in "Numerical Methods that Work" (Forman S.
Acton, 1970, published by Harper & Row), and you may be able to learn from.
Although he misunderstood your question, Gerry Thomas' responses contain
a few flowers amongst the thorny bullshit and bluster. If you can
express A like this:
A = A' - f(beta) I,
where A' is a matrix which is independent of beta, f(beta) is a scalar
function of beta, and I is the identity matrix, then solving det(A) = 0
is equivalent to finding the eigenvalues of A'. Once we have found the
eigenvalues, then finding beta comes down to solving the scalar equation
lambda_{i} = f(beta)
for each eigenvalue lambda_{i} (there will be as many of these as the
dimension of A' -- but they won't all necessarily be distinct from one
another). So your problem is reduced to two steps: find the eigenvalues
of A' (which can be done by one of the IMSL routines), and then solve
the above implicit equation for beta.
>
> In summary, as you see I'm not an expert in Mathematics and perhaps this is
> the problem. But I am very gratefully for yor help. And, please, don't get
> .
> I only was looking for a "good" fortran routine which solves me the problem.
Nobody was getting with you -- merely a disagreement between the
various personalities on this newsgroup :)
cheers,
Rich
| |
| Rich Townsend 2004-05-14, 5:32 pm |
| Rich Townsend wrote:
<snip>
> This is a general problem with root-finding algorithm's -- they don't
Boy, I certainly need more coffee -- I've committed the cardinal sin of
using an apostrophe to announce an incoming 's'. D'oh!
Rich
| |
| Rich Townsend 2004-05-14, 6:31 pm |
| Gerry Thomas wrote:
<snip>
> It's interesting to
> note that on sci.math.num-analysis his post was pretty well ignored as it's
> customary for respondents there to spot homework and at most to provide a
> hint.
For those following Gerry's transvestite drag act, pop across to
sci.math.num-analysis -- he's making just as much of a fool of himself
over there as he is in c.l.f!
cheers,
Rich
| |
| Steven G. Kargl 2004-05-14, 6:31 pm |
| In article <c837pp$jbv$1@ccserver13.unican.es>,
"Miguel Ángel Solano Vérez" <solanom@unican.es> writes:
> Thanks very much for the answers and the "discussion". I think that my
> problem is not if det(A)=constant or not. Obvioulsy, in a numerical
> computing language (fortran, for example) det(A) is constant for a value of
> beta. But I have to find the value of beta that makes det(A)=0. This is very
> common in electromagnetics: finding zeros of "complex characeristic
> equations". So, I can compute det(A) without problems (this would be the
> characteristic function), but when the matrix is a bit large (20x20 or more)
> I know that in some cases I have to find a real value of beta (which makes
> det(A)=0) and I do not. The strategy for finding the zeros, that I are
> using, is very poor, because is not based in any mathematical basis: simply
> I see where the sing of det(A) changes (both real and imaginary parts) and I
> go narrowing the interval. Also, I have a fortran routine based in Muller's
> Method for complex functions that I used some time ago for a complex
> trigonometric function, and it is not work as well as I would like: some
> solutions were not found.
I understand. The appendix in my paper describes exactly
what you want to do. You need to be able to compute the
contour integral. Let f(b) = det(A). You want to find a
value B such that f(B) = 0. The winding number theorem
and a Cauchy integral will allow you to do this. Briefly,
the winding number theorem gives
1 / f'(b)
N - P = ------ | ------- db
i*2*pi / f(b)
The integral is a contour integral around a simply closed contour.
N is the number of zeros within the contour and P is the number of
poles. You need to be able to compute f(b) and f'(b) where f'(b)
is the derivative of f(b) with respect to b.
If there are no poles, then P = 0. If N = 1, you have a single
isolated zero. You can refine the contour to guarantee N=1, P=0.
Once you have N=1, P=0, then
1 / b f'(b)
B = ------ | -------- db
i*2*pi / f(b)
The trick is to choose a circular contour around a candidate location
of for B. In my application B was not very sensitive to the radius
of contour, and one can evaluate a good estimate of B with a fairly
small number of evaluations of f(b).
--
Steve
http://troutmask.apl.washington.edu/~kargl/
| |
| Gerry Thomas 2004-05-14, 6:31 pm |
|
"Steven G. Kargl" <kargl@troutmask.apl.washington.edu> wrote in message
news:c83c4b$jb8$1@nntp6.u.washington.edu...
> In article <c837pp$jbv$1@ccserver13.unican.es>,
> "Miguel Ángel Solano Vérez" <solanom@unican.es> writes:
value of[color=darkred]
very[color=darkred]
the[color=darkred]
more)[color=darkred]
makes[color=darkred]
simply[color=darkred]
and I[color=darkred]
Muller's[color=darkred]
some[color=darkred]
>
> I understand. The appendix in my paper describes exactly
> what you want to do. You need to be able to compute the
> contour integral. Let f(b) = det(A). You want to find a
> value B such that f(B) = 0. The winding number theorem
> and a Cauchy integral will allow you to do this. Briefly,
> the winding number theorem gives
>
> 1 / f'(b)
> N - P = ------ | ------- db
> i*2*pi / f(b)
>
> The integral is a contour integral around a simply closed contour.
> N is the number of zeros within the contour and P is the number of
> poles. You need to be able to compute f(b) and f'(b) where f'(b)
> is the derivative of f(b) with respect to b.
>
> If there are no poles, then P = 0. If N = 1, you have a single
> isolated zero. You can refine the contour to guarantee N=1, P=0.
> Once you have N=1, P=0, then
>
> 1 / b f'(b)
> B = ------ | -------- db
> i*2*pi / f(b)
>
> The trick is to choose a circular contour around a candidate location
> of for B. In my application B was not very sensitive to the radius
> of contour, and one can evaluate a good estimate of B with a fairly
> small number of evaluations of f(b).
>
This sounds interesting. What's the reference on your www? Is there any
connection with pseudospectra techniques?
Thanks,
Gerry T.
> --
> Steve
> http://troutmask.apl.washington.edu/~kargl/
| |
| Steven G. Kargl 2004-05-14, 7:31 pm |
| In article <3Zapc.81451$FH5.1890949@news20.bellglobal.com>,
"Gerry Thomas" <gfthomas@sympatico.ca> writes:
>
> "Steven G. Kargl" <kargl@troutmask.apl.washington.edu> wrote in message
> news:c83c4b$jb8$1@nntp6.u.washington.edu...
>
> This sounds interesting. What's the reference on your www? Is there any
> connection with pseudospectra techniques?
>
Here's the complete citation
Steven G. Kargl and Philip L. Marston, ``Observations and modeling of
the backscattering of short tone bursts from a spherical shell: Lamb
wave echoes, glory, and axial reverberations,'' Journal Acoustical
Society America, vol 85, pp 1014-1028 (1989).
I don't have a PDF version of this paper. The code to implement the
above two integrals is on the order of 40 lines.
I don't know of any connection with pseudospectra techniques. It
is a straight forward application of contour integrals in the
complex plane.
--
Steve
http://troutmask.apl.washington.edu/~kargl/
| |
| Rich Townsend 2004-05-14, 7:31 pm |
| Steven G. Kargl wrote:
> In article <3Zapc.81451$FH5.1890949@news20.bellglobal.com>,
> "Gerry Thomas" <gfthomas@sympatico.ca> writes:
>
>
>
> Here's the complete citation
>
> Steven G. Kargl and Philip L. Marston, ``Observations and modeling of
> the backscattering of short tone bursts from a spherical shell: Lamb
> wave echoes, glory, and axial reverberations,'' Journal Acoustical
> Society America, vol 85, pp 1014-1028 (1989).
>
> I don't have a PDF version of this paper. The code to implement the
> above two integrals is on the order of 40 lines.
>
> I don't know of any connection with pseudospectra techniques. It
> is a straight forward application of contour integrals in the
> complex plane.
>
Indeed. I use something very similar to find the complex frequency
spectrum of self-excited/damped nonradial pulsation modes in stars.
cheers,
Rich
| |
| Gerry Thomas 2004-05-14, 8:31 pm |
|
"Steven G. Kargl" <kargl@troutmask.apl.washington.edu> wrote in message
news:c83ghc$amk$1@nntp6.u.washington.edu...
> In article <3Zapc.81451$FH5.1890949@news20.bellglobal.com>,
> "Gerry Thomas" <gfthomas@sympatico.ca> writes:
>
> Here's the complete citation
>
> Steven G. Kargl and Philip L. Marston, ``Observations and modeling of
> the backscattering of short tone bursts from a spherical shell: Lamb
> wave echoes, glory, and axial reverberations,'' Journal Acoustical
> Society America, vol 85, pp 1014-1028 (1989).
>
> I don't have a PDF version of this paper. The code to implement the
> above two integrals is on the order of 40 lines.
>
> I don't know of any connection with pseudospectra techniques. It
> is a straight forward application of contour integrals in the
> complex plane.
>
Thanks, I'll take a p.
--
E&OE
Ciao,
Gerry T.
| |
| Miguel Ángel Solano Vérez 2004-05-15, 2:31 pm |
| Tank you very much to all of you for the answers.
I will try to solve te problem employing your hints. Perhaps I will contact
to some of you if I need help. O.K.?
Cheers.
Miguel Angel
************
"Rich Townsend" <rhdt@barVOIDtol.udel.edu> escribió en el mensaje
news:c83ivd$n01$1@scrotar.nss.udel.edu...
> Steven G. Kargl wrote:
>
> Indeed. I use something very similar to find the complex frequency
> spectrum of self-excited/damped nonradial pulsation modes in stars.
>
> cheers,
>
> Rich
| |
|
| "Miguel Ángel Solano Vérez" wrote:
>
> problem is not if det(A)=constant or not. Obvioulsy, in a numerical
> computing language (fortran, for example) det(A) is constant for a value of
> beta. But I have to find the value of beta that makes det(A)=0.
In that case you need a function zero finder to solve, f(beta)=0. Look
for "zeroin" in netlib/go dir. btw, "go" stands for Golden Oldies.
> common in electromagnetics: finding zeros of "complex characeristic
> equations".
For finding complex zeroes use Jenkins & Traub "cpoly" in
netlib/toms/419. It's considered a gold (& industry) standard.
| |
| Rich Townsend 2004-05-15, 6:32 pm |
| bv wrote:
> "Miguel Ángel Solano Vérez" wrote:
>
>
>
> In that case you need a function zero finder to solve, f(beta)=0. Look
> for "zeroin" in netlib/go dir. btw, "go" stands for Golden Oldies.
For those who don't know, zeroin is effectively the same as the Brent
root solver appearing in Numerical Recipes -- except, I think that
zeroin is more robust, the NR version being susceptible to instabilities.
>
>
>
>
> For finding complex zeroes use Jenkins & Traub "cpoly" in
> netlib/toms/419. It's considered a gold (& industry) standard.
>
Doesn't cpoly work only with polynomials? Or can it do general complex
root finding?
cheers,
Rich
| |
|
| Rich Townsend wrote:
>
> Doesn't cpoly work only with polynomials? Or can it do general complex
> root finding?
No, it can't do general cases, however, you can easily adapt real root
finding solvers to handle complex roots. e.g. "lmdif" or "lmder" in
netlib/minpack dir.
| |
| Rich Townsend 2004-05-17, 8:31 pm |
| bv wrote:
> Rich Townsend wrote:
>
>
>
> No, it can't do general cases, however, you can easily adapt real root
> finding solvers to handle complex roots. e.g. "lmdif" or "lmder" in
> netlib/minpack dir.
>
>
Apologies for labouring the point, but these routines are for
minimization. While root finding can certainly be recast as
minimization, it is often a pretty inefficient way to approach the
problem, since one often converges to minima which aren't zeros of the
original function.
A suitable root-finding algorithm, which can be generalized to complex
variables, is the standard secant algorithm. However, I'm not sure how
one would go about generalizing a more-sophisticated algorithm, such as
Brent's method.
The problem with complex root finding comes in being able to bracket a
root in the complex plane -- and ensure that it remains bracketed. Steve
Kargl's contour integration approach is one way to deal with this
bracketing problem, since it allows one to constrain the number of roots
within a certain region of the plane.
cheers,
Rich
| |
| Gerry Thomas 2004-05-17, 9:31 pm |
|
"bv" <bvoh@Xsdynamix.com> wrote in message
news:40A94539.A0E11D5E@Xsdynamix.com...
[...]
It's been known for a while: "Ueber unendlich viele Algorithmen zur
Auslosung der Gleichungen," by E. Schroder, Math. Anal. (1870), and
translated by G.W. Stewart. Steve Kargl's approach would appear to be
related to Konig's theorem and is implicit in the use of the psuedospectra
of the det of a matrix. It's a small world afterall!
--
E&OE
Ciao,
Gerry T.
| |
| Steven G. Kargl 2004-05-18, 1:31 am |
| In article <7ccqc.13221$qJ5.352024@news20.bellglobal.com>,
"Gerry Thomas" <gfthomas@sympatico.ca> writes:
>
> It's been known for a while: "Ueber unendlich viele Algorithmen zur
> Auslosung der Gleichungen," by E. Schroder, Math. Anal. (1870), and
> translated by G.W. Stewart. Steve Kargl's approach would appear to be
> related to Konig's theorem and is implicit in the use of the psuedospectra
> of the det of a matrix. It's a small world afterall!
>
I certainly won't claim that the approach I mentioned originated with
me. Indeed, many classic text on Complex Analysis from circa 1950/70
discuss the Winding Number theorem, the calculus of residues, and
Cauchy theorem. I haven't looked for earlier literature, but I suspect
that early atomic/nuclear/quantum physics used a similar approach
(e.g., in Regge Pole analysis). What is amazing to me, is that few
people seem to be aware that one can use contour integrals in the
complex plane to locate zeros a (complicated) functions.
--
Steve
| |
| Gerry Thomas 2004-05-18, 1:31 am |
|
"Steven G. Kargl" <kargl@c-67-168-59-70.client.comcast.net> wrote in
message news:F7gqc.19148$qA.2203458@attbi_s51...
> In article <7ccqc.13221$qJ5.352024@news20.bellglobal.com>,
> "Gerry Thomas" <gfthomas@sympatico.ca> writes:
psuedospectra[color=darkred]
>
> I certainly won't claim that the approach I mentioned originated with
> me. Indeed, many classic text on Complex Analysis from circa 1950/70
> discuss the Winding Number theorem, the calculus of residues, and
> Cauchy theorem. I haven't looked for earlier literature, but I suspect
> that early atomic/nuclear/quantum physics used a similar approach
> (e.g., in Regge Pole analysis). What is amazing to me, is that few
> people seem to be aware that one can use contour integrals in the
> complex plane to locate zeros a (complicated) functions.
>
> --
> Steve
I fully agree, why settle for a lion when you can have a tiger (I recall
this from a musty text by Bromwich, Heaviside, or somesuch)? I wasn't
trying to denigrate your approach (and it's clear to me that you possess
enough ego defenses to know this), but like all good ideas it's difficult
to stake priority claim, and what the heck, to rephrase Henry Ford,
'history is still bunk'.
--
E&OE
Ciao,
Gerry T.
|
|
|
|
|