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
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
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
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
_________________________________________________________________
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
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
-------------------------------------------------
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.
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
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.
>
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
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
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
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
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
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
-------------------------------------------------------------------------------
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
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
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/
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
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
>
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
======================================================================
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
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
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