Newsgroup sci.math.num-analysis 28968

Directory

Maintenance release of MATCOM V2 and MAT
=============================================
The MATCOM V2.01 (The Matlab to C++ Compiler) and MAT V2.01
(C++ Matrix Class Library) maintenance package has been released.
Both MATCOM V2 and MAT can be downloaded freely from our home
page for a 30 days evaluation period.
This maintenance package includes -
o Bug fixes for rare indexing conditions
o Correct fscanf operation when reading more than 4000 elements
o zeros, ones, eye round of non-integer arguments
o colormap function
o "multiply defined" problem with new gcc
o Matrix^(-0.5) fixed
o HP+BLAS, IBM+ESSL problems fixed
o Other minor bug fixes
The maintenance package is free for our licensed users -
a username and password are needed to download it from our
home page, http://www.mathtools.com. If you have
not yet licensed our software, now is the time to do it and
enjoy the introductory prices of 1996.
Tips and Tricks for MATCOM V2 / MAT users
==============================================
o  EasyMex: MEX files made easy for PC/Windows users
   An EasyMex solution is available for PC users. This solution
   creates MEX files without using cmex.bat or depending on the
   WATCOM C++ compiler. This is accessible to all users as
   TechDoc #1028 (http://www.mathtools.com/techdocs/td1028.html).
   If you were having trouble using MEX files on Windows NT,
   try EasyMex.
   (When using EasyMex, be sure to change directory in Matlab to
   the mex directory before running it. Do not run it from the
   path; this might cause GPF.)
   Note: This is for Windows 95 / Windows NT only.
o  Float matrices save memory and time
   Use float matrices in program which does not require more
   precision. Your programs will use half the memory and run
   faster. To compile a program using float matrices:
   matcom myprogram /float
   If you have previously compiled it without /float, add a /force
   switch to force matcom into re-translating:
   matcom myprogram /float /force
o  WATCOM C++ & Visual C++ installation
   If you have problems installing MATCOM V2 for Visual or WATCOM
   compilers, please take a look at the following techdocs:
   Visual C++:  http://www.mathtools.com/techdocs/td1018.html
                http://www.mathtools.com/techdocs/td1027.html
   WATCOM C++:  http://www.mathtools.com/techdocs/td1026.html
 =====================================================================
 MathTools Ltd.                     Web:   http://www.mathtools.com
 P.O.Box 855                      Email:   info@mathtools.com
 Horsham, Pennsylvania              Fax:   1-888-MATHTOOLS (toll free)
 USA 19044-0855                            1-215-957-1719
 =====================================================================

