Subject: Re: Question
From: jarausch@numa1.igpm.rwth-aachen.de (Helmut Jarausch)
Date: 18 Nov 1996 13:17:17 GMT
In article <328CC813.1569@highend.com>, Michael writes:
|> I'm solving an ODE system:
|>
|> d2_x/dt_2=f1(t,x,y,z,...)
|> d2_y/dt_2=f2(t,x,y,z,...)
|> d2_z/dt_2=f3(t,x,y,z,...)
|>
|> When I use a numerical method to calculate x(t), y(t), z(t), precision
|> is proportional to dt in n-th power. However,
|> I need to calculate y,z as functions of x. How can I estimate precision
|> in this case?
At time T you get \bar x = x(T) + x_eps
\bar y = y(t) + y_eps
where x_eps and y_eps are estimated errors (by the ODE scheme) which
are (in principle) proportional to dt in n-th power.
If you like to interpret y(T) as an Y(x) := y(x(T)) then
\bar y = y(T) + y_eps = Y(\bar x - x_eps) + y_eps
= Y(\bar x) - \tilde Y' * x_eps + y_eps
where \tilde Y' is an estimate for for Y' near \bar x
Now d/dt y(x(t)) = Y'(x)* d/dt x(t)
If you have rewritten your ODE-system as a first order system then the
ODE-solver gives you estimates for d/dt y(x(t)) and d/dt x(t)
and therefore an estimate for Y'(x) .
Keeping in mind that all these are estimates (not bounds) your error
sums up to - d/dt y(x(t)) / d/dt x(t) * x_eps + y_eps
which shows that you are in trouble if d/dt y(x(t)) / d/dt x(t) is
small.
Hope this helps a bit,
Helmut Jarausch.
Subject: SUMMARY: Computation of the inverse of the error function
From: Olivier Siohan
Date: Mon, 18 Nov 1996 14:10:45 GMT
I'd like to thank all the people who responded to my post.
QUESTION: where to get an efficient routine to compute the inverse of
the error function.
ANSWERS:
1. The package "dcflib" does the job. Cumulative and inverse
cumulative of many probability density function. C language. Available
on the net, search "dcflib"
2. Another package:
Source code available in the excellent new book:
Title : "A NUMERICAL LIBRARY IN C FOR
SCIENTISTS AND ENGINEERS"
(ISBN 0-8493-7376-X), 1995
Author : H. T. Lau
Publisher : CRC Press, Inc.
2000 Corporate Blvd., N.W.
Boca Raton, FL 33431
USA
telephone: (561)-994-0555
A diskette is included with all the source code.
3. Theoretical paper about rational approximations for the inverse of
the error function:
J.M. Blair, C.A. Edwards, and J.H. Johnson, "Rational Chebyshev
approximations for the inverse of the error function", Mathematics
of Computation, 30 (1976) 827-830 + microfiche appendix
4. Usefull link:
http://www.fred.net/mandalay
--
Olivier Siohan, Rm 2D-459
AT&T; Labs - Research, 600 Mountain Avenue, Murray Hill, NJ 07974
908-582-5877 - siohan@research.att.com - fax:908-582-4271
Subject: Re: g++ compiler v.2.7.2
From: Oleg Krivosheev
Date: 18 Nov 1996 09:16:17 -0600
Hi,
"M. Seto" writes:
>
> In a program that I have written, I have no problem compiling the
> following using an older version of g++:
>
> #include
> complex i(0, 1);
> int main()
> {
> complex M[2][2], r;
> complex a;
> ....
> etc.
> }
>
> Can someone kindly let me know how to declare complex variables with
> this new g++ compiler?
i believe, complex was converted into template class by
ANSI/ISO C++ commettee, so proper usage now is
#include
complex i(0, 1);
int main() {
complex M[2][2], r;
....
etc.
}
check C++ FAQ for public C++ draft location.
OK
Subject: HELP : I need a Spagetti Code Cleaner.
From: Matt Davies
Date: Tue, 19 Nov 1996 09:57:28 +1000
Hi Guys,
I need a spag code cleaner for some FORTRAN77 code. The code is very
structured but uses DO/CONTINUE instead of DO/ENDDO etc.
I would really like to hear from someone who has found such a beast...
Regards,
Matt D.
--
*=,._.,=*'~`~'*=,._.,=*'~~'*=,._.,=*'~`~'*=,._.,=*'~`~'*=,._.,=*'~`~'*=,.
: Matt Davies "TRAIN...DOOR...SPACE..." Ph: +61 7 3864 4309 PERSONAL :
: QUT School of Mathematics, Gardens Pnt, Fx: +61 7 3864 5252 OPINION :
: PO Box 2434, Brisbane Q 4001, Australia Em: m.davies@fsc.qut.edu.au :
*=,._.,=*'~`~'*=,._.,=*'~~'*=,._.,=*'~`~'*=,._.,=*'~`~'*=,._.,=*'~`~'*=,:
COPYRIGHT (C) 1996, M J Davies, All Rights Reserved. So there ... ;P
Subject: Re: Errata for Abramowitz/Stegun?
From: harper@kauri.vuw.ac.nz (John Harper)
Date: 18 Nov 1996 23:24:29 GMT
In article <32836B43.55A8@BBN.com>, Bill Marshall wrote:
>Is there an errata list for Abramowitz & Stegun?
My copy is the 9th Dover printing. Errata to that are:
p3 ln(sqrt(2*pi)) should end with 03297 not 03296
p6 def. of temperature: true when A+S wrote it, but now take that
degree symbol off the K.
p337 8.13 All references to E and K elliptic integrals are to E(k),
K(k), not E(m). K(m), in spite of p.590 saying otherwise.
p504 13.1.1 M(a,b,z) and U(a,b,z) are not independent solutions if a=0,
but then the equation is elementary for any value of b.
13.1.2 The series fails if b=0.
13.1.3 If a is an integer <= 0 and b = integer then the solution
indep. of M is not U. Use 13.1.12 - 13.1.19 instead.
13.1.5 U(a,b,z) is not many-valued if a-b+1 is an integer <=0.
p1038 Laplacian in spherical coordinates: for spherical read spheroidal.
I tell my students they really understand special functions if they find
a previously unknown error in A+S (or G+R).
John Harper Mathematics Dept. Victoria University Wellington New Zealand
Subject: ODE23/ODE45 Question
From: "W.J.Jenkins"
Date: Mon, 18 Nov 1996 21:34:56 -0500
My apologies if this is a FAQ, but I was wondering whether the "tol"
parameter in the ode23/ode45 function calls in MATLAB is an absolute
tolerance or relative to the actual value of the function being
modeled. For example, if the function is scaled such that at t=tf
it is equal to, say, 10^6, should one use a value of "tol" more like
"1" than the default of 10^-6 (for ode45)?
For example, the call is
[t,y]=ode45('func',0,tf,y0,tol);
Thanks for your suggestions and help
--
------------------------------------------------------------------------
W.J. Jenkins email: wjenkins@whoi.edu
MS #25
Woods Hole Oceanographic Institution fax: 508-457-2193
Woods Hole, MA 02543 http://kopernik.whoi.edu
Subject: Q: derivatives using MAPLE
From: kenny@scbspk.mit.edu (Sean Patrick Kenny)
Date: 19 Nov 1996 05:59:04 GMT
Given a function that contains various powers of time derivatives,
how do I calculate the partial derivative of the function wrt one of
the time derivatives?
A simple example of the problem:
F:=diff(x(t),t)^2+diff(x(t),t)*diff(y(t),t);
I need dF/xdot, and dF/ydot
where xdot = diff(x(t),t)
The answer I'm looking for is:
2*diff(x(t),t)+diff(y(t),t);
Thanks,
Sean
--
Subject: Re: [Q] Using pseudoinverse in Bayes discriminant function?
From: Greg Heath
Date: Mon, 18 Nov 1996 20:29:12 -0500
On Thu, 14 Nov 1996, Paul T. Karch wrote:
> wrote:
-----------------------SNIP----------------------------------------
> If possible, I would like to get references on Baysian classifiers,
> etc.; by mail or post . This is for self-study. Thanks in advance.
Devijver, Pierre A. and Kittler, Josef (1982), Pattern recognition:
a statistical approach, Englewood Cliffs, N.J.: Prentice/Hall International.
Duda, Richard O. and Hart, Peter E., (1973), Pattern Classification and scene
analysis, New York: Wiley.
Fukunaga, Keinosuke (1st ed. 1972, 2nd ed. 1990), Introduction to statistical
pattern recognition, New York(1st ed.), Boston(2nd ed.): Academic Press.
James, Mike (1985), Classification algorithms, New York: Wiley.
However, I can't remember seeing any reference that uses the pseudoinverse.
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
Subject: Re: Division by 2
From: kptben@aol.com
Date: 19 Nov 1996 07:15:55 GMT
John Hench wrote:
>
> 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
If the divisor is known at compile-time to be a power of 2, then most good
compilers will substitute the division instruction with a simple shift (a
straightforward type of optimization known as 'strength reduction.') If
the divisor is variable, then there's little the compiler can do. There is
an easy test for power-of-2'ness, however... for any integer x, if (x & (x
- 1)) == 0, then x is a power of 2. So if you know in advance that the
divisor is likely (but not guaranteed to be) a power of 2, you could write
an adaptive function like: [uint32 == unsigned long]
uint32 divide(uint32 a, uint32 b) {
if((b & (b - 1)) != 0) return a / b; // general case
// ____ otherwise, shift by appropriate amount to simulate divide
uint32 zcnt = __cntlzw(b); // count leading zeros of b
return a >> (31 - zcnt);
}
This only works if you have a "count leading zeros" instruction in
hardware... PowerPC has it at least. There's no real quick software way
to do it, at least not fast enough to be really useful. Also, signed
division requires a few more tweaks to handle the various cases.
Hope this helps,
--
Ben Weiss
Imaging Scientist
MetaTools Inc.
Subject: Just try this, it will work
From: rscott@fcc.com_ (R. Scott Perry)
Date: 19 Nov 96 03:21:23 GMT
In article <135.353194369934@news.nemonet.com>, sbolting@nemonet.com (Steve
Boltinghouse of 1009 Bird St., Hannibal, MO 63401) says...
>Take five minutes to read this and it WILL change your life.
Unfortunately, the change in Steve's life is that he may be spending the
next few years in jail.
Not only are chain letters illegal in the U.S. (2 year jail term for
first offence, see below), they CAN NOT WORK. There would have to be
over 1,000,000,000,000 people participating for them to work. Here's
PROOF:
---
The chain letters claim that you will make at least $10,000 (often they
claim $20,000 or $50,000), by people sending you $1 bills.
If the person who originally wrote the chain letter made $10,000 off of
it, then there are at least 10,000 people who have joined so far. That's
obvious; if someone made $10,000 from people sending him $1, then 10,000
people have joined so far.
In order for each of those 10,000 people to make $10,000, there would
need to be 100,000,000 people joining (10,000 people already in the
program, each causing another 10,000 to join; 10,000 times 10,000 is
100,000,000). That, too, is quite obvious.
So, if you saw this chain letter and joined it, according to the way the
chain letter claims it works, there would be 100,000,000 other people
joining. For all of you to make your $10,000, you would each need to
cause another 10,000 people to sign up. How many are we at now? It's
simple: 100,000,000 people times 10,000 new people each. That's
1,000,000,000,000 (1 trillion) people joining. That's about 200 times as
many people as there are on this planet.
Now do you understand why chain letters won't work?
---
Chain letters are illegal if the chain letter or money is transferred
via U.S. Mail, and can get you a $1,000 fine and a 2 YEAR jail sentence
for your first offense. The second offence can get you in jail for
5 years. Even if transmitted via the Internet, they violate Title 18
USC Section 1302, the Postal Lottery Statute. For verification, you
can currently go to http://www.usps.gov/websites/depart/inspect/
chainlet.htm to see that it is illegal, and go to http://www.usps.gov/
websites/depart/inspect/usc18/lottery.htm for the text of the law).
---
Feel free to copy this in response to any chain letters you see,
hopefully as enough people see this, the chain letters will stop.
-Scott
RSP
Subject: Re: eigenvalues of ill-cond matrices
From: spellucci@mathematik.th-darmstadt.de (Peter Spellucci)
Date: 19 Nov 1996 09:53:57 GMT
this was clearly a bad advice, since you may completely loose
precision in the small eigenvalues, which for some applications
are the most interesting and you reallly need them to some
relative precison. for computing
small eigenvalues to some precision Jacobi's method is often
better than QR (see the work by Veselic in Numerische Mathematik)
However, in the situation given by this posting there is no help
other than brute force in high precision computation.
peterIn article , hwolkowi@orion.math.uwaterloo.ca (Henry Wolkowicz) writes:
|> 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: Just try this, it will work
From: rscott@fcc.com_ (R. Scott Perry)
Date: 18 Nov 96 23:40:38 GMT
In article <135.353194369934@news.nemonet.com>, sbolting@nemonet.com (Steve
Boltinghouse of 1009 Bird St., Hannibal, MO 63401) says...
>Take five minutes to read this and it WILL change your life.
Not only are chain letters illegal in the U.S. (2 year jail term for
first offence, see below), they CAN NOT WORK. There would have to be
over 1,000,000,000,000 people participating for them to work. Here's
PROOF:
---
The chain letters claim that you will make at least $10,000 (often they
claim $20,000 or $50,000), by people sending you $1 bills.
If the person who originally wrote the chain letter made $10,000 off of
it, then there are at least 10,000 people who have joined so far. That's
obvious; if someone made $10,000 from people sending him $1, then 10,000
people have joined so far.
In order for each of those 10,000 people to make $10,000, there would
need to be 100,000,000 people joining (10,000 people already in the
program, each causing another 10,000 to join; 10,000 times 10,000 is
100,000,000). That, too, is quite obvious.
So, if you saw this chain letter and joined it, according to the way the
chain letter claims it works, there would be 100,000,000 other people
joining. For all of you to make your $10,000, you would each need to
cause another 10,000 people to sign up. How many are we at now? It's
simple: 100,000,000 people times 10,000 new people each. That's
1,000,000,000,000 (1 trillion) people joining. That's about 200 times as
many people as there are on this planet.
Now do you understand why chain letters won't work?
---
Chain letters are illegal if the chain letter or money is transferred
via U.S. Mail, and can get you a $1,000 fine and a 2 YEAR jail sentence
for your first offense. The second offence can get you in jail for
5 years. Even if transmitted via the Internet, they violate Title 18
USC Section 1302, the Postal Lottery Statute. For verification, you
can currently go to http://www.usps.gov/websites/depart/inspect/
chainlet.htm to see that it is illegal, and go to http://www.usps.gov/
websites/depart/inspect/usc18/lottery.htm for the text of the law).
---
Feel free to copy this in response to any chain letters you see,
hopefully as enough people see this, the chain letters will stop.
-Scott
RSP
