mm-1003
===
Subject: Re: small eigenvalues
> Suppose 0<=l1 <= l2 <= ... <=ln are the eigenvalues of a
real
> symmetric matrix. Let lr be the first non-zero eigenvalue,
i.e.
> 0=l1=...=l(r-1)<lr<=...<=ln.
> Obviously, numerically, the 0 eigenvalues are not exactly 0.
> There must be a mathematical argument or a heuristic for
deciding
> what r is. Forgive my ignorance, IÕm not a numerical
analysis
> person, but trying to learn these tricks as I run into them.
A practical test for zero might be the magnitude of L(I) is
less than
the maximum magnitude of all the eigenvalues times machine
epsilon in
working precision, which is close to the definition of a
singular
matrix adopted in LINPACK.
--
===
Subject: re: intriguing numbers
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision:
1.9 primary) id i12MKK228282;
i still dont have a clear enough picture of it. so youre
telling me
those are not true?
===
Subject: Re: high order polynomial interpolation
> Hi all,
> Given: a set of polynomials defined by a recurrence
relation, no
closed-form
> available.
> Required: the coefficients a[i] of a polynomial of this set
a[N]*x^N + ...
+
> a[0], N can be 40 or more
> Reason: The polynomials have a physical meaning and I need
to calculate
all
> the roots in a certain interval (which I do by calculating
the
eigenvalues
> of the companion matrix, hence why I need the coefficients,
and then
> verifying which of them lie in the interval of interest).
> For the moment, I do, as recommended in ÔAccuracy and
Stability of
Numerical
> AlgorithmsÕ:
> - Create N+1 Chebychev points
> - Do a Leja ordering on these points
> - Evaluate the polynomial in these points by its recurrence
> - Solve the corresponding Vandermonde system (using the
algorithm in
ASNA)
> - Build the companion matrix
> - Calculate the eigenvalues of the companion matrix
> This works reasonably well for degrees up to 40 (I can
compare to exact
> calculations in Maple using 200 digits fixed precision). For
higher
degrees
> the error becomes bigger.
> Question 1 : is there another way to do this ?
> Question 1a : would a least-squares approach with a lot
more data points
> given a more stable algorithm and hence better solution?
> Question 1b : is there a polynomial root finder available
that returns
all
> roots (in an interval) using only function evaluations and
does not rely
on
> the explicit knowledge of the coefficients?
> Question 2 : how can I make the previous algorithm more
stable?
> gert
I have no opinion as to whether what I suggest below is
better than what
others have suggested, but I think it of some general
interest. Denote
your basis polynomials as p_n(x). I assume that p_n is of
degree n, and
that it can be expressed as some combination of p_i, i < n,
and that the
multiplier of p_i is a polynomial of degree < n - i + 1. Let
your given
polynomial be sum_{i=0}^n b_i p_i((x).
Start at the top. b_n p_n can be expressed in terms of lower
order p_i.
To this combination add in the appropriate coefficients
b_{n-1} to the
multiplier of p_{n-1}. Now you can do the same to p_{n-1},
working your
way down till you get to p_0 which I hope is a constant. As
you are
working you way down, each of the p_i will be multiplied by a
polynomial
in the monomial basis.
An example for using this approach are codes for converting
between the
Chebyshev and Monomial bases at
http://mathalacarte.com/cb/mom.fcg/ya54.
You may find the description there of some value in deriving
what you
need in your case.
Fred
===
Subject: Re: high order polynomial interpolation
>Hi all,
>Given: a set of polynomials defined by a recurrence
relation, no
closed-form
>available.
>Required: the coefficients a[i] of a polynomial of this set
a[N]*x^N +
... +
>a[0], N can be 40 or more
>Reason: The polynomials have a physical meaning and I need
to calculate
all
>the roots in a certain interval (which I do by calculating
the
eigenvalues
>of the companion matrix, hence why I need the coefficients,
and then
>verifying which of them lie in the interval of interest).
>For the moment, I do, as recommended in ÔAccuracy and
Stability of
Numerical
>AlgorithmsÕ:
>- Create N+1 Chebychev points
>- Do a Leja ordering on these points
>- Evaluate the polynomial in these points by its recurrence
>- Solve the corresponding Vandermonde system (using the
algorithm in
ASNA)
>- Build the companion matrix
> - Calculate the eigenvalues of the companion matrix
>This works reasonably well for degrees up to 40 (I can
compare to exact
>calculations in Maple using 200 digits fixed precision).
For higher
degrees
>the error becomes bigger.
>Question 1 : is there another way to do this ?
>Question 1a : would a least-squares approach with a lot
more data
points
>given a more stable algorithm and hence better solution?
>Question 1b : is there a polynomial root finder available
that returns
all
>roots (in an interval) using only function evaluations and
does not rely
on
>the explicit knowledge of the coefficients?
>Question 2 : how can I make the previous algorithm more
stable?
>gert
> if you have a recurrence relation for the polynomial then
you have also
> a recurrence relation for its derivatives simply by
differentiating the
given one
> with respect to x. hence you could directly use methods
which compute
real
> zeros of a polynomial in an interval using the values of
derivatives
without ever
> using its coefficients, for example a damped Newtons or
Laguerres method
> combined with MaehlyÕs trick of implicit deßation or \
some
marching
scheme
> without any deßation (e.g. estimating the polynomial from
below or above
> by parabolas as used by Barth (published in Computing in
the early 70Õs)
> hth
> peter
Bravo! This is exactly what I was going to suggest, so it must
be right ;-)
--
Julian V. Noble
Professor Emeritus of Physics
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
God is not willing to do everything and thereby take away
our free will and that share of glory that rightfully belongs
to us. -- N. Machiavelli, The Prince.
===
Subject: Re: high order polynomial interpolation
>>if you have a recurrence relation for the polynomial then
you have also
>>a recurrence relation for its derivatives simply by
differentiating the
given one
>>with respect to x. hence you could directly use methods
which compute
real
>>zeros of a polynomial in an interval using the values of
derivatives
without ever
>>using its coefficients, for example a damped Newtons or
Laguerres method
>>combined with MaehlyÕs trick of implicit deßation or \
some
marching
scheme
>>without any deßation (e.g. estimating the polynomial from
below or above
>>by parabolas as used by Barth (published in Computing in
the early 70Õs)
>>hth
>>peter
> Bravo! This is exactly what I was going to suggest, so it
must
> be right ;-)
But if not all zeros are real, LaguerreÕs method is not
guaranteed
to converge. This may be a problem at high degrees. On the
other
hand, eigenvalue algorithms nowadays converge in practice
always
(although I think rigorous convergence proofs are still
lacking
in the nonhermitian case).
Arnold Neumaier
===
Subject: Re: Dahlquist/Bjorck Numerical textbook online
Prof Bill Smith of University of Guelph once told me the
BjorckÕs book
was best for least squares solutions, especially wrt SVD. Now
here we
have newest with his co-author Dahlquist.
> It is just so good of these people to put this out on the
net!!!<<<
all the best from Canada
and as we say in Swedish
tuck saw mewcket
Paul
> At http://www.mai.liu.se/~akbjo/NMbook.html , volumes 1 and
2 of the
> 3-volume series Numerical Mathematics and Scientific
Computation, by
> available online in PostScript format.
===
Subject: Re: Dahlquist/Bjorck Numerical textbook online
> At http://www.mai.liu.se/~akbjo/NMbook.html , volumes 1 and
2 of the
> 3-volume series Numerical Mathematics and Scientific
Computation, by
> available online in PostScript format.
was economically feasible to publish more books online.
-John
===
Subject: Re: RungeKuttta
> IÕm using a runge-kutta 4th order to simulate a vehicle
moving in a
> circle. I create an acceleration vector that is ALWAYS
perpendicular
> to the velocity vector. However, I noticted that the
velocity vector
> of my vehicle is slowly getting bigger. Why is the velocity
vector
> getting bigger if the acceleration that IÕm giving to rk4
is always
> orthogonal with the velocity? Anybody know if there is a
better
> integrator to use for this problem?
Not quite enough info for suggesting a suitable fix. In the
absence of
that, this is what you could do:
a) Check integrator by running w/successive halving of the
step. If
divergence rate is not decreasing, you need a new integrator.
Or,
b) Converges, but too slowly, then you can (in this case)
introduce
external error control by renormalizing* velocity at each
integration
step.
* An old aerospace simulation trick when reducing steps, was
seldom an
option. e.g. quaternion kinematics.
--
Dr.B.Voh
------------------------------------------------------
Applied Algorithms http://sdynamix.com
===
Subject: Re: C routines for Special Functions
[Followups set to comp.lang.c]
<snip Error Handling
> IÕm not a fan of the Ôint status = \
gsl_function(...); if
(status)
> {...Õ method of error handling for functions that return
doubles.
> I am especially not a fan of having default error handlers
call
> abort(), like the assert macro in C. Of course IÕve had \
the
head
> of IT tear me a new asshole when a programmer of mine left
an assert
> in production code that got called 15 minutes before market
open,
> so I canÕt claim to be unbiased.
If the assertion was written properly (i.e. the condition was
one that
/must/ be true if the program is correctly written), then the
fault lies
not with the programmer who left an assert in production code
but with
whoever failed to test the condition that led to the
assertion being fired.
<snip>
--
Richard Heathfield : binary@eton.powernet.co.uk
Usenet is a strange place. - Dennis M Ritchie, 29 July 1999.
C FAQ: http://www.eskimo.com/~scs/C-faq/top.html
K&R answers, C books, etc: http://users.powernet.co.uk/eton
===
Subject: Approximating ,,phi using Newton method
Let sqrt[n](X):=X**{1/n}=X^{1/n} ,
(1) f(x)=x^2 -x -1 , I=(1/2,2] ,
and for j in {2,3,4,5,6} consider the functions F_j:I-->R ,
where
j | F_j(x)
2 | f(x)
---------------------------------------------
3 | f(x)/sqrt[2](2x-1)
---------------------------------------------
4 | f(x)/sqrt[3](3f(x)+5)
---------------------------------------------
5 | f(x)/sqrt[4]((2x-1)*(2f(x)+5))
---------------------------------------------
6 | f(x)/sqrt[5](f^2(x)+5f(x)+5)
For j in {2,3,4,5,6} the equations F_j(x)=0 , x in I ,
are equivalent: these equations have the same root, namely
phi = ( 1+sqrt[2](5))/2 .
In order to approximate phi we construct the iterative
functions
N_j(x) = x - F_j(x)/F_jÕ(x) , (Newton) , j in {2,3,4,5,6} ,
and then the iterations ( x_j(n) )_{n>=0} , where
(1) x_j(k+1) = N_j(x_j(k)) , k=0,1,... , (j in {2,3,4,5,6} )
with the same start-point x_j(0) = 2 . Therefore we have five
methods to approximate phi.
Which sequence converges faster to phi ?
===
Subject: Re: Can this be solved?
>Peter,
>Say I have this: y=x/(x+exp(c1-c2*x)) where c1 and c2 are
parameters
>that set the shape of the curve (see link
>www.gardenwithinsight.com/help100/00000467.htm for more
information)
>Using Excel, you may assume values for x and plot y.
>Now, how can I predict c1 and c2 values so I get a
satisfactory curve
>(utility curve)to my needs; instead of randomly changing
these two
>parameters?
> Mike
mike, depends on the data and how ambitious you are in doing
applied
statistics
right.
i assume: data(x(i),y(i)), i=1 to N
x(i) error free , errors in y(i) random, independent
from each other for any i, expectation value zero and of the
same variance
(in words: no systematic measurement error and all
measurements of the same
precision)
then standard least squares applies (is in EXCEL for example)
you minimize
sum_{i form 1 to N} (y(i)-x(i)/(x(i)+exp(c1-c2*x(i)))) ^2
with respect to c1, c2. there are nonlinear least squares
solvers of good
quality in the public domain, see
http://plato.la.asu.edu/topics/problems/nlolsq.html
you need reasonable initial values for this. use
log( (x(i)/y(i)) - x(i) ) = c1-c2*x(i)
for those i for which the argument of the log is positive:
either take 2 such i and solve the linear system in two
unknowns c1 and c2,
or
better, use all such i and fit a least squares line through
these points.
again there are linear squares solvers at the source cited
above, but for
this
task you need no initial guess.
hth
peter
===
Subject: Numerical integration of Rayleigh distribution
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision:
1.9 primary) id i13Df4B05491;
Hey,
I should numerically integrate function of the form int[0,inf]
(x/P*exp(-x^2/P2)*f(x))dx by using Gaussian quadrature. I am
not
familiar with numerical integration and I am quite unsure how
to
select the nodes. I am using Matlab.
Helka
===
Subject: Re: Numerical integration of Rayleigh distribution
>Hey,
>I should numerically integrate function of the form
int[0,inf]
>(x/P*exp(-x^2/P2)*f(x))dx by using Gaussian quadrature. I am
not
>familiar with numerical integration and I am quite unsure
how to
>select the nodes. I am using Matlab.
>Helka
you must use
(x/P)*exp(-(x/P)^2)
as a weight function on [0,inf] and compute the quadrature
rule from this
one.
I donÕt know whether GautschiÕs code ORTHPOL \
(nr 726 in
http://www.netlib.ornl.gov/toms
) allows infinite intervals.
It might be easier for you to use a ready to use code from
quadpack for
this. (DQAGI) code 691 in the same source.
this one transforms [0,inf] to [0,1] by
z=1/(x+1)
you get a boundary singularity and this is resolved by high
order
Gauss-Kronrod
integration. nothing easy for a beginner.
there is also a lecture notes book describing quadpack in
detail.
hth
peter
===
Subject: Re: Intel Fortran for Windows: student edition $29
[...]
>Has anyone tried Visual Fortran? What is it like - anything
like Delphi,
>Jbuilder, etc?
Not even close. The appelation ÔVisualÕ is so \
the product can
piggyback on
the appeal of MicrosoftÕs .NET IDE even though Intel Fortran
doesnÕt fully
integrate into that IDE: no drag and drop, context-sensitive
help,
intellisence, mixed-language dllÕs, etc. It can be used from
the command
line if you prefer. Both Salford and Lahey have Visual
Fortran products for
.NET that are more at home in .NET than Intel is ever likely
to be.
> I havenÕt tried the Intel Visual Fortran, but \
IÕve been
using Compaq
> Visual Fortran for 5 years or so, and Intel purchased the
compiler
> division of Compaq, and IVF is the follow-on.
> The answer to your question is that its not really like
Delphi.
> Basically, VF is the name brand Fortran 90/95 compiler for
windows.
> The VF product tightly integrates the compiler with the MS
Visual
> Studio IDE. Together they are a very nice package,
frankly--there are
> more than a few times that IÕm working in \
MatlabÕs IDE that
I wish I
> were working in CVF (conditional breakpoints for example)
> Anyway, its a very nice product, and the service is truly
> second-to-none.
MattÕs comments are valid wrt Compaq Visual Fortran and its
integration
into Visual Studio 98 (not .NET) although it doesnÕt supply \
a
mechanism for
containment of components with a visual interface. Alas, CVF
is no longer
under development.
In contrast to CVF, support for Intel Visual Fortran leaves a
lot to be
desired and youÕre now expected to pay for it.
--
Ciao,
Gerry T.
===
Subject: HELP!!! a set theory problem.
IÕve got a problem:
Given R is the set of real numbers.
define: (x1, y1) < (x2, y2) if x1 < x2 OR ( x1 = x2 AND y1 <
y2)
define: f is monotone increasing if f((x1, y1)) < f((x2, y2))
for every
(x1, y1) < (x2, y2)
Question: Is there such a monotone increasing function f
which is also a
bijection from R x R to R ??? Prove it.
----------
Any Idea? I think this could be as same as to prove you can
insert a set
of real numbers with the same cardinality of R between any
two real
numbers(with all values bigger than the first and smaller than
the
second). How can I prove this???
Song
===
Subject: Re: overdetermined linear system
>My program need to solve some overdetermined linear system,
and the
>coefficient matrix could or not be a compatible matrix. I
need solution
in
>both situation.
>Anyone know an efficient method to solve this systems?
>Someone tell me about an Eremin iterative algorithm who try
to minimize
the
>Chebyshev norm. If you know something about this topic you
could tell me
>somthing like...where to find texts who explain that or also
only give me
>the algorithm.
>..and sorry for my english!
>Alberto Gori
if you need only a solution of an overdetermined linear
system then
linear least squares is the irst idea and algorithms can be
found at
http://plato.la.asu.edu/topics/problems/nlolsq.html
but you write about minimize the chebyshev norm. with nothing
better
available you can do this using linear programming:
minimize x(n+1) subject to
-x(n+1) <= sum_{j=1 to n} a(i,j)*x(j)-b(i) <= x(n+1) for
i=1,....,M
but there are better approaches, like the (discrete)
Remez-Algorithm
(in case of a specail matrix A representing a Haar system) or
the specailized codes COCA and CAPROX
at
http://plato.la.asu.edu/topics/problems/approx.html
hth
peter
===
Subject: Re: overdetermined linear system
> My program need to solve some overdetermined linear system,
and the
> coefficient matrix could or not be a compatible matrix. I
need solution
in
> both situation.
> Anyone know an efficient method to solve this systems?
> Someone tell me about an Eremin iterative algorithm who try
to minimize
the
> Chebyshev norm. If you know something about this topic you
could tell me
> somthing like...where to find texts who explain that or also
only give me
> the algorithm.
> ..and sorry for my english!
> Alberto Gori
Try Singular Value Decomposition (SVD). I found lots of hits
on Google
for SVD and overdetermined systems, e.g.,
http://cloe.ucsd.edu/claudia/saratov96/node14.html
Good luck,
OUP
===
Subject: Re: overdetermined linear system
One Usenet Poster <me@my.computer.org> ha scritto nel
messaggio
> My program need to solve some overdetermined linear system,
and the
> coefficient matrix could or not be a compatible matrix. I
need
solution
in
> both situation.
> Anyone know an efficient method to solve this systems?
> Someone tell me about an Eremin iterative algorithm who try
to minimize
the
> Chebyshev norm. If you know something about this topic you
could tell
me
> somthing like...where to find texts who explain that or also
only give
me
> the algorithm.
> ..and sorry for my english!
> Alberto Gori
> Try Singular Value Decomposition (SVD). I found lots of
hits on Google
> for SVD and overdetermined systems, e.g.,
> http://cloe.ucsd.edu/claudia/saratov96/node14.html
> Good luck,
> OUP
SVD is not my solution because it couldnÕt work with an
incompatinle
coefficient matrix!
I need solution in every situation!
===
Subject: Faultisolating Kalman GPS/INS Multisensor Navigation
5 day 40 hour course instructed by Dr. Sam C. Bose of
Technalytics
http://www.technalytics.com/course.htm
for detailed course outline, instructor qualifications,
previous
attendee comments, lodging and registration info.
Ms Vlasta Stastna
Course Coordinator
===
Subject: Re: 2D discrete hartley transform
> IÕm trying (unsuccessfully) to find some C \
code that can
perform a 2D-DHT
on
> real-valued data obtained from a 2D image. Does anyone have
workable code
> they are willing to share? IÕve been plugging away with \
the
implementation
> of fftw_dht
> in fftw, but IÕm finding it more than \
challenging to sort
through their
> code.
This page has C code from two different authors. The
first set appears to have been written for use with BMP
image files. I havenÕt actually run this code, \
so I canÕt
comment on itÕs quality. Good luck.
http://www.jjj.de/fft/fftpage.html
- Brian
===
Subject: Re: 2D discrete hartley transform
> This page has C code from two different authors. The
> first set appears to have been written for use with BMP
> image files. I havenÕt actually run this code, \
so I canÕt
> comment on itÕs quality. Good luck.
> http://www.jjj.de/fft/fftpage.html
Ah, good, it looks like I mis-spoke; both of these codes
implement the
non-separable transform (by a post-processing pass after
computing the
1D DHTs of the rows and columns).
(IÕm still curious to learn why people find this \
transform
useful.)
===
Subject: Re: 2-D convolution using FFTW3
I have found the solution for this convolution.
It is simple.
step 1: padding the two input matrix with zero.
nExt = nx + ny -1
step 2: convert the two array using FFT forward,
step 3: multiply the two converted array component by
component.
step 4: IFFT the multiply array
> Hi everyone,
> I want to do a 2D real discrete convolution using FFTW3,
can anybody
hel me?
> Jason Deng
===
Subject: Re: Nonlinear curve fit for saturation data
>> I need an automated method to fit a curve of the form y =
>> c*x*(1+a*e^(b*x)) to test data.
>> Voc(PU) Igf(A)
>> 0.196 5.00
>> 0.304 7.85
>> 0.399 10.43
>> 0.616 16.20
>> 0.797 21.54
>> 0.906 25.89
>> 0.942 29.01
>> 1.014 32.75
>> 1.051 37.25
>> 1.105 44.80
>> 1.196 61.18
>> 1.304 105.61
>> This is approximately fit with:
>> a = 0.000232
>> b = 7.005
>> c = 25.74
>Fitting to lowest sum of squared relative error, I get:
>a = 2.45216969e-004
>b = 6.96242410e+000
>c = 2.56948865e+001
>which seemed to me the best result of the fitting
>targets I tried. Fitting to lowest sum of squared
>absolute error, I get:
IsnÕt this just the usual least squares fit, as \
opposed to your
immediately previous
example, where used the relative error for the fit (I got the
same numbers
you did for
this case, using relative error). If it is, your results seem
to be
different from what
myself and others got for a standard least squares fit. Have I
misunderstood what you
mean by fitting to lowest sum of squared absolute error?
>a = 0.02081207
>b = 4.03580003
>c = 15.44149065
>which is rather a bit different. IÕm adding this equation
>to my curve fitting site now.
> James Phillips
> http://zunzun.com