In article <3299C149.3A9659B2@asu.edu>,
Hans D Mittelmann   wrote:
#Jaroslav Stark wrote:
#> 
#> Can anyone point me to efficient ways of computing the derivative of
#> det(A) for non-invertible A, or alternatively an efficient way of
#> calculating the matrix of co-factors of A.
#> 
#> Thus for in genral we have
#> 
#> D det(A) = trace(B.DA)
#> 
#> where B is the transpose of the matrix of co-factors of A. When A is
#> invertible, B is just det(A).A^-1, but what about the genral case?
#> 
#> Answers by e-mail would be appreciated.
#> 
#> J. Stark
#> 
#> E-Mail j.stark@ucl.ac.uk
#Hi,
#here is my theorem. Unless it can be found in the literature, I'd like
#to be quoted as the source.
#
#    D det(a) = sum(i=1,n) product(j.ne.i) lambda_j(A)
#
#This requires one call of, say, the QR algorithm and is thus a O(n^3)
#method.
#
\documentstyle{article}\parindent0pt\parskip1ex\def\dia{{\rm diag}}
\def\bmx#1{\left(\begin{array}{@{}#1@{}}}\def\emx{\end{array}\right)}
\def\trp{^{\rm T}}\def\inv{^{-1}}\def\itp{^{\rm -T}}\def\unv{{\bf e}}
\def\fro{_{\rm F}}
\begin{document}
The problem of calculating the derivative of the determinant of a
matrix $A$ is an interesting one.  As has been pointed out the
derivative with respect to the $(i,j)$-element is the signed cofactor
of that element. (To see this, expand the determinant along the $i$th
row).  When $A$ is nonsingular the signed cofactors are the elements
of $\det(A)A\inv$.  However this approach to computing the derivative
will yield inaccurate results when $A$ is ill conditioned, and the
process breaks down entirely when $\det(A)=0$.
Another possibility is to use the fact that the determinant is the
product of the eigenvalues of $A$. Thus its derivative may be found by
combining derivatives of eigenvalues.  However, this process breaks
down when $A$ has multiple eigenvalues. (The consequences for a zero
multiple eigenvalue are especially nasty\,---\,it is defective
whenever $\det'(A)\neq 0$.)
A better approach is through the singular value decomposition of
$A$.  To fix our notation let
\[
 U\trp A V = \Sigma
\]
where $U$ and $V$ are orthogonal and
$\Sigma=\dia(\sigma_{1},\ldots,\sigma_{n})$ with $\sigma_{1}\geq
\cdots \geq \sigma_{n}\geq 0$.  We will assume that
$\det(U)=\det(V)=1$ so that $\det(A)=\sigma_{1}\cdots\sigma_{n}$.  We
will denote the $j$th columns of $U$ and $V$ by $u_{j}$ and $v_{j}$.
First consider the general problem of assessing the effects of a
perturbation $A+E$ on the determinant.  We have
\[
\det(A+E) = \det(U\trp AV + U\trp EV) \equiv \det(\Sigma + F),
\]
where the $(i,j)$-element of $F$ is $u_{i}\trp Ev_{j}$.
Now the derivative of $\det(\Sigma)$ with respect to its $(i,j)$-element
is zero if $i\neq j$ and is
\begin{equation}
\label{deldef}
\delta_{i} = \prod_{k\neq i}\sigma_{k}
\end{equation}
if $i=j$.  Consequently, up to first order terms
\begin{equation}
\label{dape}
\det(A+E) \cong \sum_{k=1}^{n} (u_{k}\trp Ev_{k})\delta_{k}.
\end{equation}
To determine the derivative of $\det(A)$ with respect to the
$(i,j)$-element of $A$, let $E = \epsilon\unv_{i}\unv_{j}\trp$, where
$\unv_{i}$ and $\unv_{j}$ are the $i$th and $j$th unit vectors.  This
represents a perturbation of $\epsilon$ in the $(i,j)$-element of $A$.
By (\ref{dape}) the corresponding perturbation in $\det(A)$ is
$\epsilon\sum_{k=1}^{n} u_{ik}v_{jk}\delta_{i}$.  Hence
\[
\frac{\partial \det(A)}{\partial a_{ij}} =
    \sum_{k=1}^{n} u_{ik}v_{jk}\delta_{k},
\]
where $\delta_{k}$ is defined by (\ref{deldef}).  In particular if
$\det(A)=0$, then $\sigma_{n}=0$ and
\[
\frac{\partial \det(A)}{\partial a_{ij}} = u_{in}v_{jn}\delta_{n}.
\]
It should be noted that the above approach works with any orthogonal
reduction of $A$ to a form where its cofactors can be determined by
inspection.  Bidiagonal form is a possibility, but in this case the
expression for the derivative has $O(n^{2})$ terms.
Hope this helps.
Pete Stewart\\
University of Maryland\\
College Park\\
Dec  1 1996
\end{document}

Hi.
I wonder if there is anyone who have a speedy solution to the
following problem;
Integrate the function exp(-a/sin(x)) for x between 0 and pi/2.
The parameter a is a positive constant.
So far I have been using a self-made trivial solution, a sort
of an adaptive Simpson approximation. But is there a faster
and more reliable way ?
Best regards
Peter Jansson
mailto:jansson@tsl.uu.se

Well, as I mailed to Mister McClure, a 'complex' 2x2 matrix is simply
         l a -b l
         l b  a l
which is another way to mean a+ib.
-- 
M.TIBOUCHI
>See you soon on the Global Village

In article (Dans l'article) <57n8nk$tq@news1.mnsinc.com>, Manuel Alfonso
 wrote (écrivait) :
> Given the following:
> 
> -An ellipse centered at the origin
> -A radius R1 from the origin to some point A on the ellipse (R1 length known)
> -A radius R2 from the origin to some point B on the ellipse (R2 length
unknown)
> -The arclength from A to B of the ellipse
> 
> How would one find the angle between R1 and R2?
> 
> Please e-mail response(s) to malfonso@logicon.com.
> 
> Thanks
Just an element of answer : you should maybe consider the polar equation
of the ellipse or the parametric one.
Say you use parametric representation.
x = a sin (omega1 t + phi1)
y = b sin (omega2 t + phi2) is an ellipse centered on origin of which
you've got the precise coords for each point, and also the radius (sqrt
(x^2+y^2)).
Working on such equalities, you'll maybe find out.
-- 
M.TIBOUCHI
>See you soon on the Global Village

I heard about the exponential and logarithm of complex numbers.
They said a complex num has infinetely many logarithms, unless it's real.
Has anyone infos on that ?
-- 
M.TIBOUCHI
>See you soon on the Global Village

In article (Dans l'article) <329F248E.1350@asci.fr>, Jeanfaivre Laurent
 wrote (écrivait) :
> Bonjour a tous,
> 
> Je recherche des paradoxes scientifiques ! Je m'explique : la science
> permet souvent d'etablir des resultats qui defient completement le sens
> physique voire le bon sens commun. Par exemple, quel non scientifique
> irait s'imaginer qu'un avion est plus petit en vol qu'au sol ?
> Par consequent, quel que soit votre domaine, si vous connaissez de pres
> ou de loin, des "phenomenes", quels qu'ils soient, qui peuvent paraitre
> surprenants au commun des mortels, je vous remercie de me les decrire.
> Je recherche d'autre part des pieds-de-nez scientifiques, du type de la
> demonstration mathematique (truquee heureusement) de 2=1.
> 
> N'hesitez pas a m'envoyer vos suggestions, tout m'interesse ...
> Merci d'avance.
Faites gaffe. Ecrire en anglais dans un Newsgroup francais, c'est deja
dangereux, alors le contraire, c'est prendre le risque de ne pas obtenir
de reponse! ;-)
Pour les paradoxes, je peux vous proposer un tuyaux prolifique : la
theorie des ensembles. Bien que je ne sois pas un eminents mathematicien
(juste un collegien de 13 ans), je sais qu'on peut y trouver de petites
merveilles : "les ensembles paradoxaux".
Un exemple : "Soit E l'ensemble des personnes d'une ville qui se taillent
elles-meme la barbe. Le barbier de cette ville coupe la barbe a toutes les
personnes ne faisant pas partie de E (c-a-d qui ne se taille pas la
barbe). En fait-il lui-meme partie ?"
Amusez-vous bien!
A+
-- 
M.TIBOUCHI
>See you soon on the Global Village

In article  edb@syd.dms.csiro.au "Ed Breen" writes:
>>It is consistent in that passing an apparent actual array arrgument and
>>takeing a formal array parameter both imply that the argument/parameter
>>is really a pointer. How would you propose to change this without making
>>a *major* change to the language to make arrays first class objects and
>>passable directly as function arguments?
>
>I am not certain what, if any, type of backlash the following
>options will generate, but here goes.
>
>I have no objections to arrays and functions being 2nd class
>citizens. Although, given that structures and unions have a higher
>status, I would have thought, that it would be a simple matter to raise
>arrays to the same level.
>
>However, the following rules should work, while keeping arrays and
>functions in 2nd class:
>
>        A formal array parameter will be treated as its specified type
>when used as an operand to the sizeof operator.
As things stand sizeof tells you the size of the actual argument passed.
In your proposal the size needn't correspond to anything meaningful
in the program (e.g. the actual argument or the array pointed to by the
actual argument). There is also a problem with something like:
int main(int argc, char *argv[])
since char *[] is a type with no specified size. However the parameter
still exists as an object of type char ** and with the corresponding size.
For consistency shouldn't sizeof(*&x;) be the same as sizeof(x) where x
is an lvalue?
>or even:
>
>        A formal parameter will be treated as its specified type, but
>arrays and functions will be passed by reference only.
As things stand 'array' parameters are really pointers and, as such,
can be assigned to. What you are suggesting would prevent that. It also
means a major new concept, pass by reference, be added to the language.
>>Arrays always "evaluate" to a pointer to their first element (& and
>>sizeof don't take the value of their arguments, they take the address and
>>size respectively). 
>
>As this information can be obtained from each parameter's
>specification, why lose it? Either be honest with the programmer and
>make such parameter declarations illegal or accept what the programmer
>has declared as the true type of the parameter.
The problem is that such a parameter type has no use since an array
type could never be passed - any array 'evaluates' to a pointer to
its first element. Also array access is typically handled through pointers
so a pointer may be all that is available, for example with malloc'd
arrays.
Eliminating the rewrite rule (at least for array parameters, and making
array parameter types illegal) is a possibility (except that it breaks
existing code). I could certainly live without it when writing new code.
>The intention of the above I believe, is to ensure that the scope of
>the parameters are the same as the identifiers declared in the top
>level of the function. I see no problems with a parameter retaining
>its original specification. Why should information about a type be
>lost? FAPP, the programmer should not be concerned about the rewriting
>rule, all she needs to be concerned with is that arrays and functions
>are passed, in essence, by reference only.
Except, as I noted, references to C arrays are often in the form of
pointers to an element of the array. This isn't a suitable value to
pass as a reference function argument. Functions where you could only
pass a genuine array lvalue wouldn't be very useful.
-- 
-----------------------------------------
Lawrence Kirby | fred@genesis.demon.co.uk
Wilts, England | 70734.126@compuserve.com
-----------------------------------------

Could someone please comment on whether the following is a (theoretically) 
valid method of computing the Lobatto points of degree n+1?
1.  The Lobatto points of degree 3 are {-1,0,1}.
2.  If x and y are consecutive Lobatto points of degree n, then 
    there is exactly one Lobatto point `z' of degree n+1 in the 
    interval [x,y].
3.  Each Lobatto point is a simple root of P'_n, so P'n changes sign
    at z.
4.  The secant method can be used with initial endpoints [x,y] to 
    approximate z.
5.  The other two Lobatto points are -1 and 1.
This method seems to work, but I have been unable to justify it.
I've also come across an algorithm that uses Newton's method to find
each Lobatto point with initial guess of cos(j*Pi/n), j=1..n-1, but
have not found justification for that, either.
Can anybody point me towards some references regarding such justification?
Thank you.
  ==Jake

In article , 
medtib@club-internet.fr says...
>
>I heard about the exponential and logarithm of complex numbers.
>
>They said a complex num has infinetely many logarithms, unless it's 
real.
>
>Has anyone infos on that ?
>
>-- 
>M.TIBOUCHI
That depends.  Abramowitz and Stegun (Handbook of Mathematical Functions) 
define the "logarithm function" "ln(z)" as a single valued function whose 
imaginary part is in the range -pi < y <= pi, but the "general logarithm 
function" "Ln(z)" is the corresponding many valued function whose 
imaginary part differs from that of ln(z) by an integer multiple of 2*pi 
(even when z has no imaginary part). 
ln(z) is the principal branch of Ln(z), and is discontinuous over the 
negative real axis, while Ln(z) is continuous everywhere (except the 
origin).
See Abramowitz and Stegun for more details.
Damien Holloway

Hi!
I'm looking for FFT procedures (C or Fortran) with these
characteristics:
    Real to complex (and vice versa)
    3-dimensional
    Any NX, NY, NZ possible 
I'm using the SGI complib.sgimath library procedure "sfft3du",
which exactly matches my needs. But, of course, the source code
is not available.
Any help in finding source code to replace "sfft3du" is very
much appreciated!
Ralf

Please help me understand the following homework problems.
I need to find the standard coordinatization of 1+x+x^3 in P3.  The answer
at the back of the book says:
                                [1]
                                [1]
                                [0]
                                [1]
basically a 4x1 matrix.  My question is how can there be a 4x1 matrix in a
3 dimensional space.
-- 
Robert Gelb
Senior Systems Analyst
Data Express
Garden Grove, California USA
(714)895-8832

We are looking for everything dealing with mathematics concerning billards.
We are especially interested in trajectories in rectangular billards.
We also look for informations about approximating transcendant and irrationnal 
numbers with continuous fractions and 
also about ergodic theorem (Poincare's theorem)

Robert Gelb wrote:
> 
> Could someone please explain the Gram-Schmidt Process in regards to the
> orthogonalization process (inner spaces)?
> 
> Thanks in advance.
> 
> --
> Robert Gelb
> Senior Systems Analyst
> Data Express
> Garden Grove, California USA
> (714)895-8832
As a Mathematics graduate from CSULB, I seem to remember that the Math
Department teaches a class on Approximation Theory. 
1) Look up who teaches it and get them to explain it to you.
2) OR check with Dr. Warner or Dr. Maryfield in the Math Dept.
3) OR check out from the library:
   Introduction to Numberical Analysis by F.B. Hildebrand (QA297.H54 
   1974) and read Chap 7 Least-Squares Polynomial Approximation and
   note 7.13 Gram Approximation. There are other books, but this should
   get you started.
I hope this is of some help.
Gordon Talge
-- 
+------------------------------------------------------------+
|  Gordon Talge WB6YKK            e-mail: gtalge@pe.net      |
|  Department of Mathematics      QTH: Loma Linda, CA        | 
|  Mt. San Jacinto College        Lat.  N  34° 03' 03.6720"  | 
|  San Jacinto, CA                Long. W 117° 15' 09.5760"  |
+------------------------------------------------------------+

In article <329F1EEC.41C6@Sorry.Trying.to.avoid.spam>,
Arjendu Pattanayak   wrote:
>Given that I know the value of the integrals (all notation in TeX):
>
>\int f(x) dx
>
>and
>
>\int x^2 f(x) dx
>
>what can I say about the integral
>
>\int exp(-d.x^2)f(x) dx
>
>where d is some positive constant and the range of all the integrals are
>from -\infty to \infty ?
Nothing much, unless you know that f >= 0.  In that case obviously 
int(-d.x^2) f(x) dx <= int f(x) dx 
but that's all you can say (you can approach this bound arbitrarily
closely with f(x) mostly concentrated near x=0, with a little bit
far away to make int x^2 f(x) dx whatever you want).
Robert Israel                            israel@math.ubc.ca
Department of Mathematics             (604) 822-3629
University of British Columbia            fax 822-6074
Vancouver, BC, Canada V6T 1Y4

In article (Dans l'article) <329F1EEC.41C6@Sorry.Trying.to.avoid.spam>,
Arjendu Pattanayak  wrote (écrivait) :
> what can I say about the integral
> 
> \int exp(-d.x^2)f(x) dx
> 
> where d is some positive constant and the range of all the integrals are
> from -\infty to \infty ?
> 
Just one remark. I'm not sure but if I can remember, \int exp(-x^2) isn't
writable in terms of algebraic functions. I think it's Ei (x).
-- 
M.TIBOUCHI
>See you soon on the Global Village

In article <57i7vi$42d@lastactionhero.rs.itd.umich.edu>,
Hyun Min Peck  wrote:
> In article <576rfd$lqg@electra.saaf.se>, pausch@electra.saaf.se says...
> 
>>Even single-dimensional arrays are "second-class citizens" in C
>>-- try this example:
>> 
>>    void sub( int arr[10] )
>>    {
>>        int i;
>>        for( i = 0; i < sizeof(arr)/sizeof(arr[0]); i++ )
>>            printf( " %d", arr[i] );
>>    }
>> 
>>    main()
>>    {
>>        int arr[10] = { 1,2,3,4,5,6,7,8,9,10 };
>>        sub( arr );
>>    }
>> 
>>This program does not output  " 1 2 3 4 5 6 7 8 9 10" as one naively
>>might expect.  Instead it outputs " 1" or " 1 2", depending on the
>>environment used.....
> 
> The reason why you are getting the results is that 
> sizeof (arr)/sizeof (arr[0]) = 1
Not always -- on some compilers sizeof(int*)/sizeof(int) will be 2, not 1...
> not 10 as you might expect.
> Try instead i<10, You will get the output "1 2 3 4 5 6 7 8 8 10".
Of course..... now why does C treat array parameters this way?
You could also try:
    typedef struct { int arr[10]; } ARRAY;
    void sub( ARRAY a )
    {
        int i;
        for( i = 0; i < sizeof(a.arr)/sizeof(a.arr[0]); i++ )
            printf( " %d", a.arr[i] );
    }
    main()
    {
        ARRAY a = { {1,2,3,4,5,6,7,8,9,10} };
        sub( a );
    }
Now you'll get the output   1 2 3 4 5 6 7 8 9 10     ......
-- 
----------------------------------------------------------------
Paul Schlyter,  Swedish Amateur Astronomer's Society (SAAF)
Grev Turegatan 40,  S-114 38 Stockholm,  SWEDEN
e-mail:  pausch@saaf.se     psr@home.ausys.se    paul@inorbit.com

Please help me with my schoolwork.  I don't understand the example below
out of my school book.
This example is on Linear Transformation and 
Matrix Representation.
***************************************************
Example 8. Chapter 5.3, page 330
"First Course In Linear Algebra" by Moore, Yaqub
Let F:P2 -> P2 be the linear transformation
	F(c+bx+ax^2) = (c-a) + (a+b)x + (2c-b-3a)x^2
The standard ordered basis for P2 is {1,x,x^2}, so the standard matrix
representation of F is found by computing F(1)=1+2x^2; F(x)=x-x^2; and
F(x^2)=-1+x-3x^2.  
***************************************************
I don't understand how F(1), F(x), F(x^2) were found.  Can someone
explain, how these were found?
 -- 
Robert Gelb
Senior Systems Analyst
Data Express
Garden Grove, California USA
(714)895-8832