mm-1006
===
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?
>Eric
>>YouÕre probably using an explicit RK4.
>>Explicit RK schemes do not conserve energy.
>>That is probably the reason for your troubles.
>>You should try a symplectic implicit RK.
> Simplectic methods do not necessarily
> conserve energy. Neither condition implies
> the other.
Could you precise?
What about simplectic AND symmetric schemes?
===
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?
Eric
>YouÕre probably using an explicit RK4.
>>Explicit RK schemes do not conserve energy.
>>That is probably the reason for your troubles.
>>You should try a symplectic implicit RK.
> Simplectic methods do not necessarily
> conserve energy. Neither condition implies
> the other.
> Could you precise?
> What about simplectic AND symmetric schemes?
Symmetry is not enough. Counterexample: the trapezoidal rule
is
symmetric but not symplectic.
Of the methods that are symplectic and directly preserve
energy,
without
projection onto the Hamiltonian (H) surface, the fixed step
implicit
midpoint rule (IMR) is the simplest for *linear* problems (ie,
quadratic H). We use a nonlinear variant of it for
geometrically
nonlinear structural mechanics (e.g, simulation of aircraft
maneuvers) and it works quite well. I can imagine the IMR
will do
fine in the OP problem.
Recommended reading: Chapter 8 of Dynamical Systems and
Numerical
Analysis, by Stuart and Humphries, Cambridge, especially 8.5.
The basic theorem is: it is not possible to find *general *
methods
that
preserve both properties for arbitrary H. Proof: Zhong and
Marsden,
1998.
===
Subject: Re: Which finite elements would you recommend for
this problem?
> For linear advection-diffusion equation, an
Eulerian-Lagrangian
> type of Galerkin method should be suitable. These will
rather
BjÀrn-Ove, John,
Yours andy.
===
Subject: FFT with normal distribution
I get a N-bit pseudo-random sequence from a very good
generator.
I transform the sequence using the FFT.
I discard the first component (which is always real) and the
second half of
the components (which are identical to the first half). So I
have N/2-1
components from 1 to N/2.
Both components seem normally distributed. I seen that by a
graphical
display (binning the numbers) and by the Kolmogorov-Smirnov
test.
I calculate the KS test for M N-bit sequences (N=16 kbit).
The M p-values gotten from the imaginary components are
uniformly
distributed and this is what it is expected because the M
sequences are
random.
But the M p-values gotten from the real components are not
uniformly
distributed. Why?
Cristiano
===
Subject: high order polynomial interpolation
http://www.pandasoftware.com/PAA
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
===
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).
Alan Miller has quadruple precision Fortran 90 code to solve
polynomial equations -- see q_pzeros.f90 at
http://www.users.bigpond.net.au/amiller/quad.html .
===
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
===
Subject: Re: high order polynomial interpolation
> 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).
This is quite unstable numerically. Computing the coefficients
of a power representation usually deteriorates the condition
of a problem significantly.
Rather evaluate the polynomial at N points x i in the region
of
interest and create the Newton interpolant.
The zero-finding problem can then be written as an eigenvalue
problem with a sparse Hessenberg matrix (analogous to deriving
the companion matrix). Find its zeros z i approximately,
and repeat the same procedure with the z i in place of the x
i.
If necessary, repeat again.
This gives you all zeros at the accuracy that can be expected
from a sensitivity analysis of your original representation
(whatever it is).
If N is large and only some eigenvalues are wanted, it is
enough
to recompute k function values at the approximate zeros in
the region
of interest, keeping N-k other function values unrefined.
To save work in the eigensolver, make sure that your
Hessenberg
matrix is of the form that the solver tries to establish in
the
first place (else transpose the matrix).
Other methods (not using eigenvalue techniques) can be found
at
A. Neumaier, Enclosing clusters of zeros of polynomials,
http://www.mat.univie.ac.at/~neum/papers.html#polzer
[~ is a tilde]
A. Sch.8afer and A. Neumaier, Divided differences,
shift transformations, and LarkinÕs root finding \
method,
Math. Comput. 45 (1985), 181-196.
(download N42 in
http://www.mat.univie.ac.at/~neum/publist.html )
Arnold Neumaier
===
Subject: exact solution of an inverse matrix
Distribution: inet
I am looking for an exact solution of an inverse matrix (not
from Hilbert
matrix
nor diagonal ones ) and which not ill conditioned
If someone has an idea
Karim.
===
Subject: Re: exact solution of an inverse matrix
Distribution: inet
> I am looking for an exact solution of an inverse matrix
And youÕre not just looking for the solution to a linear
system?
Try linking lapack to an infinite precision arithmetic
package. I doubt
that the outcome would be useful, but thatÕs what \
IÕd try.
V.
--
email: lastname at cs utk edu
homepage: cs utk edu tilde lastname
===
Subject: Re: exact solution of an inverse matrix
Distribution: inet
> I am looking for an exact solution of an inverse matrix
(not from
> Hilbert matrix nor diagonal ones ) and which not ill
conditioned
> If someone has an idea
> Karim.
of order N
===
Subject: Re: exact solution of an inverse matrix
> I am looking for an exact solution of an inverse matrix
(not from Hilbert
> matrix nor diagonal ones ) and which not ill conditioned
> If someone has an idea
> Karim.
what do you mean with exact? which dimension?
Michael
--
Remove the sport from my address to obtain email
www.enertex.de - Innovative Systeml.9asungen der Energie- und
Elektrotechnik
===
Subject: Re: exact solution of an inverse matrix
> I am looking for an exact solution of an inverse matrix
(not from
Hilbert
> matrix nor diagonal ones ) and which not ill conditioned
> If someone has an idea
> Karim.
> what do you mean with exact? which dimension?
> Michael
> --
> Remove the sport from my address to obtain email
> www.enertex.de - Innovative Systeml.9asungen der Energie-
und
Elektrotechnik
Create a matrix of rank 1, that is, of the form
M = I - v otimes v^dag
where v is an N-component vector and I is the NxN unit matrix.
Then
M A = I
becomes
A - v v^dag A = I
so that
A = I + v v^dag A
However,
(1 - v^2) v^dag A = v^dag
so that
v v^dag
A = I + ----------
1 - v^2
--
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: exact solution of an inverse matrix
Distribution: inet
>I am looking for an exact solution of an inverse matrix (not
from
Hilbert
>matrix nor diagonal ones ) and which not ill conditioned
>If someone has an idea
>Karim.
>>what do you mean with exact? which dimension?
>>Michael
>>--
>>Remove the sport from my address to obtain email
>>www.enertex.de - Innovative Systeml.9asungen der Energie-
und
Elektrotechnik
> Create a matrix of rank 1, that is, of the form
> M = I - v otimes v^dag
> where v is an N-component vector and I is the NxN unit
matrix.
> Then
> M A = I
> becomes
> A - v v^dag A = I
> so that
> A = I + v v^dag A
> However,
> (1 - v^2) v^dag A = v^dag
> so that
> v v^dag
> A = I + ----------
> 1 - v^2
===
Subject: Re: Intel Fortran for Windows: student edition $29
>>While browsing ProgrammerÕs Paradise, I saw that a student
edition of
>>Intel Visual Fortran for Windows 8.0 is available for $29.
The link is
>>http://www.programmersparadise.com/Product.pasp?txtCatalog=
Paradise&txtCat
e
>gory=&txtProductID=I23+0D03
>> This product should give Matlab a run for its money, and
>> run your problems in seconds instead of days!
>Has anyone tried Visual Fortran? What is it like - anything
like Delphi,
>Jbuilder, etc?
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.
===
Subject: Re: real analysis
>A polynomial P of degree n (n> k)has the property that
> P(c)=...=P(k)(c)=0, where P(j)(c) is the jth derivative of
P at c.
>Show that
>P(c)=(x-c)(k+1)Q(x)
>{(x-c)(k+1)represents (x-c)to the power (k+1)}
>whereQ(x) is a polynomial of degree n-k-1.
homework
hint: Taylor
hth
peter
===
Subject: Re: Can this be solved?
>Can this equation be solved?
>y=a+b*exp(- c * x)
>At x=L, y=0 and
> x=M, y=Alpha (where 1> Alpha>0)
> x=H, y=1
>where H>M>L>0
>My objective is to find a, b, c
>Any reference would be appreciated? Online link is fine \
too?
>Mike
> homework?
> three equations in three unknowns, well , nonlinear
> insert x=L, get a dependent on b and c,
> insert x=M , expressing now b as a function of c
> insert x=H
> 1/Alpha = (exp(-c*H)-exp(-c*L))/(exp(-c*M)-exp(-c*L))
> solve for c numerically: http://www.netlib.org/fmm/fzero.f
> c>0 very small gives on the right hand side approximately
> (H-L)/(M-L) > 1
> and c to infinity gives 1 on the right hand side. so at
least for
> (H-L)/(M-L) > 1/Alpha
> there is a solution
> hth
> peter
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
===
Subject: Re: Can this be solved?
>Can this equation be solved?
>y=a+b*exp(- c * x)
>At x=L, y=0 and
> x=M, y=Alpha (where 1> Alpha>0)
> x=H, y=1
>where H>M>L>0
>My objective is to find a, b, c
>Any reference would be appreciated? Online link is fine \
too?
>Mike
> homework?
> three equations in three unknowns, well , nonlinear
> insert x=L, get a dependent on b and c,
> insert x=M , expressing now b as a function of c
> insert x=H
> 1/Alpha = (exp(-c*H)-exp(-c*L))/(exp(-c*M)-exp(-c*L))
> solve for c numerically: http://www.netlib.org/fmm/fzero.f
> c>0 very small gives on the right hand side approximately
> (H-L)/(M-L) > 1
> and c to infinity gives 1 on the right hand side. so at
least for
> (H-L)/(M-L) > 1/Alpha
> there is a solution
> hth
> peter
Other than that, are you aware of an s-curve formula that
gives
y(a,b,c,x) value? I went into the above formula after I gave
up on
finding an s-curve formula?
Mike
===
Subject: Re: C routines for Special Functions
>> GSL ties you into a number of bad design decisions
>> that make it clumsy to use.
>Would you care to elaborate?
> Of course this is just my opinion. I think the GSL group is
doing
> a great service and since the code is GPL I should really
be rolling
> up my sleeves and coding instead of complaining. :-)
I second this.
> When the first example program in their documentation \
defines
Ôint
> main(void) { ... }Õ you have to wonder a bit. Why the
ÔvoidÕ? It
> is perfectly legal, just odd looking and unnecessary. When
I was
> evaluating this for use in a production library at a large
Equity
> Derivatives firm a couple of years ago, I kept \
finding ÔoddÕ
things.
> Ultimately I chose a Fortran library because GSL was not
sufficiently
> mature at that time. It has definitely come a long way since
then.
> If you belive a modern library should be written in C
instead of C++,
> you can stop reading now. Many of my qualms are due to
issues relating
> to the gyrations GSL goes to because it does not use C++.
IÕve taken
> ßak over the years for insisting on writing certain
libraries in C (and
> even more ßak for using Fortran) but I have come to the
conclusion that
> the small audience that insists on C is, well, small. IÕm
open to data
> indicating otherwise.
<SNIP>
There are two G95 projects (long story), which are trying to
build a
free Fortran 95 compiler. The sites are
http://gcc.gnu.org/fortran/
and http://g95.sourceforge.net/ . Once they exist, I think
Fortran 95
would be a good language for something like the Gnu Scientific
Library. There would be no need to define special vector,
matrix, and
complex number classes, since these features are already part
of the
interoperability with C will be standardized. Most of the
algorithms
in the GSL already exist as public domain Fortran 77 code at
Netlib.
If that code were modernized by using features of Fortran 95,
it might
be used even more than it already is.
===
Subject: Re: C routines for Special Functions
> Of course this is just my opinion. I think the GSL group is
doing
> a great service and since the code is GPL I should really
be rolling
> up my sleeves and coding instead of complaining. :-)
> When the first example program in their documentation \
defines
Ôint
> main(void) { ... }Õ you have to wonder a bit. Why the
ÔvoidÕ?
> It is perfectly legal, just odd looking and unnecessary.
Hmmm, odd looking and unnecessary and ... part of the C
standard ;-)
e.g. see ANSI/ISO/IEC 9899:1990 section 5.1.2.2.1 for int
main(void)
> If you belive a modern library should be written in C
instead of C++,
> you can stop reading now.
Fair enough. ItÕs a GNU project so the main reason we used C
is
because of the GNU coding standards (which specify that C be
used
where possible). At the time GSL was started (over seven
years ago)
C++ techniques such as template metaexpressions and traits
were just
being developed. I am still waiting to see a good general
numerical
library in C++ however.
(We recommend using C with a very high-level object oriented
language
such as Lisp or Python.)
> IÕm not a fan of the Ôint status = \
gsl_function(...); if
(status)
> {...Õ method of error handling for functions that return
doubles.
IÕm not a particular fan of it either, but it is the
conventional way
to handle it in C. We worked on a very ßexible
dynamically-scoped
error context system for a while but it was not popular.
> I am especially not a fan of having default error handlers
call
> abort(), like the assert macro in C.
Without the abort() on error default we had many users
blithely
ignoring errors and complaining that GSL was producing
incorrect
results. So we had to make abort() the default.
It is done with a single configurable gsl_error() function so
it can
be changed globally.
> GSL does have Ônatural formÕ calling \
conventions available
(more
> code on their side to maintain) but they do not permit error
> checking.
These should return NaN/Inf for error/overßow conditions (if
there
are any that donÕt, it is bug).
> I prefer returning NaNÕs that have embedded error
information.
We have to be portable so this is not really an option for a
core
feature (?).
> Excessive use of C macros in source code. Makes the code
hard to
> understand/fix/maintain.
I agree the macros for the implementation of the different
types of
vectors and matrices are inelegant. Possibly it would be
better to
use a scripting language or m4 to generate the code, rather
than cpp.
> Not modular. Dependencies between modules should be
explicit.
Tried, but turns out to be difficult in practice (originally
all
modules were independent) and users seem not particularly
interested.
Probably one could write a script to determine the
dependencies
automatically, from the source, #includes or the .a file.
> At any rate, thatÕs my 2 cents. Feel free to disagree, \
itÕs
an internet
> newsgroup after all. Comments are welcome.
I would say that we did have to make some design decisions,
but
hopefully they are not intrinsically bad, just a question of
not being
able to please all of the people all of the time.
--
Brian Gough
(GSL Maintainer)
Network Theory Ltd
http://www.network-theory.co.uk/
===
Subject: Re: C routines for Special Functions
[fÕups to comp.lang.c]
> Of course this is just my opinion. I think the GSL group is
doing
> a great service and since the code is GPL I should really
be rolling
> up my sleeves and coding instead of complaining. :-)
> When the first example program in their documentation \
defines
Ôint
> main(void) { ... }Õ you have to wonder a bit. Why the
ÔvoidÕ?
> It is perfectly legal, just odd looking and unnecessary.
> Hmmm, odd looking and unnecessary and ... part of the C
standard ;-)
> e.g. see ANSI/ISO/IEC 9899:1990 section 5.1.2.2.1 for int
main(void)
More so, C99 has deprecated ()...
6.11.6 Function declarators
The use of function declarators with empty parentheses (not
prototype-
format parameter type declarators) is an obsolescent feature.
--
Peter
===
Subject: 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
===
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.
this is the decision concerning the rank of your matrix, so
you need to
know
how much larger the eigenvalue lr will be above the roundoff
level.
you can assume that the inßueneces of rounding errors on the
determination
of
eigenvalues of a symmetric mtrix is like the inßuence of a
symmetric
eps-perturbation in the original matrix. and there is a
theorem wich says
if A and B are symmetric and the eigenvalues of A and B are
ordered in
magnitude
then the difference
abs(eigenvalue_i(A)-eigenvalue_i(B)) <= norm(A-B)
hence in the order of some units in the last digit, that means
n*eps*norm(A)
for example in this application.
hth
peter
===
Subject: 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.