Newsgroup sci.math.num-analysis 28553

Directory

Subject: Re: discrete hankel transforms -- From: mathar@qtp.ufl.edu (Richard Mathar)
Subject: Using Fortran-77 to solve IVPs -- From: "Ritesh Patel"
Subject: Kernel Regression -- From: Mark Leung
Subject: eigenvalues of ill-cond matrices -- From: alan koivunen
Subject: mesh generation program needd -- From: bihu@newstand.syr.edu (Bin Hu)
Subject: Re: Directed rounding on the Pentium -- From: "Jeffery J. Leader"
Subject: Directed rounding on the Pentium -- From: "Stephen W. Hiemstra"
Subject: Re: (1/n)! as integral expression -- From: sven@milliway.shnet.org (Sven Koehler)
Subject: leading dimension -- From: Zoltan Bus
Subject: Re: PLEASE HELP ME: I need to factor/simplify an algebraic expression. -- From: John Hudson
Subject: Re: eigenvalues of ill-cond matrices -- From: Paul Abbott
Subject: Re: Using C for number-crunching (was: Numerical solution to -- From: tydeman@tybor.com
Subject: Re: Directed rounding on the Pentium -- From: checker@netcom.com (Chris Hecker)
Subject: Re: mesh generation program neededf -- From: robert@informatik.rwth-aachen.de (Robert Schneiders)
Subject: Division by 2 -- From: John Hench
Subject: Curve through tree points -- From: Ricky Helgesson
Subject: Re: Curve through tree points -- From: kfroenigk@mmm.com (Karl F. Roenigk)
Subject: Re: Curve through tree points -- From: jeppe@dfi.aau.dk (Jeppe Burchhardt)
Subject: Re: Directed rounding on the Pentium -- From: n8tm@aol.com
Subject: Re: Directed rounding on the Pentium -- From: "Stephen W. Hiemstra"
Subject: Re: Division by 2 -- From: n8tm@aol.com
Subject: Accurate exp() -- From: n8tm@aol.com
Subject: Re: Curve through tree points -- From: jeppe@dfi.aau.dk (Jeppe Burchhardt)
Subject: Re: how do you generate a sine wave with simple adds and subtracts. -- From: Stan Armstrong
Subject: Re: Need help with integration -- From: Stan Armstrong
Subject: Re: leading dimension -- From: Konrad Hinsen
Subject: Re: eigenvalues of ill-cond matrices -- From: hwolkowi@orion.math.uwaterloo.ca (Henry Wolkowicz)
Subject: numerical integration of non linear equ -- From: cooder
Subject: Re: [Q] Need to plot in spherical coords. -- From: dwells@nrao.edu (Don Wells)
Subject: Need your help -- From: "Brad Marts"
Subject: complex Newton's method -- From: trscavo@forbin.syr.edu (Tom Scavo)
Subject: Quadratic constrained LS minimisation -- From: Moritz Harteneck
Subject: Re: Accurate exp() -- From: "Dann Corbit"
Subject: Re: Directed rounding on the Pentium -- From: "Stephen W. Hiemstra"
Subject: running FDTD fortran codes on SGI -- From: Tim Black
Subject: Re: Need your help -- From: "Brad Marts"
Subject: Curl and divergence estimation from irregularly spaced data -- From: Olaf Amm
Subject: Re: Directed rounding on the Pentium -- From: "Jeffery J. Leader"
Subject: Re: Directed rounding on the Pentium -- From: "Jeffery J. Leader"
Subject: Re: [Q] Using pseudoinverse in Bayes discriminant function? -- From: Greg Heath

Articles

Subject: Re: discrete hankel transforms
From: mathar@qtp.ufl.edu (Richard Mathar)
Date: Sun, 10 Nov 1996 14:20:50 EST
In article <5579s1$kb5@scotsman.ed.ac.uk>, msh@holyrood.ed.ac.uk ( Mark Higgins) writes:
|> I currently have the hankel transform that uses lagged and related
|> convolutions form NETLIB, (TOMS/588) but this seems to rely on using a
|> external genarating funtion for the data to be transformed.  I can
|> create a function wich fitted some spline to my data thus simulating a
|> continuons function to be transformed
|> 
|> My questions are --> 
|> 
|> anyone have experience of doing this for discrete
|> data, and does it work ok ?
|> 
|> and anyone have a discrete hankel transform that will work from a
|> discrete series ?
|> 
Because the Hankel transform is a convolution between 1/x and f(x),
the discrete version would use the multiplication of the FTs [that
is a function proportional to a theta-function times the FT of f(x)]
and transform them back. Mapping the convolution onto two FTs with
an intermediate truncation seems to be the most natural way, and avoids
any difficulties from the pole in 1/x .
cf: T Hasegawa, T Torii, I Ninomiya, "Generalized Chebyshev
Interpolation and Its Application to Automatic Quadrature", Math
Comp 41 (164) (1983) 537-53
T Hasegawa T Torii, "An automatic quadrature for Cauchy Principal value
integrals", ibida 56 (194) (1991) 741-54
P Luchini, "Fourier analysis of numerical integration formulae",
Comp Phys Comm 83 (2&3) (1994) 227-35
-- 
mathar@qtp.ufl.edu
Return to Top
Subject: Using Fortran-77 to solve IVPs
From: "Ritesh Patel"
Date: 10 Nov 1996 21:22:23 GMT
Hi all,
I need some help with Fortran-77 code to solve 4th order Runge-Kutta
method, Euler's method, modified Euler's method, Search and bound method,
and general IVPs.
Any help would be greatly appreciated.
Regards,
Jeenal Patel.
Email :   jeenal@phase.u-net.com
        OR k951066@crystal.kingston.ac.uk
Return to Top
Subject: Kernel Regression
From: Mark Leung
Date: Sun, 10 Nov 1996 18:01:30 -0500
Does anyone know if there is any computer package for
running kernel regression (or parzen window estimation)?
Thanks.
Mark
Return to Top
Subject: eigenvalues of ill-cond matrices
From: alan koivunen
Date: Sun, 10 Nov 1996 23:31:21 -0500
Hi - 
	I'm trying to find eigenvalues and eigenvectors of some ill-conditioned
	matrices. They are real, symmetric Toeplitz matrices, an example of which
	is coded below:
	do i=1,N
	  do j=i,N
	    A(i,j) = exp(-2.0*(j-i)*(j-i)/((N-1)*(N-1)))
	    A(j,i) = A(i,j)
	  enddo
	enddo
	With N=13, the eigenvalue magnitudes range in value from about 
	10^(-16) to about 10. The eigenvalues should be real and positive,
	but using MATLAB's eig() function and IMSL's DEVCSF routine produces
	negative (small) eigenvalues - so things aren't quite right.
	So, is there anyone out there who can suggest a way to get 
	eigenvalues/vectors for these troublesome matrices ? Maybe some
	sort of iterative refinement of these initial estimates ?
	Any suggestions are appreciated.
		Thanks
		Al Koivunen
Return to Top
Subject: mesh generation program needd
From: bihu@newstand.syr.edu (Bin Hu)
Date: 11 Nov 1996 05:02:00 GMT
Hi
I am doing some work on FEM problem. However, I donot have a good program
on 3-D mesh generation. CAn anyone recomend one?
thanks
Bin
Return to Top
Subject: Re: Directed rounding on the Pentium
From: "Jeffery J. Leader"
Date: Sun, 10 Nov 1996 21:14:38 -0800
Kahan gave an interesting talk on this type of thing at the Gragg
Conference in Monterey last weekend--all the directed rounding,
overflow/underflow, etc. stuff that's built into the chips but not
reasonably accessible from most compilers.  Interesting but frustrating.
-- 
I turn away with fright and horror from this lamentable evil of 
functions without derivatives.
 -Charles Hermite
Return to Top
Subject: Directed rounding on the Pentium
From: "Stephen W. Hiemstra"
Date: 11 Nov 1996 04:32:59 GMT
I have implemented a series of math functions which use the directed
rounding commands available for the floating point unit (fpu) on the
Pentium.  In checking through these functions, I divided 1 by 3 rounding up
then rounding down to compare the results.  In both cases, I obtained the
same result (.333333 out to the number of places that can be observed in my
debugger (Visual C++ 4.2) which presumably is the full representation of
real10 precision.  (My routines are actually written in MASM 6.11d
assembler).
While I might have expected the same result for numbers outputted to the
double precision variables (they were also the same value, but not carried
out to as many places), I expected to see a 2 in the last place for
rounding down and a 3 in the last place for rounding up.  Is this the
correct expectation?  Have I errored in my computations or does the FPU
simply not display the last digit being rounded?
Stephen
Return to Top
Subject: Re: (1/n)! as integral expression
From: sven@milliway.shnet.org (Sven Koehler)
Date: 07 Nov 1996 13:32:00 +0100
The last integral in this message may easily be evaluated by substituting
\sqrt[n]{t} for x and remembering the well-known functional equation
\Gamma(x+1)=x\Gamma(x).This approach seems to be much more perspicious  
than yours.
Return to Top
Subject: leading dimension
From: Zoltan Bus
Date: Mon, 11 Nov 1996 09:35:57 +0100
Hi!
Somebody, please, define me what is the meaning of the expression
'leading dimension' which can be found in almost FORTRAN matrix
subroutine. Excuse me for this silly question.
Thank you very much
Zoltan Bus
Return to Top
Subject: Re: PLEASE HELP ME: I need to factor/simplify an algebraic expression.
From: John Hudson
Date: Mon, 11 Nov 1996 09:29:19 +0000
Tai wrote:
> 
> Hi,
> 
> I need to factor/simplify: (a^5 - b^5)/(a^3 - b^3)
> 
Being a matter of general and vital interest I typed this into MathCad 
6 and simplified it symbolically and it came back with
a^5 - b^5     =     a^4 + b.a^3 + b^2.a^2 + b^3.a + b^4
_________           ___________________________________
a^3 - b^3                 a^2 + b.a + b^2
so it has found one common factor which appears to be a-b
Return to Top
Subject: Re: eigenvalues of ill-cond matrices
From: Paul Abbott
Date: Mon, 11 Nov 1996 18:58:28 +0800
alan koivunen wrote:
>I'm trying to find eigenvalues and eigenvectors of some ill-conditioned
>matrices. They are real, symmetric Toeplitz matrices, an example of 
>which is coded below:
> 
>         do i=1,N
>           do j=i,N
>             A(i,j) = exp(-2.0*(j-i)*(j-i)/((N-1)*(N-1)))
>             A(j,i) = A(i,j)
>           enddo
>         enddo
> 
>With N=13, the eigenvalue magnitudes range in value from about
>10^(-16) to about 10. The eigenvalues should be real and positive,
>but using MATLAB's eig() function and IMSL's DEVCSF routine produces
>negative (small) eigenvalues - so things aren't quite right.
> 
>So, is there anyone out there who can suggest a way to get
>eigenvalues/vectors for these troublesome matrices ? 
Any of the Computer Algebra systems should be able to do this without
any trouble.  Here is a solution using Mathematica 3.0:
In[1]:= Toeplitz[n_]:=Table[Exp[-2(j-i)^2/(n-1)^2], {i, n}, {j, n}]
In[2]:= Timing[Eigenvalues[N[Toeplitz[13], 35]]]
Out[2]= {1.8 Second, 
  {9.714289658177802948116981413446968039217238874`35, 
   2.8087047226601903588888427367223952406986690664`35, 
   0.4318113387746175560652633047322733083888887423`35, 
   0.04211529596485480041989551847174304301749342335`35, 
   0.0029201614706262355299008707132439452221134816464`35, 
   0.00015247429079448480011310165012714444499654618885`35, 
   6.150970244115566358292942430866977377286994027`35*^-6, 
   1.9294232139272084071757708173097648232541969225`35*^-7, 
   4.662923268040823404382745324585306129887057093`35*^-9, 
   8.452676823084502937530316623177950895236479696`35*^-11, 
   1.0890952782993752514438522541130430812858525505`35*^-12, 
   8.941025867763620061655233278333984110022684667`35*^-15, 
   3.531596847111418268275223468242875950227707268`35*^-17}}
Cheers,
	Paul 
_________________________________________________________________ 
Paul Abbott
Department of Physics                       Phone: +61-9-380-2734 
The University of Western Australia           Fax: +61-9-380-1014
Nedlands WA  6907                         paul@physics.uwa.edu.au 
AUSTRALIA                           http://www.pd.uwa.edu.au/Paul
          God IS a weakly left-handed dice player
_________________________________________________________________
Return to Top
Subject: Re: Using C for number-crunching (was: Numerical solution to
From: tydeman@tybor.com
Date: 11 Nov 1996 10:25:36 GMT
>: First, I recently posted a disclaimer to two of my earlier postings,
>@> but now I want to disclaim the disclaimer. :-) Someone had posted
>@> saying that const qualification (in C) does not convey non-alias
>@> information, because the const qualifier guarantees only that no
>@> attempt will be made to write to the underlying addresses *through the
>@> const-qualified pointer*.
C9X, the revision of the C language that we (the C standards committee)
hope to have done by year end 1999, has added a new type qualifier: restrict.
It applies to pointers and indicates that there are no aliases to the 
pointed at object (unless the programmer gives up that property by doing
an assignment of the restricted pointer to another pointer).  In effect,
a restricted pointer acts like a pointer to malloc'ed memory.  A good place
to use restricted pointers is the parameters to a function.  The compiler
vendors who have already implemented this new feature report great performance
improvements.
Fred Tydeman          +49 (7031) 288-964  Tydeman Consulting
Meisenweg 20          tydeman@tybor.com   Programming, testing, C/C++ training
D-71032 Boeblingen                        Voting member of X3J11 (ANSI "C")
Germany               Sample FPCE tests:  ftp://ftp.netcom.com/pub/ty/tydeman
Return to Top
Subject: Re: Directed rounding on the Pentium
From: checker@netcom.com (Chris Hecker)
Date: Mon, 11 Nov 1996 10:27:26 GMT
"Jeffery J. Leader"  writes:
>Kahan gave an interesting talk on this type of thing at the Gragg
>Conference in Monterey last weekend--all the directed rounding,
>overflow/underflow, etc. stuff that's built into the chips but not
>reasonably accessible from most compilers.  Interesting but frustrating.
That sounds great.  Is the paper on the web anywhere?
Chris
Return to Top
Subject: Re: mesh generation program neededf
From: robert@informatik.rwth-aachen.de (Robert Schneiders)
Date: 11 Nov 1996 09:51:18 GMT
Hi, 
the page 
http://www-users.informatik.rwth-aachen.de/~roberts/meshgeneration.html
has a section on meshing software. 
Hope it helps 
 Robert 
Robert Schneiders
Lehrstuhl fuer Angewandte Mathematik, insb. Informatik 
RWTH Aachen 
Ahornstr. 55
52056 Aachen 
Germany 
email: robert@feanor.informatik.rwth-aachen.de 
WWW: http://www-users.informatik.rwth-aachen.de/~roberts/
Tel.: +241-804558
Fax: +241-8888215
Return to Top
Subject: Division by 2
From: John Hench
Date: Mon, 11 Nov 1996 11:53:09 +0000
I have a quick question:  Do most computers treat
the division of a (real) floating point number by 
2 differently than division by an arbitrary scalar?
It seems to me more related to addition to me, 
since the exponent of the floating point number 
is decremented.  I guess what I'm driving at is
this: is it faster to divide by two than to 
divide by an arbitrary scalar?  I suppose that
every computer is different, but I'm just trying
to get a feel for this...
Thanks, John
------------------------------------------------
 Dr. J.J. Hench  
 Dept. of Mathematics, Univ. of Reading, England   
 Institute of Informatics and Automation, Prague
-------------------------------------------------
Return to Top
Subject: Curve through tree points
From: Ricky Helgesson
Date: Mon, 11 Nov 1996 13:12:38 +0100
I will draw a curve through three points.
The three coordinates are "indata" and whish to get the equation for the
curve through these points so that I can plot that curve.
The equations are (of course):
y1 = a*(x1)^2 + b*x1 + c
y2 = a*(x2)^2 + b*x2 + c
y3 = a*(x3)^2 + b*x3 + c
The thing is that I have no program for calculating this, like matlab. I
know that in matlab, this is not too hard, so my question is:
Could some good person out there help me calculating this curve or give
me a hint where I can find matlab for Windows 3.11/95/NT?
Thanks in advance,
Ricky Helgesson
Com In Center AB
SWEDEN
Return to Top
Subject: Re: Curve through tree points
From: kfroenigk@mmm.com (Karl F. Roenigk)
Date: Mon, 11 Nov 1996 06:57:54 -0600
Ricky Helgesson wrote:
> 
> I will draw a curve through three points.
> 
> The three coordinates are "indata" and whish to get the equation for the
> curve through these points so that I can plot that curve.
> 
> The equations are (of course):
> 
> y1 = a*(x1)^2 + b*x1 + c
> y2 = a*(x2)^2 + b*x2 + c
> y3 = a*(x3)^2 + b*x3 + c
> 
> The thing is that I have no program for calculating this, like matlab. I
> know that in matlab, this is not too hard, so my question is:
> 
> Could some good person out there help me calculating this curve or give
> me a hint where I can find matlab for Windows 3.11/95/NT?
> 
> Thanks in advance,
> 
> Ricky Helgesson
> Com In Center AB
> SWEDEN
You might express your independent variable relative to the curve at
point 1, this reduces the problem to 2 unknowns in an equation as a
function now of a transformed independent variable w=xi-x1.  Then you
have:
y2-y1=a'*(x2-x1)^2+b'*(x2-x1)=a'*w2^2+b'*w2=r1
y3-y1=a'*(x3-x1)^2+b'*(x3-x1)=a'*w3^2+b'*w3=r2
which is easily solved by ratio of determinants. e.g.
Au=r has a solution if |A|*0,
u(i)=|A(i)|/|A|
where |A(i)| is the determinant of the coefficient matrix having the ith
column replaced by the right hand side vector, and u(i) is the ith
unknown variable of interest.  This is a convenient way to solve in
closed form any nonsingular 2 or 3 dimension problem.
Return to Top
Subject: Re: Curve through tree points
From: jeppe@dfi.aau.dk (Jeppe Burchhardt)
Date: 11 Nov 1996 13:13:22 GMT
Ricky Helgesson (rickyhelgesson@hotmail.com) wrote:
: I will draw a curve through three points.
: The three coordinates are "indata" and whish to get the equation for the
: curve through these points so that I can plot that curve.
: The equations are (of course):
: y1 = a*(x1)^2 + b*x1 + c
: y2 = a*(x2)^2 + b*x2 + c
: y3 = a*(x3)^2 + b*x3 + c
: The thing is that I have no program for calculating this, like matlab. I
: know that in matlab, this is not too hard, so my question is:
: Could some good person out there help me calculating this curve or give
: me a hint where I can find matlab for Windows 3.11/95/NT?
: Thanks in advance,
: Ricky Helgesson
: Com In Center AB
: SWEDEN
--
----
Jeppe Burchhardt, Ph.D.
Centre of Magnetic Resonance, Aarhus University Hospital, Denmark
Return to Top
Subject: Re: Directed rounding on the Pentium
From: n8tm@aol.com
Date: 11 Nov 1996 13:42:24 GMT
I expect that it would be necessary to display the entire result (in
hexadecimal, for instance) in order to check the effects of directed
rounding.  The changes in the low order bit will not display neatly in
decimal, although they will show up if enough places are printed (17 for
the usual double).
Tim
Return to Top
Subject: Re: Directed rounding on the Pentium
From: "Stephen W. Hiemstra"
Date: 11 Nov 1996 13:31:07 GMT
Jeffery,
Thanks for the note.  I have seen references to this problem that refer to
high level languages like FORTRAN and Pascal.  However, the functions are
accessible in assembler and seem to work for the most part.  I am able, for
example, to record when numerical exceptions occur and flip the right bits
using the Pentium control word.  Being able to observe the actual outcomes
as they occur may, however, be a different story.
Stephen
Jeffery J. Leader  wrote in article
<3286B63E.5BB6@WorldNet.att.net>...
> Kahan gave an interesting talk on this type of thing at the Gragg
> Conference in Monterey last weekend--all the directed rounding,
> overflow/underflow, etc. stuff that's built into the chips but not
> reasonably accessible from most compilers.  Interesting but frustrating.
> 
Return to Top
Subject: Re: Division by 2
From: n8tm@aol.com
Date: 11 Nov 1996 13:52:55 GMT
On most modern RISC architectures, there is no special way to divide by 2.
As there are separate registers for floating point and integers, there is
no quick way to short cut. However, I fixed up my copy of gcc for the
68020 so that it recognized such cases and used the special instruction
which scales by powers of 2.  Most compilers will recognize cases where
division can be changed to multiplication without changing the result.  On
the VAX architecture, there was a small list of special constants,
including .5 and 2, which could be buried in the instruction, and their
compilers recognized those.  Some of the older architectures, including
Harris 24-bit and GE 600/Honeywell 6000 had an exponent add instruction
which was used for this, but had to include a separate check to skip the
operation if it would cause underflow.
So, you're right, the answer to this question is machine-dependent and
historical.
Tim
Return to Top
Subject: Accurate exp()
From: n8tm@aol.com
Date: 11 Nov 1996 13:57:57 GMT
I have nearly completed revision of my libm functions for float, double,
and long double (107-113 bit mantissa).  I get satisfactory elefunt
results everywhere, except that I have seen slightly better results for
exp().  I am using an implementation similar to the 4.3BSD with the
coefficients recommended by Cody & Waite.  My elefunt tests turn out
better (and faster) than BSD, but there may be a possibility of a 10%
improvement in RMS error.  On the other hand, my sinh() results are as
good as any, and those are sensitive to exp().  Can anyone reveal the last
secret of accuracy for exp()?
Tim
Return to Top
Subject: Re: Curve through tree points
From: jeppe@dfi.aau.dk (Jeppe Burchhardt)
Date: 11 Nov 1996 13:45:31 GMT
Jeppe Burchhardt (jeppe@dfi.aau.dk) wrote:
: Ricky Helgesson (rickyhelgesson@hotmail.com) wrote:
: : I will draw a curve through three points.
: : The three coordinates are "indata" and whish to get the equation for the
: : curve through these points so that I can plot that curve.
: : The equations are (of course):
: : y1 = a*(x1)^2 + b*x1 + c
: : y2 = a*(x2)^2 + b*x2 + c
: : y3 = a*(x3)^2 + b*x3 + c
: : The thing is that I have no program for calculating this, like matlab. I
: : know that in matlab, this is not too hard, so my question is:
: : Could some good person out there help me calculating this curve or give
: : me a hint where I can find matlab for Windows 3.11/95/NT?
: : Thanks in advance,
: : Ricky Helgesson
: : Com In Center AB
: : SWEDEN
: --
: ----
: Jeppe Burchhardt, Ph.D.
: Centre of Magnetic Resonance, Aarhus University Hospital, Denmark
--
Well, something went wrong with my newsreader (I guess, it wasn't actually
the newsreader but the operator :-() so I try again:
First of all, your three points must be distinct in the sense that
x1 != x2, x1 != x3, x2 != x3. If you define the quotients
q1 = y1/((x1-x2)(x1-x3))
q2 = y2/((x2-x1)(x2-x3))
q3 = y3/((x3-x1)(x3-x2))
you will find that the coefficients are given by
a = q1 + q2 + q3
b = -(x2+x3)q1-(x1+x3)q2-(x1+x2)q3
c = x2*x3*q1+x1*x3*q2+x1*x2*q3
Hope this helps
----
Jeppe Burchhardt, Ph.D.
Centre of Magnetic Resonance, Aarhus University Hospital, Denmark
Return to Top
Subject: Re: how do you generate a sine wave with simple adds and subtracts.
From: Stan Armstrong
Date: Mon, 11 Nov 1996 10:48:39 +0000
In article <560vdv$qns@mars.efn.org>, Stephen Dunbar
 writes
>
I had thought that that this was a special case of the algorithms that I
use for rotations, from x,y to X,Y which are stated as follows:
1)      X=x+T(y-Y)
2)      Y+y-T(x+X)
Where the X,Y terms on the right are at time t(n-1) and those on the
left at t(n) and T is the tan half angle.
If it were so setting T fixed at say .01 and making the new X and Y into
the old x and y for the next iteration would give X and Y as the sine
and cosine functions.  Unfortunately, as it turns out, those relations
are valuable where the vector length is given and are convergent for
rotations in the range +/-45 degrees.  (For larger angles you can rotate
the number of right angles to get the t into the x-axis quadrant.)
This may not help your current problem but could be a useful process for
later.  
-- 
Stan Armstrong
Return to Top
Subject: Re: Need help with integration
From: Stan Armstrong
Date: Mon, 11 Nov 1996 10:52:19 +0000
In article <55o568$6td@bach.convex.com>, Dave Dodson 
writes
>In article <327ED624.3093@worldnet.att.net>,
>Brian Gross   wrote:
>>What is the intergral of the following:
>>
>>                 X
>>           --------------- dx
>>              3
>>             X   -  1
>
>Use the method of partial fractions to write
>
>               x/(x^3-1) = a/(x-1) + b/(x^2+x+1)
>
>and integrate the pieces.
>
>Dave
Dave,  You b term needs to be b(2x+1)+c, the b part of which, like the a
fraction gives a log and the c part an arctan, according to the back of
my envelope.
Regards
-- 
Stan Armstrong
Return to Top
Subject: Re: leading dimension
From: Konrad Hinsen
Date: 11 Nov 1996 14:48:17 +0100
Zoltan Bus  writes:
> Somebody, please, define me what is the meaning of the expression
> 'leading dimension' which can be found in almost FORTRAN matrix
> subroutine. Excuse me for this silly question.
Not silly at all. The "leading dimension" is simply the length
of a column of an array, i.e. if you have
   real a(5, 7)
then the leading dimension is 5. Because of Fortran's memory layout
for arrays, subroutine have to know the leading dimension to
correctly access array elements. Knowing the second dimension is
not necessary unless the subroutine wants to do range checking
on indices.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                          | E-Mail: hinsen@ibs.ibs.fr
Laboratoire de Dynamique Moleculaire   | Tel.: +33-76.88.99.28
Institut de Biologie Structurale       | Fax:  +33-76.88.54.94
41, av. des Martyrs                    | Deutsch/Esperanto/English/
38027 Grenoble Cedex 1, France         | Nederlands/Francais
-------------------------------------------------------------------------------
Return to Top
Subject: Re: eigenvalues of ill-cond matrices
From: hwolkowi@orion.math.uwaterloo.ca (Henry Wolkowicz)
Date: Mon, 11 Nov 1996 16:00:53 GMT
In article <328706CB.48D1@physics.uwa.edu.au>,
Paul Abbott   wrote:
>alan koivunen wrote:
>
>>I'm trying to find eigenvalues and eigenvectors of some ill-conditioned
>>matrices. They are real, symmetric Toeplitz matrices, an example of 
>>which is coded below:
>> 
>>         do i=1,N
>>           do j=i,N
>>             A(i,j) = exp(-2.0*(j-i)*(j-i)/((N-1)*(N-1)))
>>             A(j,i) = A(i,j)
>>           enddo
>>         enddo
>> 
>>With N=13, the eigenvalue magnitudes range in value from about
>>10^(-16) to about 10. The eigenvalues should be real and positive,
>>but using MATLAB's eig() function and IMSL's DEVCSF routine produces
>>negative (small) eigenvalues - so things aren't quite right.
>> 
>>So, is there anyone out there who can suggest a way to get
>>eigenvalues/vectors for these troublesome matrices ? 
>
>Any of the Computer Algebra systems should be able to do this without
>any trouble.  Here is a solution using Mathematica 3.0:
>
>In[1]:= Toeplitz[n_]:=Table[Exp[-2(j-i)^2/(n-1)^2], {i, n}, {j, n}]
>
>In[2]:= Timing[Eigenvalues[N[Toeplitz[13], 35]]]
>
>Out[2]= {1.8 Second, 
>  {9.714289658177802948116981413446968039217238874`35, 
>   2.8087047226601903588888427367223952406986690664`35, 
>   0.4318113387746175560652633047322733083888887423`35, 
>   0.04211529596485480041989551847174304301749342335`35, 
>   0.0029201614706262355299008707132439452221134816464`35, 
>   0.00015247429079448480011310165012714444499654618885`35, 
>   6.150970244115566358292942430866977377286994027`35*^-6, 
>   1.9294232139272084071757708173097648232541969225`35*^-7, 
>   4.662923268040823404382745324585306129887057093`35*^-9, 
>   8.452676823084502937530316623177950895236479696`35*^-11, 
>   1.0890952782993752514438522541130430812858525505`35*^-12, 
>   8.941025867763620061655233278333984110022684667`35*^-15, 
>   3.531596847111418268275223468242875950227707268`35*^-17}}
>
>Cheers,
>	Paul 
>_________________________________________________________________ 
>Paul Abbott
>Department of Physics                       Phone: +61-9-380-2734 
>The University of Western Australia           Fax: +61-9-380-1014
>Nedlands WA  6907                         paul@physics.uwa.edu.au 
>AUSTRALIA                           http://www.pd.uwa.edu.au/Paul
>
>          God IS a weakly left-handed dice player
>_________________________________________________________________
Ill-conditioned matrices mean that you will have trouble solving systems
of equations, but it does not mean that you have trouble finding the
eigenvalues. One simple fix to avoid worrying about negative
eigenvalues, is to add
    20*I
to all the matrices. Then you will have a range for the eigenvalues from
   20 to 30
-- 
||Henry Wolkowicz                |Fax:   (519) 725-5441
||University of Waterloo         |Tel:   (519) 888-4567, 1+ext. 5589
||Dept of Comb and Opt           |email:  henry@orion.math.uwaterloo.ca
||Waterloo, Ont. CANADA N2L 3G1  |URL: http://orion.math.uwaterloo.ca/~hwolkowi
Return to Top
Subject: numerical integration of non linear equ
From: cooder
Date: Mon, 11 Nov 1996 18:22:02 +0100
I want to get some informations about the numerical integration
of non linear equations (with algorithms if possible);
for example RK, Adams, etc.
If you know some URL about it...
thanks...
 -- Cooder@easynet.fr --
Return to Top
Subject: Re: [Q] Need to plot in spherical coords.
From: dwells@nrao.edu (Don Wells)
Date: 11 Nov 1996 17:10:17 GMT
"CR" == Craig Ross  writes:
CR> I need to plot a surface in spherical coordinates ..
CR> ..  The data is spaced on a regular
CR> grid, from 0 <= phi <= 90 degrees and 0 <= theta <= 90 degrees..
CR> I currently have access to Matlab, Maple, Excel, and Axum but haven't
CR> found a simple solution..
There are an infinity of possible ways to map the sphere onto the
plane.  These are map projections with names like Mercator,
Aitoff-Hammer, Lambert, etc. Some map projections are good for
displaying the entire sphere, while others work well for limited areas
of the sphere. Some have the property that shapes aro preserved
locally, others can preserve areas while sacrificing shape
preservation. Navigators have often liked the Mercator because it
projects "rhumb lines" (lines of constant bearing) as straight lines,
whereas for other applications the 'gnonomic' projection is preferred
because it renders great circles as straight lines.
In the astronomy community this problem appears under the subject
"world coordinate systems" [WCS] in our interchange format, FITS
[Flexible Image Transport System]. A page containing an overview and
links to WCS documents is at
    http://fits.cv.nrao.edu/documents/wcs/wcs.html
The first link in that page, file wcs.all.ps (520 KB), is a paper
titled "Representations of celestial coordinates in FITS"; this is the
current draft of the proposed conventions for WCS in FITS. It is a
43-page definitive treatment of the mathematics and geometry of this
problem. NOTE: send this file to your fastest, most powerful
Postscript printer, and don't be dismayed if it takes many minutes to
come out -- there are 35 _complicated_ figures in this paper, which
make it a compute-intensive aerobic exercise for any Postscript printer!
Three sets of WCS source code are available in the directory
    http://fits.cv.nrao.edu/src/wcs/
     idl_wcs_code.tar.gz    04-Oct-94 16:36    28K  
     wcslib-2.4.tar.gz      10-Sep-96 10:46   285K  
     worldpos.tar.gz        13-Oct-94 17:04    16K  
My bet is that you could make a lot of progress by copying one or more
algorithms from worldpos.tar.gz, which supports 8 of the most popular
projective mappings; probably the 'AIT' (Hammer-Aitoff) case would be
a good start if you want to see the whole-sphere behavior of your
function. 'STG' (Stereographic, zenithal orthomorphic) or 'SIN'
(Orthographic) would be good choices if you want to see a 'polar cap'
region. The 'wcslib' and IDL packages support all 25 of the
projections discussed in the big WCS paper.
--
  Donald C. Wells         Associate Scientist         dwells@nrao.edu
                    http://fits.cv.nrao.edu/~dwells
  National Radio Astronomy Observatory                +1-804-296-0277
  520 Edgemont Road,   Charlottesville, Virginia       22903-2475 USA
Return to Top
Subject: Need your help
From: "Brad Marts"
Date: 11 Nov 1996 18:14:50 GMT
I was recently given the problem of finding the area under 1/(e^(x^2) + x)
from 0 to oo.  The answer only has to be accurate within .005.  This was
not so difficult for 1 to oo because 1/x^2 can be used.  However, I was
having trouble finding a function that has finite area from 0 to 1 and is
greater than the function in question.  If anyone has any ideas or
suggestions I would love to hear from you.  Either post to the newsgroup or
mail directly to me at 
bmarts@gchs.com
Thanks
Brad
-- 
Any unsolicited e.mail recieved at this address that is not relevant to the
post made can and will be prosecuted to the full extent of the law
Any e.mail recieved as a result of removing my e.mail address from this
post and adding it to a list can and will be prosecuted to the full extent
of the law
Return to Top
Subject: complex Newton's method
From: trscavo@forbin.syr.edu (Tom Scavo)
Date: 11 Nov 1996 17:28:03 GMT
Hi,
Anyone who has studied Newton's method on the real line knows its 
geometric interpretation, that is, a sequence of tangent lines whose 
zeros converge to a root of a function.  What is the corresponding 
geometric interpretation of Newton's method in the complex plane?
Explanations or pointers into the literature will be most appreciated.
If you post, please cc trscavo@syr.edu as well.
Cheers,
-- 
Tom Scavo
mailto:trscavo@syr.edu
http://web.syr.edu/~trscavo/
Return to Top
Subject: Quadratic constrained LS minimisation
From: Moritz Harteneck
Date: Mon, 11 Nov 1996 15:51:57 +0000
Hi,
I want to minimise a least squares problem which is quadratically 
constrained, i.e.
minimise( vAv ) under (vBv=p)
where v is a n-dimensional vector of the free parameters, A describes the 
performance criterion and B and p the quadratic constraints.
Does anybody know any good pointers?
Thanks
	Moritz
===========================================================================
Moritz Harteneck
SPD / EEE	                    email: moritz@spd.eee.strath.ac.uk
University of Strathclyde	    http://www.spd.eee.strath.ac.uk/~moritz
Glasgow G1 1XW, Scotland, U.K.
Tel: ++44-(0)141-552 4400 ext 2808  Fax: ++44-(0)141-552 2487
Return to Top
Subject: Re: Accurate exp()
From: "Dann Corbit"
Date: 11 Nov 1996 19:07:20 GMT
You might like to look at the implementation in MPFUN
for high accuracy versions.  MPFUN is a FORTRAN
library of multiple precision routines, and is available on
the net.  The implementation of exp() is quadratically
convergent.
n8tm@aol.com wrote in article <19961111140000.JAA28726@ladder01.news.aol.com>...
> I have nearly completed revision of my libm functions for float, double,
> and long double (107-113 bit mantissa).  I get satisfactory elefunt
> results everywhere, except that I have seen slightly better results for
> exp().  I am using an implementation similar to the 4.3BSD with the
> coefficients recommended by Cody & Waite.  My elefunt tests turn out
> better (and faster) than BSD, but there may be a possibility of a 10%
> improvement in RMS error.  On the other hand, my sinh() results are as
> good as any, and those are sensitive to exp().  Can anyone reveal the last
> secret of accuracy for exp()?
> Tim
> 
Return to Top
Subject: Re: Directed rounding on the Pentium
From: "Stephen W. Hiemstra"
Date: 11 Nov 1996 18:26:15 GMT
Tim,
Good point.  I will check it out.
Stephen
n8tm@aol.com wrote in article
<19961111134500.IAA28551@ladder01.news.aol.com>...
> I expect that it would be necessary to display the entire result (in
> hexadecimal, for instance) in order to check the effects of directed
> rounding.  The changes in the low order bit will not display neatly in
> decimal, although they will show up if enough places are printed (17 for
> the usual double).
> Tim
> 
Return to Top
Subject: running FDTD fortran codes on SGI
From: Tim Black
Date: Mon, 11 Nov 1996 14:17:52 +0800
I am getting ready to buy a box to develop/run some Finite Difference
Time Domain fortran codes for modeling antenna near-field
radiation/characteristics. I'm looking at SGI's O2 entry level
workstations(because of the power/price) and was wondering if anyone had
any advice on choosing a workstation for such computationally intensive
tasks. I have no experience with SGI, and would like to find some
benchmarks comparing Pentium Pros, SGI O2's SGI Indigo's, etc. Any
clues?
Thanks,
Tim Black
trblac00@pop.uky.edu
University of Kentucky
EE Dept.
Return to Top
Subject: Re: Need your help
From: "Brad Marts"
Date: 11 Nov 1996 18:42:12 GMT
Ummmm... Never mind.  I figured out my mistake.  I probably sounded like a
complete fool but such is life.
Brad
Return to Top
Subject: Curl and divergence estimation from irregularly spaced data
From: Olaf Amm
Date: Mon, 11 Nov 1996 16:37:10 -0800
Hi,
I am looking for an algorithm that is able to estimate the curl and
divergence of a vector field from measurements of that field, given on
irregularly spaced points.
It would be best if that algorithm could do it without gridding the data
prior to the estimation; however, i would also be thankful if somebody
can recommend a good gridding routine for my purpose,
Thanks in advance,
			Olaf
======================================================================
  Olaf Amm                            Tel: + 358 - 0 - 1929 4630
  Finnish Meteorological Institute    Fax: + 358 - 0 - 1929 4603
  Department of Geophysics
  P.O.BOX 503                         E-mail: olaf.amm@fmi.fi
  FIN-00101 Helsinki
  Finland
======================================================================
Return to Top
Subject: Re: Directed rounding on the Pentium
From: "Jeffery J. Leader"
Date: Mon, 11 Nov 1996 10:53:26 -0800
Chris Hecker wrote:
> "Jeffery J. Leader"  writes:
> >Kahan gave an interesting talk on this type of thing at the Gragg
> >Conference in Monterey last weekend
[...]
> That sounds great.  Is the paper on the web anywhere?
Not to the best of my knowledge.  I think there's supposed to be a
proceedings published--or maybe a special issue of LAA or something. 
Carlos Borges (mailto:borges@math.nps.navy.mil) has the relevant info.
on that.  Kahan may have it somewhere on the web--the title was
something like "Things I can do that you can't" referring to the
relatively few compilers that support this stuff.  (Elseweher in this
thread the original poster points out that it's all available in
assembly, of course.  Just not in something useful, with a few
exceptions here and there--some Borland stuff, for example.)
-- 
The tools of the mathematician's trade are pencil and paper: in 
consequence, no mathematician ever carries either, and they always 
have to borrow a pen and scribble on a napkin.
 -Ian Stewart
Return to Top
Subject: Re: Directed rounding on the Pentium
From: "Jeffery J. Leader"
Date: Mon, 11 Nov 1996 10:53:26 -0800
Chris Hecker wrote:
> "Jeffery J. Leader"  writes:
> >Kahan gave an interesting talk on this type of thing at the Gragg
> >Conference in Monterey last weekend
[...]
> That sounds great.  Is the paper on the web anywhere?
Not to the best of my knowledge.  I think there's supposed to be a
proceedings published--or maybe a special issue of LAA or something. 
Carlos Borges (mailto:borges@math.nps.navy.mil) has the relevant info.
on that.  Kahan may have it somewhere on the web--the title was
something like "Things I can do that you can't" referring to the
relatively few compilers that support this stuff.  (Elseweher in this
thread the original poster points out that it's all available in
assembly, of course.  Just not in something useful, with a few
exceptions here and there--some Borland stuff, for example.)
-- 
The tools of the mathematician's trade are pencil and paper: in 
consequence, no mathematician ever carries either, and they always 
have to borrow a pen and scribble on a napkin.
 -Ian Stewart
Return to Top
Subject: Re: [Q] Using pseudoinverse in Bayes discriminant function?
From: Greg Heath
Date: Mon, 11 Nov 1996 14:53:56 -0500
On Fri, 8 Nov 1996, Dukki Chung wrote:
> Hi. Reently, I had to use Bayes classifier for a pattern classification
> problem. The Bayes discriminant function is:
> 	di(x) = - [ ln|Ci| + (x-mu)^t Ci^-1 (x-mu)]
        di(x) = - [ ln|Ci| + (x-mui)^t Ci^+ (x-mui) -2 lnPi]
> The problem was, the covariance matrix Ci was near singular, so the
> inverse could not be calculated. So, I used pseudoinverse instead of real
> inverse. What I'm wondering is whether this is a valid, justifiable 
> mathematical or statistical approach.
Yes. I've always used the pseudoinverse. The ill-conditioning of the 
covariance matrix results in near zero eigenvalues corresponding to 
directions in space for which the distribution has nearly a constant 
value(i.e., nearly a zero variance).
> I would be appreciated for any comments, suggestions, references, or any 
> pointers.
Check the eigendirections associated with the near-zero eigenvalues.   
Classes with near constant values in those directions might be able to be 
classified quite easily based on that fact alone.
Hope this helps.
Gregory E. Heath     heath@ll.mit.edu      The views expressed here are
M.I.T. Lincoln Lab   (617) 981-2815        not necessarily shared by 
Lexington, MA        (617) 981-0908(FAX)   M.I.T./LL or its sponsors
02173-9185, USA
Return to Top

Downloaded by WWW Programs
Byron Palmer