Newsgroup sci.math.num-analysis 28636

Directory

Hi everybody!
I am looking for a PC routine for solution differetial equation by
orthogonal collocation method on finite elements.
 Thanks in advance!
e-mail:Darka.Ban@ijs.si

Hi to all number crunchers,  
I'm glad to announce the first test release of my pentium optimized  
BLAS level 1 library. The most important functions of the BLAS level 1 
(xAXPY, xCOPY, xDOT, xNRM2, in both single and double precision)  
have been completely rewritten in GNU assembler for achieving full 
performance on intel pentium processors. The library is available in 
both binary and source form for DOS/DJGPP and UNIX/LINUX at 
        http://cip.physik.uni-wuerzburg.de/~mlkessle/blas1.html 
See there for further informations. Suggestions, bug reports,  
comments etc. are very welcome, email adress see below. 
Thanks,  
        Manuel 
-- 
Manuel Kessler
Graduate Student at the University of Wuerzburg, Germany, Physics Department
SNAIL: Zeppelinstrasse 5, D-97074 Wuerzburg, Germany
EMAIL: mlkessle@cip.physik.uni-wuerzburg.de
WWW: http://cip.physik.uni-wuerzburg.de/~mlkessle

Anyone know where to find C source code to solve for the eigenvalues and 
eigenvectors of a complex self-adjoint matrix, or a C version
of EISPACK or LAPACK?

Hello to you all,
I don't know if this is the right newsgroup to post in,
but I have a question about odds on a casino table.
The situation is the following:
The Table has 3 rows, which each have 15 different numbers in it, and
there is the 0.
So in total there are 46 numbers, so the odds for getting a specific
number is 1/46, the odd for getting a specific row is 15/46.
Now when row 1 and 3 are thrown (that means that only one number in
that row had to be thrown) 15 and 20 times, and row 2 hasn't been
thrown yet.
Now are the odds now better to bet on row 2, cause the other two have
been throw more, and row 2 hasn't been thrown at all.
The throw device is neutral, so has no number one choice.
Any help would be appreciated, I think the odd for row 2 are better
since the other two have been thrown more, and because the device is
neutral ALL numbers are thrown as much times in the end.
Thanks
Arno Zwegers
amzweg@cistron.nl
http://www.cistron.nl/~amzweg
Arno Zwegers
amzweg@cistron.nl
http://www.cistron.nl/~amzweg

In article <56bc1t$58r@hcunews.hiroshima-cu.ac.jp>,
Hideo Hirose   wrote:
>In Japan, many researchers pronounce LaTeX as "latef." Is it correct? How do you 
>pronounce TeX and LaTeX actually, especially in the united states?
Quite close.
 Knuth's intention was to have it pronounced with the "X" sounding 
like the throat-clearing "ch" in "Loch Ness", from the mouth of a 
Scottish speaker, or as in Spanish "Juan". (X played this role in 
the classical Greek alphabet, and it's the same in Cyrillic as 
used by Russians and many others.)
The most frequently heard US pronunciation, by those who see the
Latin letter X (as I would see it if I didn't know the 
accompanying trivium), is of course "Leytecks".
Hope it helps, ZVK (Slavek).

In article <56bc1t$58r@hcunews.hiroshima-cu.ac.jp> Hideo
Hirose, hirose@cs.hiroshima-cu.ac.jp writes:
>In Japan, many researchers pronounce LaTeX as "latef." Is it correct? How do you 
>pronounce TeX and LaTeX actually, especially in the united states?
>
>
I pronounce it so it rhymes with  'luxury yacht'
Mike Yukish
Applied Research Lab
may106@psu.edu
http://elvis.arl.psu.edu/~may106/

In article ,
Moritz Harteneck   wrote:
>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
>
>===========================================================================
>
Here are some comments from an old matlab program I wrote:
% constr_opt_eig.m                         Mike Elgersma     October 19, 1995 
%--------------------------------------------------------------------------
% Quadratically constrained optimization, with Lagrange multiplier lam.
%
% Minimize    (x-x0)'*H*(x-x0)     subject to     x'*Q*x = eps'*eps
%
%    where  '  is conjagate transpose  and  eps is a scalar
%-----------------------------------------------------------------------------
% Blaise Morton pointed out that a shrinking ellipsoid is tangent to a fixed 
% ellipsoid with outward pointing normal vector for only one value of the 
% shrink factor.  This results in the Lagrange multiplier having one real 
% solution with one sign, while all other real solutions have the other sign.
%
%       . .
%      .   .
%      .   . . . . . . 
%      .   ..    x0    .   As ellipse expands, it is tangent to a fixed ellipse
%      . 0 . . . . . .     for four different values of the expansion factor,
%      .   . expanding     but only the first time the normal points outward.
%      .   .    
%       . .
%
%      fixed
%-----------------------------------------------------------------------------
% Set  f(x,lam) = (x-x0)'*H*(x-x0) + lam*( x'*Q*x - eps'*eps )
%
% If H and Q are real, then x and lam will be real at minima of f.
% If H or Q are not real, f may not be real, so x and lam may be complex when
% df/fx and df/dlam are zero.
%-------------
% f is minimum when  0 = .5*df/dx = H*(x-x0) + lam*Q*x
%                    0 = df/dlam = x'*Q*x - eps'*eps
%-------------
%   x = inv(H + lam*Q) * H*x0
%
% eps'*eps = x0'*H'*inv(H' + lam'*Q') * Q * inv(H + lam*Q)*H*x0
%          = x0'*H'*inv(H' + lam'*Q') * Q * (    x            )
%          = x0'*H'*(               y                         )
%-------------
% Let G(lam) = (H + lam*Q)
%
% then:
%
% eps'*eps = (H*x0)'*y
%     H*x0 = G*x
%      Q*x = G'*y       Q may not be invertible?
%--------------
%
%  (H*x0 * (H*x0)')*y = H*x0*((H*x0)'*y) = eps'*eps*H*x0 = eps'*eps*G*x
%
%  so      [eps*G  H*x0*(H*x0)']*[eps'*x] = 0
%                                [  -y  ]
%
% Combining with   Q*z = G'*y  gives:
%
%  [    Q      eps'*G(lam)']*[eps'*x] = [0]   Almost an eigenproblem for lam,
%  [eps*G(lam) H*x0*(H*x0)'] [ -y   ]   [0]   but G'(lam) has lam'.
%
% If the data matrices are real, then replace conjagate-transpose with 
% transpose so a true eigenproblem.
% Hermitian, so full set of eigenvectors, so can deflate in power iteration.
% But could have repeated eigenvalues which slow down convergence.
% If lam(1) = lam(2)+epsilon with distinct eigenvectors, then slow convergence.
%
% If the data matrices are comlex:
%    1) OK if some eigenvalues are real.
%    2) Does the optimum occur at some real eigenvalue?
%----------------------------------------------
% If x0 = 0
% 
% then
% 
% f is minimum when 0 = .5*df/dx = H*(x-x0) + lam*Q*x
%                                = (H + lam*Q)*x      a generalized eig problem
% but is   x'*Q*x = eps'*eps  for the solutions of the eig problem?
%------------------------------------------------------------------------------
% CONVERGENCE
%
% If only one real positive eig (and all complex eigs have neg real part)
% and we have a bound:
%   0 < lam_pos < lam_max
% then can guarantee convergence of the real inverse-power iteration by setting
% the shift equal to lam_max.
% If we use Dale's bi-section algorithm to get an upper bound for lam_pos
% that is less than a factor of two too large, then the convergence factor is:
% error(i+1)/error(i) <  |lam_pos - lam_max|/|lam_neg - lam_max|
%                     <  |lam_pos - lam_max|/|lam_max|
%                     <  1/2
% so convergence will be at the same rate as for Dale's bisection algorithm.
% 
% We can also speed up convergence by choosing a good initial value for
% the eigenvector in the inverse iteration.  
%  Set initial  x = x0/sqrt(x0'*Q*x0)      so x'*Q*x = 1 
% 
%------------------------------------------------------------------------------

Here is the problem: 
        Let K = {0, 1, ..., k - 1} (k > 1) be a set of k logic values.
Let Union and Intersection be two operators defined on K. Union is
defined as the bitwise OR operation between two elements represented in
binary numbers (having each log(k) bits, the base of the log is 2). Thus
for instance, for k = 8 we have, 1 Union 2 = 001 Union 010 = 011 = 3.
Intersection is defined as the bitwise AND operation between two
elements represented in binary numbers. Thus for instance, for k = 8 we
have, 5 Intersection 6 = 101 Intersection 110 = 100 = 4.
        Let k be a power of 2 (i.e. k = 2^r, with r > 0). Unary central
relations are the non-empty and proper subsets of K. A unary central
relation R is closed under Union and Intersection if x Union y and
x Intersection y are in R whenever x and y are in R. In other words:
(x in R and y in R) implies (x Union y is in R and x Inter y is in R).
        Now the problem statement: For k = 2^r, how many unary central
relations are closed under Union and Intersection ?
                                 ---
        I was not able to find a closed-form formula or even a recursive
formula.
        Using a computer program I find that, for k = 8 and k = 16,
there are 72 and 730 such unary central relations, respectively.
        There are eleven (11) unary central relations for k = 4:
        {0} {1} {2} {3} {01} {02} {03} {1,3} {2,3} {0,1,3} and {0,2,3}.
                        !!!GOOD LUCK!!!
With best regards.
A.N.

>==========A Montvay {TRACS}, 11/4/96==========
>Hi all,
>
>I want to use a Runge-Kutta (or similar) method on an IBVP on the
>three-dimensional wave equation
>
>d2p/dx2+d2p/dy2+d2p/dz2-d2p/dt2=0
it's second-order in time, so you need another variable for the
first derivative in time: use
    dp
q = --
    dq
Use Runge-Kutta on the vector:
   [    ]     [                             ]
 d [ p  ]     [ q                           ]
-- [    ]  =  [                             ]
dt [ q  ]     [ d2p/dx2 + d2p/dy2 + d2p/dz2 ]
   [    ]     [                             ]
That makes it a first-order ordinary differential equation, but it's
for a vector. The vector is not just two elements, it's for as many
pairs of elements as there are points in the space you're modelling.
So, for each point in space, you maintain two variables: p and q
calculate the derivatives d2p/dx2  etc. by a linear method: ie it's
just the second difference of p, divided by the x- y- or z- distance^2.
Simon

Dear all,
we're facing a problem with a simple Gaussian function.
Does anybody know the closed form of:
(LATEK) \frac{d^m}{dx^m}(e^{-x^2})
(Visual)
.             2
.        m  -x
.       d  e
.      --------
.           m
.        d x
for any m in N?
Any help will be greatly appreciated.
Kind regards,
Marco magix@dibe.unige.it
Giorgio nemo@dibe.unige.it

In article <3286C503.83F@pop.uky.edu>, Tim Black   wrote:
>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?
The limitation on FDTD is basically floating point speed.  The exception
is for very large problems where you may use virtual memory, if you go the
easy route, or asynchronous I/O, if you are clever.  In those cases bus 
speed, disk speed, and I/O bandwidth become issues.  
In any case, there is a large amount of information (including floating 
point performance) about many machines, from Gateway PC's to massively 
parallel machines at:
   http://performance.netlib.org/
Larry Seftor
seftor@mrj.com

I'm looking for an efficient routine to compute the inverse of the
error function (or the inverse of the cumulative Normal distribution).
This routine is available in Matlab (ERFINV), but I need a C (or
fortran) implementation.
Any pointer appreciated.
-- 
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 

MathTools Ltd. is pleased to announce MAT, a Matlab Compatible 
C++ Matrix Class Library, designed for development of advanced 
scientific high-level C++ code. Evalution version of the MAT 
can be downloaded from our home page, http://www.mathtools.com.
The library includes Complex math, Binary and unary operators, 
Powerful indexing capabilites, Signal processing, File I/O, Linear
algebra, String operations and Graphics. Over 300 mathematical 
functions are included in MAT.
MAT supports matrices of doubles, floats, ints and chars
mixed in the program. Images can be stored in matrices of chars,
using 1/8 memory storage. On many applications, where 8 digits of
precision are sufficient, float-precision matrices can save half
the memory usage. Memory allocation and  de-allocation is managed
automatically.
2D and 3D Graphics are supported in standalone applications using
the freely available graphics package, Gnuplot.
Fully functional, time limited (30 days) evaluation version of
MAT can be downloaded freely from our home page,
http://www.mathtools.com.
An introductory price for MAT is 250$ for a single user
license, with discounts on quantity and for academic institutes.
Prices are valid until Jan 1, 1997.
Sincerely,
The MathTools.
 =====================================================================
 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 <3287C6B6.224@sumppu.fmi.fi>, Olaf Amm  wrote:
> 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
Could do a radial basis function fit of each component, then differentiate
the fits.  
For scalar data u_i at n points x_i (in m-dimensions) fit
u(x)=\sum_{j=1}^n a_jf(|x-x_j|)  +b.x +c
where f(r)=r^2\log r or f(r)=r^3.
There are n+m+1 unknowns: a_j, vector b, c:
n equations come from u(x_i)=u_i
m equations from \sum_{j=1}^na_jx_j=0
and 1 from   \sum_{j=1}^n a_j=0
This is a symmetric system of equations.  For economy chose about a dozen
data points in the neighbourhood of the place you want to estimate the
derivatives.
                                                   Tony
---------------------------------------------------------------------
Professor A.J. Roberts     
Dept of Mathematics & Computing     E-mail: aroberts@usq.edu.au
University of Southern Queensland   Phone:  (076) 312943
Toowoomba, Queensland 4350          Fax:    (076) 312721
AUSTRALIA                           WWW: http://www.sci.usq.edu.au
                                    /pub/MC/staff/robertsa/home.html
---------------------------------------------------------------------

I have an array of data, about 27 by 36, each row and column in the array
is not necessarily evenly spaced (i.e. the axes could be 0, 10, 20, 35, 50, 60, ...
and 0, 2.5, 5, 10, 15, 20, 25,...). (This data is sampled from the output of
a street light on various angles.)
I need to be able to get an average for a give 'square' defined
points (upper left and lower right) which are arbitrary. The current method
is to use Simpson's rule on the middle row of data that crosses through the
square, but this is not accurate enough, I want to use a method that takes advantage
of all the data I have in 2 dimensions. I don't really know even what sort of
algorithm to look for, this is not my area of expertise.
Please help, I'm really stuck! :-)
Thanks,
Hamish

Does anyone know of any good and simple arbitrary tetrahedral mesh
generators.  My vision for the best mesher for my purposes would be
one that produces a high quality mesh (matches geometry very closely)
and can take the geometry (surface definitions) from a CAD program.
Should also be able to quickly refine the mesh in areas of interest.
The mesher should be able to write standard format output to feed into
my analysis program . My goal is to be able to mesh in an hour (define
geometry and mesh). Hopefully someone besides Santa Claus can help me!
Thanks in advance.
Rod Miller
rlmiller@afit.af.mil

Hi,
	Can someone refer me to a book or source where I can find how to 
solve PDEs using the method of orthogonal collocation.  I am attempting
to solve an Adsorption equation that is one dimensional in space and
time, second order in space and nonliear. Does anyone know of a package
that will solve such a problem?  I have IMSL package, but I could not
find a suitable routine for solving this problem.
Thanks
Venu

Mike Johnston (johnston@dimensional.com) wrote:
: Hi,
: Is there any free (or almost free) software available (preferably
: source code) to calculate the thermodynamic properties of steam and
: water?
: I have found several places that charge $200 and up. I just want this
: for entertainment and can't justify that kind of $$.
Years ago we published a paper with simple continuous polynomial definitions
for major thermodynamic properties of saturated steam and water in
SI and US units as a function of pressure. Temperature, enthalpy, density,
viscosity, and thermal conductivity are correlated.
Residuals are plotted over the range of validity. The reference is:
"Saturated-Steam-Property Functional Correlations for Fully Implicit
Thermal Reservoir Simulation" by WS Tortike and SM Farouq Ali,
SPE Reservoir Engineering, November 1989 pp 471-474.
These are easy to implement in a spreadsheet or otherwise program.
--
Dr. W.S. (Simon) Tortike, CPEng.
Senior Reservoir Engineer
    & Manager, Engineering Business Development
Mincom Petroleum Technology | WWW: http://www.mincom.com
Mincom Pty Ltd              | Tel: +61 7 3364-9821  (UTC/GMT +10)
138 Juliette Street         | Fax: +61 7 3364-9802
Greenslopes Qld 4120        | e-mail: tortike@mincom.com      | | | | | |
PO-Mail:  PO Box 72, Stones Corner, Qld 4120, Australia       |M|I|N|C|O|M

Mike Johnston (johnston@dimensional.com) wrote:
: Hi,
: Is there any free (or almost free) software available (preferably
: source code) to calculate the thermodynamic properties of steam and
: water?
: I have found several places that charge $200 and up. I just want this
: for entertainment and can't justify that kind of $$.
Years ago we published a paper with simple continuous polynomial definitions
for major thermodynamic properties of saturated steam and water in
SI and US units as a function of pressure. Temperature, enthalpy, density,
viscosity, and thermal conductivity are correlated.
Residuals are plotted over the range of validity. The reference is:
"Saturated-Steam-Property Functional Correlations for Fully Implicit
Thermal Reservoir Simulation" by WS Tortike and SM Farouq Ali,
SPE Reservoir Engineering, November 1989 pp 471-474.
These are easy to implement in a spreadsheet or otherwise program.
--
Dr. W.S. (Simon) Tortike, CPEng.
Senior Reservoir Engineer
    & Manager, Engineering Business Development
Mincom Petroleum Technology | WWW: http://www.mincom.com
Mincom Pty Ltd              | Tel: +61 7 3364-9821  (UTC/GMT +10)
138 Juliette Street         | Fax: +61 7 3364-9802
Greenslopes Qld 4120        | e-mail: tortike@mincom.com      | | | | | |
PO-Mail:  PO Box 72, Stones Corner, Qld 4120, Australia       |M|I|N|C|O|M

Konrad Hinsen  wrote:
>The question has already been adressed by Donald Knuth himself, who
>explained that the X in "TeX" is actually the Greek letter chi, and
>therefore should be pronounced like "chi" in Greek, "ch" in German,
>Irish, or Scottish, or like the corresponding sounds in Russian,
>Arabic, etc. That seems to be too difficult for most speakers
>of English, so what I hear in practice is either "tek" (in analogy
>to the English pronunciation of words like "technical") or "tex",
>especially "latex", i.e. totally ignoring the intended etymology.
You are right. If I recall correctly, he utilizes the example of
the word "loch" (as in "loch-ness") to give a close impression of
the correct pronounciation of the "ch" in "technology" and thus
in "TeX".
Tassos Machaerides
PS. Did you know that the word "art", in greek is "TEXNH"
    and the word "technology" is "TEXNOLOGIA".

On 13 Nov 1996, Michael Courtney wrote:
> Numerical recipes has algorithms for non-power-of-two numbers of points
> but not nonequispaced points.  Nonequispaced points is tough.  There are
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  in Numerical Recipes:
  13.6 Spectral Analysis of Unevenly Sampled Data,  page 569 
Zoltan Bus
trev@sas.seismology.hu

I have prepared a draft (chapter for proposed book)? paper concerning
performance issues of C vs FORTRAN, which anyone is welcome to look at and
perhaps criticize constructively.  About 45kB text.
Tim

Hi!
Suppose I have a Singular Values Decomp. of matrix A,
  A = U W V^T
now, I would like to find SVD of
 A + a b^T = U'W'V'^T
(where a, b are vector of an appropriate dimension.)
Is there any updating method to obtain U',W',V' by
using old U,W,V (or without starting the whole SVD again)??
Thank you.
s.s.
--
   .
<<<.><><><><><><><><><><><><><><><><><><>---=Sermsak Uatrongjit=--<><><>
  ___________________}{          [] Department of Physical Electronic []
<<>__________________|*O======{) () Tokyo Institute of Technology     ()
                     }{          []                                 . []
<><><>- boy@ss.titech.ac.jp ---<><><><><><><><><><><><><><><><><><><.>>>

Hello,
You may try the steam table software that comes along with the ASME steam 
tables book. It should be available in any library. If you are willing to
shell out some money, NIST provides a more accurate steam table software
for $200.00.
Good Luck,
Sreekanth

In article <328A0C2F.1489@ccgate.hac.com>, Dennis Braunreiter1  writes:
|> Anyone know where to find C source code to solve for the eigenvalues and 
|> eigenvectors of a complex self-adjoint matrix, or a C version
|> of EISPACK or LAPACK?
yes: CLAPACK in netlib
cheers, peter

The 1.02 version of EULER for Windows 95 is now available from 
am.ku-eichstaett.de by anonymous ftp. Installs from an EXE file, and can be 
de-installed, if desired. OS/2 and Unix/XWindows (esp. Linux) is available 
on the same site. The Windows 95 was reported to run on NT.
EULER started as a MatLab clone. It is now a program, which can handle 
real, complex and interval numbers and matrices, has a 2D/3D graphics, a 
builtin modern programming language (extension of MatLab's), an exact 
scalar product, and the Windows 95 version can call functions in an 
external DLL. The OS/2 and Windows versions interact nicely with the GUI, 
and have a notebook style interface.
The Unix version is free, teh OS/2 version free for educational use, and 
the Windows version cheap shareware.
R. Grothmann
grothm@ku-eichstaett.de

Alioune Ngom wrote:
> 
> Here is the problem:
> 
>         Let K = {0, 1, ..., k - 1} (k > 1) be a set of k logic values.
> Let Union and Intersection be two operators defined on K. Union is
> defined as the bitwise OR operation between two elements represented in
> binary numbers (having each log(k) bits, the base of the log is 2). 
> Intersection is defined as the bitwise AND operation between two
> elements represented in binary numbers. 
>       Let k be a power of 2 (i.e. k = 2^r, with r > 0). Unary central
> relations are the non-empty and proper subsets of K. A unary central
> relation R is closed under Union and Intersection if x Union y and
> x Intersection y are in R whenever x and y are in R. In other words:
> (x in R and y in R) implies (x Union y is in R and x Inter y is in R).
> 
>         Now the problem statement: For k = 2^r, how many unary central
> relations are closed under Union and Intersection ?
>                                  ---
> 
>    I was not able to find a closed-form formula or even a recursive
> formula.
> 
Nice problem, indeed. I've no solution yet, but a more general 
formulation:
"How many "sub-lattices", closed under + and *, 
 can be obtained from a Boolean lattice over r elements?"
For instance, the same problem arises if you consider 
a set of r elements and the lattice of all its susbsets, 
with standard set Union and Intersection.
Just a partial hint for a recursive formula: 
with r elements you can at least include all the
cases obtained for 1 <= s < r, thus retaining
only s out of r elements (=bit positions in your formulation), and
these elements can be chosen in comb(r s)=r!/((r-s)!s!) different ways.
The other (r-s) elements have fixed values.
E.g. In your example r=2. With s=1 (k=2), we have, in binary notation:
{0}{1}{0,1}, which can be extended to yield:
{0} -> {00}{01}{10} (these are trivial cases)
{1} -> {10}{11}{01}  
{0,1} -> {00,10},{01,11},{00,01},{10,11}
(it seems that principle of inclusion-exclusion should be used here
to avoid counting twice or more a same configuration) 
Let me know if you find a solution!
Paolo
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Paolo Ciaccia    
DEIS - CSITE-CNR         
University of Bologna - ITALY 
mailto:ciaccia@cs.unibo.it
http://www.cs.unibo.it/~ciaccia

Dear all!
This integral is so hard that systems
like Maple already need hours to
solve the following integral (with 4 digits precision)
for small values of n!
evalf(n*int(a*exp(-0.5*(a-100)^2/225)/(sqrt(2*Pi)*15)*
 (int(exp(-0.5*(x-100)^2/225)/(sqrt(2*Pi)*15),x=0..a))^(n-1),a=0..200));

"0" is an approximation for -oo
and "200" an approximation for +oo
These values need to be changed to lower or higher values
with increasing n. For the beginning try n=50 or so.
Desired values for n are 1e10 !!!
The intergral is - roughly speaking - a double integral
over two normal distributions, where the upper integration
border of the inner integration depends on the integration
variable of the outer integration. 
I need a fast implementation (in C) of a fast integration algorithm
for this integral.
Who can help me? Thanks in advance!
Bye, Daniel Gembris
-----------------------------------------------------------------------
This article was posted to Usenet via the Posting Service at Deja News:
http://www.dejanews.com/           [Search, Post, and Read Usenet News]

In article <5604gi$ga2@sanson.dit.upm.es>, a25467  wrote:
]-Thanks, I need some help to make a university work about anything based in the
]-linear resolution of matrix with an original !algoritmo! sorry but I am from
]-spain .
]-My work must be original (only a little) and i need some ideas 
]-I am in the fourth course of maths and My teacher only told us to do something
]-original for example :
]-- in pascal or ...
]-- new algoritms or a better solution to a theorem 
]-Of couse I need !bibliografy!.
]-Thanks and sorry for my english an for everything
]-
Here are some sources of numerical analysis code in Pascal...
Part 1.  Stuff I Know Something About...
----------------------------------------
[ Books ]
"Compact Numerical Methods for Computers" by J.C. Nash
    The source code that accompanies this book is available as freeware
    thru NetLib (http://www.netlib.org/).  This version was written in
    Turbo Pascal 5.0.  The book (2nd edition, copyright 1990) is also
    available (see http://www.ioppublishing.com).  The subtitle is "Linear
    Algebra and Function Minimization", which is a very accurate summation
    of the contents.  Of all the listings in Part 1, I've gotten the most
    mileage out of this one.
        [ ] FFT   [x] matrix methods   [x] nonlinear least squares
"Numerical Analysis, Third Edition", by Burden & Faires
    This is a standard textbook for undergraduate numerical analysis
    courses.  The code for all the algorithms presented in the book was at
    one time available in Pascal.  It's old stuff, written in Turbo Pascal
    3.0, and I don't know where you'd find it anymore...maybe write the
    publishers?  (Prindle, Weber, & Schmidt, ISBN 0-87150-857-5.  But I'm
    sure there's a later edition out by now, and I suspect that the code
    disk that accompanies it is not written in Pascal.)
        [x] FFT   [x] matrix methods   [x] nonlinear least squares
"Numerical Recipes in Pascal" by Press et. al.
    The code to accompany the text comes in two flavors:
        [1] a freeware version is available from several sites (SimTel =>
    ftp.coast.net, and Garbo => garbo.uwasa.fi come immediately to mind;
    look for nrpas113.zip)  This is very ugly code, not native Pascal, but
    a translation from the original FORTRAN, with no comments.  Personally
    I'd rather translate the stuff from "Numerical Recipes in C" than work
    with the Pascal freeware.
        [2] the native Pascal version of the first edition of Numerical
    Recipes is still available.  Check the Numerical Recipes homepage at
    http://nr.harvard.edu/numerical-recipes/ for more information.
        [x] FFT   [x] matrix methods   [x] nonlinear least squares
"Scientific Pascal" by Harley Flanders
    This is new (1996), and has received at least one very favorable review
    (http://users.aol.com/catambay/pascal-central/scipascal.html).  The code
    is written in Borland Pascal 7.0.  The book can be ordered through
    Birkhauser (http://www.birkhauser.com)
        [x] FFT   [x] matrix methods   [ ] nonlinear least squares
[ Software Packages ]
"Numerical Methods Toolbox"
    Another oldie, this one dates from about the release of Turbo Pascal
    4.0.  Borland put out a bunch of ToolBoxes -- Editor, Database,
    Graphics, etc. -- and this was one of them.  They've since dropped them
    from their product line.  I don't know where you'd get a copy...
        [x] FFT   [x] matrix methods   [x] nonlinear least squares
[ Downloads ]
SWAG => The SourceWare Archival Group
    This is a large collection of Pascal code fragments (mostly Borland) on
    almost every topic.  The software runs under DOS & Windows.  Freeware,
    available at SimTel, Garbo, and elsewhere.  At last count (swag9608)
    the math section had 116 entries, including FFTs, matrix algebra,
    complex numbers, etc.
        [x] FFT   [x] matrix methods   [ ] nonlinear least squares
Various collections...
    Several anonymous ftp sites offer collections of Pascal code and
    mathematical tidbits -- once again, SimTel and Garbo are good starting
    points.
Part 2.  Stuff I Know Nothing About...
--------------------------------------
[ Books ]
The following references are courtesy of Amara Graps via the
sci.math.num_analysis newsgroup.
    OK. Add to that:
    _Data Reduction and Error Analysis for the Physical Sciences_ by Philip
    R. Bevington and D. Keith Robinson, McGraw-Hill, 1992. This is an
    update to the original 1969 version by Bevington, which has code in
    Fortran. The Turbo Pascal 4.0 code for this book is on a floppy that's
    included in the book. (Note: the original poster was asking about
    multiple regression, and there's everything you ever wanted to know
    about multiple regression here.)
    _Scientific Programming with Macintosh Pascal_ (1986) by Richard E.
    Crandall and Marianne M. Colgrove. This book is somewhat dated, so the
    hardware-related items should be ignored. But there's some interesting
    statistics,  mathematical physics and signal processing examples here.
    The code would need to be typed in by hand. The author is the same
    Richard Crandall who wrote the more recent Telos books: _Projects in
    Scientific Computation_ and _Topics in Advanced Scientific Computation_
    (those Telos books use C and Mathematica for their examples).
    Amara
    --
    ***************************************************************
    Amara Graps                            amara@quake.stanford.edu
    Solar Oscillation Investigations       Stanford University
    http://quake.stanford.edu/~amara/amara.html
    ***************************************************************
During occasional forays into Net-Land I've also uncovered these...
"Numerical Methods for Mathematics, Science, and Engineering, 2e"
    John H. Mathews
    Prentice Hall, 1992 - ISBN 0-13-624990-6
    A table of contents can be found at
        http://www.prenhall.com/divisions/esm/esm_home.html
    (click on "Browse Titles by Course", then successively on "Math and
    Science", "Mathematics", "Math Courses for Engineers and Scientists",
    and "Numerical Analysis".)  Algorithm source is available in Pascal, C,
    FORTRAN, and MatLab.
I saw these next two while browsing through the Institute of Physics site
(http://www.ioppublishing.com).  Both include (optionally?) on-disk code
written in Turbo Pascal
"Computational Techniques in Physics" by P K MacKeown and D J Newman
    ISBN: 0 85274 548 6 (illustrated paperback),0 85274 429 3 (Software)
    Pages: 240
    Price: illustrated paperback/Software £19.00/US$37.00 £26.00/US$52.00
    Date: 1987
"Scientific Programmer's Toolkit, Turbo Pascal Edition" by M H Beilby,
    R D Harding, and M R Manning
    ISBN: 0 7503 0127 9
    Pages: 438
    Price: Boxed spiral-bound manual and software paperback £99.00/US$198.00
    Date: 1991
[ Software Packages ]
Matrix Math Toolkit v3.0 for Delphi 2
    Vector and matrix classes for Delphi 2.0, written by Peter N. Roth
    (peteroth@erols.com) of Engineering Objects International.  He also
    sells a finite element package called TProteus, which I believe is
    available in both Turbo Pascal and Delphi 1.0 versions.
        [ ] FFT   [x] matrix methods   [ ] nonlinear least squares
[ Downloads ]
The Turbo Pascal Math Library by Jean Debord
    available at http://ourworld.compuserve.com/homepages/JDebord/tpmath.htm
    What follows is a quote from the TP Math homepage...
    
    Welcome to TP Math, a mathematical library for Turbo Pascal (version
    5.0 or higher).  This library is distributed as Freeware.  You may use
    it anyway you like, as long as you don't charge money for it, remove
    the copyright notices, or hold anyone liable for its results.
    This library provides routines and programs for :
        o  Mathematical functions (including trigonometric, hyperbolic,
                    special, and some probability distributions)
        o  Vector and matrix operations (including dynamic array allocation)
        o  Function minimization (Marquardt, simplex)
        o  Statistical functions
        o  Random numbers (Marsaglia generator)
        o  Linear and nonlinear regression
    
        [ ] FFT   [x] matrix methods   [x] nonlinear least squares
----
The question of where to get numerical analysis code in Pascal pops up
occasionally in the various Pascal (and now Delphi) newsgroups.  If anyone
has some other sources I haven't mentioned here I'd appreciate hearing
about them.  Freeware, shareware, commercial...nothing is too obscure.
Mark Vaughan
m.a.vaughan@larc.nasa.gov

George Klinich wrote:
> 
> Hi,
> 
> I just read a paper where they solved the NLSE by a
> split-step spectral method.  The NLSE is
> 
>        dU         h^2
>  ih  -----  = - ------ grad^2 U + (v|U|^2 + w)U
>        dt          2m
> 
> and the method described solves the non-linear part in
> real-space and the Laplacian oprerator in Fourier space.
> Where can I find information on this method?
> 
Below is Mathematica code which outlines the basic method quite simply
For  i u_t = u_xx + 2 u |u|^2 + u v(x)
SplitStep[u_, v_, dt_, D2_] := NN[.5 dt, v][LL[D2][NN[.5 dt, v][u]]]
does one step split-step method, which is a composition of a linear
differential part (LL) and a nonlinear/potentail part (NN), whcih can
be defined as follows.
LL[D2_][x_] := InverseFourier[Fourier[x]*D2];
(Fourier returns a normalized DFT)
NN[ dt_, v_][x_] := x*Exp[-I*dt*(v - 2.Abs[x]^2)];
The differential vector can be generated by:
Dgenerate[n_, dt_, S_] := 
		Block[{s = -I*dt*S^2,t},
		t = Table[Exp[s*j^2],{j,1,n/2-1}];
		Join[{1},t,{Exp[s*(n/2)^2]},Reverse[t]]];
With these definitions, the following commands give you the
approximate solution for a stationary soliton for up to t = 1.  
n = 128; 
xmin = -8; xmax = 8;
S = 2 Pi/(xmax-xmin);
u = Table[Sech[xmin + (xmax-xmin)*(j/n)],{j,1,n}];
v = Table[0.,{j,1,n}];
dt = .01;nt = 100;
D2 = Dgenerate[n, dt, S];
unew = u;
Do[SplitStep[unew,v,dt,D2],{nt}];
The Mathematica code is adequate for small simulations and can be
useful as a prototype for developing a compiled code using your
favorite FFT.  
There are many refinements that can be made to improve the accuracy
and speed of the basic code, including skipping bitreversal in the
FFT, updating change to avoid accuracy problems, and higher temporal
order compostitions.

Hamish Hubbard  wrote in article
<56e202$rrf@cantuc.canterbury.ac.nz>...
> I have an array of data, about 27 by 36, each row and column in the array
> is not necessarily evenly spaced (i.e. the axes could be 0, 10, 20, 35,
50, 60, ...
> and 0, 2.5, 5, 10, 15, 20, 25,...). (This data is sampled from the output
of
> a street light on various angles.)
> 
> I need to be able to get an average for a give 'square' defined
> points (upper left and lower right) which are arbitrary. The current
method
> is to use Simpson's rule on the middle row of data that crosses through
the
> square, but this is not accurate enough, I want to use a method that
takes advantage
> of all the data I have in 2 dimensions. I don't really know even what
sort of
> algorithm to look for, this is not my area of expertise.
MATLAB has a function called "GRIDDATA" that I've used for this type of
problem before.  In their "HELP" file it claims to use an "Inverse Distance
method."  If you have MATLAB, then use the GRIDDATA command.  If not, then
maybe you could search for the "Inverse Distance method." 
Good luck,
John