mm-1100
===
Subject: Re: Computational Evolution without velocities
Cc: p.bangert@iu-bremen.de
>
> I have heard the folklore that when one wants to simulate
some structure
> Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
> update formula
>
> x_new = x_old - Force DeltaTime / 2 Mass
>
> holds for any simulation and we do not need to take care of
velocities.
This equation is quite obviously wrong, for the simple fact
that it
does not even have consistent units. The left hand side has
dimensions
of [Length]; the right hand side is a length minus a
velocity. You cant
subtract a velocity from a length and get anything meaningful.
-- Gordon D. Pusch
===
Subject: Need help with steepest descent
Would someone care to explain how to use the steepest descent
algorithm to minimize a function, say, F(x) = x1 * x1 + 25 *
x2 * x2,
-5 < x1,x2 < 5 :-). Links to web tutorials would also be
appreciated.
And no, this is not my homework, my homework is to minimize
Rosenbrocks function :-). (I, uh, fell aspeep during the
class :-(
and now I cant find any references on the net.
Andrey
===
Subject: Re: Small sample size and logit model
>I have a sample of 6 observations and due to the
conditions I cannot
>make it bigger. Can I use any statistical models for such
a sample? I
>donot want to make a prediction, I just want to describe
the past
>period for wihch I have data. Can I use regresion? Can I
use logit?
>Perhaps you can recommend some books where this issue
would be
>discussed (small sample size, but I am NOT a
mathematician, so I need
>a rather easy-to-read book).
>
> in the extreme case, with 6 parameters (what is logit?)
> you might desribe the six data points exactly, neglecting
the data
errors.
> the rule of thumb is sqrt(of number of data
points/2)>=number of
parametrs.
> here : 1, maybe 2.
> software: see
> http://plato.la.asu.edu/topics/problems/nlolsq.html
>
> hth
> peter
variable takes only one of two possible values - zero or one.
You
estimate the model using the maximum likelihood method.
To ensure my understanding of your post I will ask one more
question.
Can I apply least squares method (preferably maximum
likelihood
method) and just ignore t-test? Can I such an estimation put
into my
ph.d.?
===
Subject: Re: Computational Evolution without velocities
> I have heard the folklore that when one wants to simulate
some structure
> Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
> update formula
> x_new = x_old - Force DeltaTime / 2 Mass
Youve got to stop reading those textbooks on cartoon
physics.
Tom Davidson
Richmond, VA
===
Subject: Re: Computational Evolution without velocities
>
> I have heard the folklore that when one wants to simulate
some structure
> Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
> update formula
>
> x_new = x_old - Force DeltaTime / 2 Mass
>
> holds for any simulation and we do not need to take care of
velocities.
What
> I would like to know is does anyone know whether this holds
only under
> certain conditions (evolution at equilibrium, etc.) and/or
knows of
research
> papers that deal with this issue and perhaps give a proof
of this claim as
a
> theorem. The strange thing is that apparently this is very
well known and
> widely practised but no one Ive spoken to knew where to
find the
details.
> Pat
>
> Free your mind. There is no spoon.
> ************************************************
> Dr. Patrick Bangert
> http://www.knot-theory.org
> Research Instructor for Mathematics
> International University Bremen
>
>
If F is constant, the above formula will give a constant
incremental
distance, i.e a constant velocity. In the absence of drag,
that aint
newtonian. Perhaps there was an assumption of a viscous
retarding force.
--
Joe Legris
===
Subject: Re: Need help with steepest descent
> Would someone care to explain how to use the steepest
descent
> algorithm to minimize a function, say, F(x) = x1 * x1 + 25
* x2 * x2,
> -5 < x1,x2 < 5 :-). Links to web tutorials would also be
appreciated.
Try adding the word gradient to your Google searches, and
youll
probably come up with better hits. The basic idea is that
from some
starting point near your minimum (or from some arbitrary
starting point if
you have a nice convex function like the above) you
repeatedly find the
functions gradient and choose a new point in the opposite
direction.
---
Roy Stogner
===
Subject: Help solving equation.
I am trying to solve the following equation.
E(z) +(H/(z^2) - B).E(z)=(2*H)/(Z^3)* ...
L
S E(t).dt. (S is the integral sign).
z
Where H and B are both constants.
My Boundary conditions are E(1.1) = -0.45 and E(500) = 0.
Any assistance would be greatly appreciated. If you could
e-mail me at
andrewl at global.net.au I would appreciate it.
Andrew Lindsay
===
Subject: Re: Need help with steepest descent
[[ This message was both posted and mailed: see
the To, Cc, and Newsgroups headers for details. ]]
>
> Would someone care to explain how to use the steepest
descent
> algorithm to minimize a function, say, F(x) = x1 * x1 + 25
* x2 * x2,
> -5 < x1,x2 < 5 :-). Links to web tutorials would also be
appreciated.
> And no, this is not my homework, my homework is to minimize
> Rosenbrocks function :-). (I, uh, fell aspeep during the
class :-(
> and now I cant find any references on the
net.
For a picture of a descent path down the Rosenbrock banana,
see the
picture in the upper left corner of
Clicking on it brings up a larger version (without the
descent path),
and a link to a PDF file which contains still higher-detail
pictures
(including the descent path), and a description of the
Mathematica code
which generated it.
--Ron Bruck
===
Subject: Re: Help solving equation.
>
> I am trying to solve the following equation.
>
> E(z) +(H/(z^2) - B).E(z)=(2*H)/(Z^3)* ...
>
> L
> S E(t).dt. (S is the integral sign).
> z
Im assuming that your equation reads (in Mathematica
notation)
eqn = (h/z^2 - b) e[z] + z e[z] ==
(2 h Integrate[e[t], {t, z, l}])/z^3
By differentiating this equation, one can eliminate the
integral from
the resulting pair of equations leading to a _third_ order
differential
equations.
> Where H and B are both constants.
>
> My Boundary conditions are E(1.1) = -0.45 and E(500) = 0.
I assume that l = 500 then?
You obtain a third boundary condition by evaluating your
equation for
z=l (where the integral vanishes):
(h/l^2 - b) e[l] + l e[l] == 0
Paul
--
Paul Abbott Phone: +61 8 9380 2734
School of Physics, M013 Fax: +61 8 9380 1014
The University of Western Australia (CRICOS Provider No
00126G)
35 Stirling Highway
Crawley WA 6009 mailto:paul@physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul
===
Subject: Grouping, non-repeating
Lets say youre going to have a dinner party
group of 24
people (12
couples). The group plans to get together monthly to have
parties.
The group has 2 rules: 1. Every couple must host at least
once. 2.
Any coulples dining together one month should not be paired
again.
For instance, in January
couple A hosts couples B, C, D
E hosts F, G, H
I hosts J, K, L
In February (and there after) couple A should never be in the
same
house as couples B or C, and A should never have to host
again.
So my question is this. Is there a generic equation (or
relationship)
that could tell me the minimum number of months and/or group
sizes
that must exist in order for the rules to stand up?
===
Subject: Re: Help solving equation.
>
>
> I am trying to solve the following equation.
>
> E(z) +(H/(z^2) - B).E(z)=(2*H)/(Z^3)* ...
>
> L
> S E(t).dt. (S is the integral sign).
> z
>
> Im assuming that your equation reads (in Mathematica
notation)
>
> eqn = (h/z^2 - b) e[z] + z e[z] ==
> (2 h Integrate[e[t], {t, z, l}])/z^3
That should have read
eqn = e[z] + (h/z^2 - b) e[z] == 2 h/z^3
Integrate[e[t],
{t, z, l}]
> By differentiating this equation, one can eliminate the
integral from
> the resulting pair of equations leading to a _third_ order
differential
> equations.
>
> Where H and B are both constants.
The equation reads
b (3e[z] + z e[z]) == (3h e[z])/z^2+(h
e[z])/z+ 3e[z]+ z
e[z]==0
Mathematica returns a (formal) closed-form solution for this
equation in
terms of 1F2 hypergeometric functions involving roots of the
cubic
z^3 + h z - z + 3 h
> My Boundary conditions are E(1.1) = -0.45 and E(500) = 0.
>
> I assume that l = 500 then?
>
> You obtain a third boundary condition by evaluating your
equation for
> z=l (where the integral vanishes):
>
> (h/l^2 - b) e[l] + l e[l] == 0
And that should have read
e[l] + (h/l^2 - b) e[l] == 0
Paul
--
Paul Abbott Phone: +61 8 9380 2734
School of Physics, M013 Fax: +61 8 9380 1014
The University of Western Australia (CRICOS Provider No
00126G)
35 Stirling Highway
Crawley WA 6009 mailto:paul@physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul
===
Subject: Re: Grouping, non-repeating
Related problems are the progressive party problem and
the social golfer problem. Google will help you find papers
on these problems. These are actually very difficult
problems to solve.
-------------------------------------------------------------
---
Erwin Kalvelagen
erwin@gams.com, http://www.gams.com/~erwin
-------------------------------------------------------------
---
> Lets say youre going to have a dinner
party group of 24
people (12
> couples). The group plans to get together monthly to have
parties.
> The group has 2 rules: 1. Every couple must host at least
once. 2.
> Any coulples dining together one month should not be paired
again.
> For instance, in January
> couple A hosts couples B, C, D
> E hosts F, G, H
> I hosts J, K, L
> In February (and there after) couple A should never be in
the same
> house as couples B or C, and A should never have to host
again.
> So my question is this. Is there a generic equation (or
relationship)
> that could tell me the minimum number of months and/or
group sizes
> that must exist in order for the rules to stand up?
===
Subject: Re: Computational Evolution without velocities
* Gordon D. Pusch:
>> Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
>> update formula
>>
>> x_new = x_old - Force DeltaTime / 2 Mass
>>
>> holds for any simulation and we do not need to take care
of velocities.
> This equation is quite obviously wrong, [...]
He may be talking about Stoermer-Verlet algorithms,
which have update formulae that are indeed relatively
simple and dont contain the velocity. Also, they are
used for n-body problems and hard sphere ßuids and have
other interesting properties (look up symplectic algorithms).
x_new = 2 x_old - x_older + DeltaTime^2 Force / Mass +
O(DeltaTime^4)
would be a correct formula. That is not far off from his
version, and given that he heard it as folklore, it may
well be it.
(Of course, the velocity is not neglected here. It is
implicitly
contained in the ODE of 2nd order that is discretized
directly.
It could also be rewritten into two ODEs of 1st order and
discretized
afterwards, which would lead to better known schemes like
Runge-Kutta.)
--
Rupert
===
Subject: Re: Computational Evolution without velocities
Sorry, there was a typo! I should have squared DeltaTime.
This is by no means cartoon physics. People are using this in
research
simulations and Im wondering if there is a sound basis or
not.
Free your mind. There is no spoon.
************************************************
Dr. Patrick Bangert
http://www.knot-theory.org
Research Instructor for Mathematics
International University Bremen
> I have heard the folklore that when one wants to simulate
some structure
> Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
> update formula
> x_new = x_old - Force DeltaTime / 2 Mass
> holds for any simulation and we do not need to take care of
velocities.
What
> I would like to know is does anyone know whether this holds
only under
> certain conditions (evolution at equilibrium, etc.) and/or
knows of
research
> papers that deal with this issue and perhaps give a proof
of this claim
as
a
> theorem. The strange thing is that apparently this is very
well known and
> widely practised but no one Ive spoken to knew where to
find the
details.
> Pat
> Free your mind. There is no spoon.
> ************************************************
> Dr. Patrick Bangert
> http://www.knot-theory.org
> Research Instructor for Mathematics
> International University Bremen
===
Subject: Wavelets package in Java
Has anyone ever used a wavelet (de)composition package under
Java? Or
a program that can be used with Java (looked at Jython, but I
didnt
find a wavelet package for that)?
Please mail replies also to gez_75@hotmail.com
Gez
===
Subject: Re: Computational Evolution without velocities
The underlying assumptions allowing to neglect the velocity
term is that it always remains negligible with respect to the
force
term. This situation occurs effectively in very viscous ßuids,
to random ßuctuations.
Dan
> Sorry, there was a typo! I should have squared DeltaTime.
>
> This is by no means cartoon physics. People are using this
in research
> simulations and Im wondering if there is a sound basis or
not.
>
>
>
> Free your mind. There is no spoon.
> ************************************************
> Dr. Patrick Bangert
> http://www.knot-theory.org
> Research Instructor for Mathematics
> International University Bremen
>
>
>>I have heard the folklore that when one wants to simulate
some structure
>>Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
>>update formula
>>x_new = x_old - Force DeltaTime / 2 Mass
>>holds for any simulation and we do not need to take care of
velocities.
>
> What
>
>>I would like to know is does anyone know whether this holds
only under
>>certain conditions (evolution at equilibrium, etc.) and/or
knows of
>
> research
>
>>papers that deal with this issue and perhaps give a proof
of this claim
as
>
> a
>
>>theorem. The strange thing is that apparently this is very
well known and
>>widely practised but no one Ive spoken to knew where to
find the
details.
>>Pat
>>Free your mind. There is no spoon.
>>************************************************
>>Dr. Patrick Bangert
>>http://www.knot-theory.org
>>Research Instructor for Mathematics
>>International University Bremen
>
>
>
===
Subject: Coupled system of non-linear PDE
Hi everybody,
Im trying to solve a non-linear PDE system in u(x,t) e
v(x,t) of
the form:
-------------------------------------------------------------
---------------
--------
diff(u(x,t),t)-a(u,v)*diff(u(x,t),x,x)-b(u,v)*diff(v(x,t),x,x
);
diff(v(x,t),t)-c(u,v)*diff(u(x,t),x,x)-d(u,v)*diff(v(x,t),x,x
);
with Neumann BC on x=0 e x=L. The initial conditions are:
u(x,0) = {u_1 for x R^3 be a continuous path and a < c < d <
b. Let
C={p(t) | c < t < d }. Must p^-1 (C) be path-conneted?
===
Subject: analysis-connected
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision:
1.9 primary) id h8NFubw11851;
. Let A is a subset of R^2 be path- conneted. Regarding A as
a subset
of the xy-plane in R^3, show that A is still path-conneted.
Can you
make a similar argument for A connected?
===
Subject: suggest?
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision:
1.9 primary) id h8NKTq430748;
I dont know for sure, but the way I represented primes in
the
posting prime number distribution - solved may help you -
Ive
been planning on trying to use it to attack the zeta fcn
myself.
Ive kind of hit a roadblock, so please write something
there
when
you check it out. It only uses (maybe) fourth grade math
until I
give a more rigorous treatment, and even that part isnt
really that
hard. Heres a taste of the way it represents primes:
105-2^x=y ,(integer x, y<121)
y=103,101,97,89,73,41,-23
105+2^x=y
y=107,109,113
for larger y, the eqn is different. To get the rest of the
primes<121, need a bit more, see the posting.
Aaron
===
Subject: Program code for finding eigenvector in Casio
Calculator
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision:
1.9 primary) id h8OJQWL23555;
hi,
anyone know how to write or solve the eigenvector in Casio
9850GB
model. i know hp and TI do. but not in casio 9850GB.
if someone know,pls email to me.
that will be very helpful.
thx
===
Subject: Confidence interval for intra class correlation
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum,
$Revision:
1.9 primary) id h8S4V3x01889;
I am analyzing inter-observer variabilty of a Continuous
variable by
using intra-class correlation. Can somebody tell me how to
calculate
the 95% confidence interval for this coefficient.
basu
===
Subject: Re: Solving Nonlinear Equation
>
> Hi
> I am chemical engineer and I am solving some nonlinear
equations
> F = F(x)
To solve F(x)=0, apply an optimization routine to f(x) = sum
F_i(x)^2.
Many excellent routines are available on the web, see, e.g.,
http://www.mat.univie.ac.at/~neum/glopt/software_l.html
[~ is a tilde]
Arnold Neumaier
===
Subject: Finding angles using program
How do we find angles in (-180 deg, +180 deg] domain?
(If we use dot product, then +45 deg and -45 deg can not be
distinguished)
Are there any particular functions in Matlab using which I
can find
the angles in the (-180 deg, +180 deg] domain?
Chetan
===
Subject: Information about Wavelet Analysis
I am looking for information (book name internet source)
about Wavelet
analysis.
===
Subject: Re: Information about Wavelet Analysis
S. Mallat, A Wavelet Tour of Signal Processing, Academic Press
is said being the bible of wavelet analysis.
Axel
===
Subject: Re: Finding angles using program
>How do we find angles in (-180 deg, +180 deg] domain?
>(If we use dot product, then +45 deg and -45 deg can not be
>distinguished)
>Are there any particular functions in Matlab using which I
can find
>the angles in the (-180 deg, +180 deg] domain?
>Chetan
It depends. In 2d there is a simple method: for vectors
a=(a1,a2) and
b=(b1,b2), do the dot product first, then calculate the
determinant
a1b2 - a2b1. If the result is >0 then a is clockwise relative
to b,
if < 0 then a is ccw from b (and if =0 then a and b are
aligned but
then the angle will be 0 anyway). You then have to adjust
your angle
accordingly.
(you will notice that this determinant is really the z
component of
the 3d cross-product).
In 3d I think you would have to arbitrarily define an
orientation
relative to the plane that contains the 2 vectors.
===
Subject: Re: Solving Nonlinear Equation
>Hi
> I am chemical engineer and I am solving some nonlinear
equations
> F = F(x)
> I was using some solvers for nonlinear equation. But as
these
>solvers calculate numerical Jacobian, it perturbs each
variable and
>calculates jacobian. This makes it very time consuming.
Thats why I
>used approximate jacobian (It is standard technique in many
chemical
>problems) and used Newton raphson method to solve the
equations.
> Now my problem is :
> Many time I dont have very good initial guess and it is
sometimes
>far away from the solution. So when i use Jacobian (either
approximate
>or numerically calculated exact one), and solve for
> J * delX = -F
> for delX, the step delX is quite big.
> If I use X(k+1) = X(k) + delX
> then every time the soution goes somewhere and it
fails..NEVER
>converges.
> So i use some sort of relaxation factor
> X(k+1) = X(k) + w * delX where w = 0 to 1
> So in the beginning iterations, i use small value of w (say
0.1 o
>0.2), and in later sateg increase it towarss 0.8, 0.9 when
it is close
>to convergence criteria.
> This works and i get solution. But value of w is manual. I
am not
>using any standard techniqus to get it. So is there any
justification
>for this factor? is there any way (fast and reliable) to find
this
>factor?
>
>
dont reenvent the wheel. what you are doing is the damped
Newton method
with a handcrafted damping. there is a complete and sound
theory for this
developed in the seventies of the last century and even
downloadable code,
the nleq-series of codes from the codelib/elib.
(deußhards group). the essential idea is to choose w (in your
notation)
such
that
norm( F(X(k)+w*delx(k)) )^2 <= (1-factor*delx)*norm(F(X(k)))^2
where
0 < factor < 1/2 for example factor =0.01
and norm(.)^2 is a scalar product, for example norm=euclidean
length
or norm(F)=euclidean length of (A*F) with a fixed matrix A
for example A=approximate Jacobian. One can show that this
converges to a
solution
if on the set of xs with F-values norm(F(x)) <=
norm(F(x(0))) the Jacobian
becomes never singular.
for codes see
http://plato.la.asu.edu/topics/problems/zero.html
hth
peter
===
Subject: Re: Computational Evolution without velocities
>
>
I have heard the folklore that when one wants to simulate some
structure
> Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
> update formula
x_new = x_old - Force DeltaTime / 2 Mass
holds for any simulation and we do not need to take care of
velocities.
>
> This equation is quite obviously wrong, for the simple fact
that it
> does not even have consistent units. The left hand side has
dimensions
> of [Length]; the right hand side is a length minus a
velocity. You cant
> subtract a velocity from a length and get anything
meaningful.
>
> -- Gordon D. Pusch
>
He obviously meant
x_new = x_old - Force * (DeltaTime)^2 / (2 * Mass)
However, that does not seem right. Newtons second law is
x = F/m
If you wish to eliminate velocities, then in an obvious
notation,
x(+) = 2x(0) - x(-) + (F/m) * (delta_t)^2
will step forward. But you then have to have the values x(0)
and x(-).
I dont know what algorithm the OP could have been talking
about since
doesnt make sense physically.
--
Julian V. Noble
Professor Emeritus of Physics
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
Science knows only one commandment: contribute to science.
-- Bertolt Brecht, Galileo.
===
Subject: Re: Need help with steepest descent
>
> Would someone care to explain how to use the steepest
descent
> algorithm to minimize a function, say, F(x) = x1 * x1 + 25
* x2 * x2,
> -5 < x1,x2 < 5 :-). Links to web tutorials would also be
appreciated.
> And no, this is not my homework, my homework is to minimize
> Rosenbrocks function :-). (I, uh, fell aspeep during the
class :-(
> and now I cant find any references on the
net.
>
>
> Andrey
check out mathworld.
http://mathworld.wolfram.com/MethodofSteepestDescent.html
hope you have a mathematica. if not, matlab.
http://www-math.cudenver.edu/~aknyazev/teaching/98/4660/fp/r1/
the idea is to perturb a parameter by an increment i.( rather
arbitrary amount, you get to decide) by perturb, I mean add or
subtract. it doesnt matter whether you add
first or subtract
first.
lets add incrememnt i to the x1 and x2, then look at what
happens to
the F[x1, x2] that is look at it;s numerical derivative.
so dfx1 = (f[x1+i,x2] - f[x1, x2])/ i and dfx2 = (f[x1,x2+i]
- f[x1,
x2])/ i
these are numerical derivatives. what happens if the
derivative is
negative? this means you are moving downhill towards the
local minima.
you want to keepp this move.
if this is positive than youre moving uphill then you want
to make
the opposite move.
ßank the routine with while loop with conditions to stop it
when it
reaches a certain minima criteria.
good luck
sean
good luck
===
Subject: Re: Bernstein-B.8ezier coefficients?
> I have to implement a specific algorithm.
> Just one question: if one talks about a Bernstein-B.8ezier
coefficient
> of a tensor-product surface one means a control point with
x,y,z
> -coordinates, right?
I think it means v = fn(x,y,z)
===
Subject: Re: analysis of a large ODE system. as per Dr.
Sspelluccis
suggestion.
to someone in the forum and perhaps in my school.
I went to the library and checked out the book Dr Spellucci
has
suggested.( hahns stability of motion in conjunction with
the
jacobians and eigenvalues) and guess what?
its been recalled by another person in my school. I wonder
who else
from my school is reading these forums... I think Ill
recall
it when
I return it.
what are the odds of two people from same school looking for
the same
book on the same week right after the book was suggested in a
newsgroup.
specially when the suggestions for it came after eigenvalues
comments.
perhaps the book contains some information not found anywhere
in the
literature. ( cant imagine what that will be) Though, if
you
need the
book like me, then you cant be all that great. LOL
-----
To Dr. Spellucci,
thats a great book btw, Dr. Spelluci. (its
very
mathematically
formalized.) perhaps thats why you have suggested it.
I have acces,now to some nice introductory materials into
nonlinear
dynamics. Im sure I would learn a lot based on your
comments
and
those of others. It really gives me what to focus on when
reading the
materials.
I really appreciated your comments as well as others. You
have always
been so helpful. I wish people in academics were more like
you. (I
sound bitter. LOL)
respectfully yours,
sean from UCIrvine
===
Subject: Re: Wavelets package in Java
>
> Has anyone ever used a wavelet (de)composition package
under Java? Or
> a program that can be used with Java (looked at Jython, but
I didnt
> find a wavelet package for that)?
>
> Please mail replies also to gez_75@hotmail.com
>
>
> Gez
This looks promising:
http://www.bearcave.com/software/java/wavelets/. Also
http://www.bearcave.com/misl/misl_tech/wavelets/index.html
(same site).
===
Subject: Re: question about 3d plottin software data table
limit
> Hy, I have a quite noise problem, by a wavelet analysis
Ive obtained
> a data table of frequency, time and energy, the problem is
that the
> table is very long (3 columns and somthing like 13 milion
of rows) I
> need to plot this table in a 3d graphic, but at now I
dont
find a
> program (quite simple to use) that has a so high limit in
is data
> table lenght (usually are of 2^15 or 2^16 of rows). Do you
know is
> exist a program with a no limit data table lenght? Also a
shareware or
> trial version is good to me.
> Kio
Just remark that a standart computer display has between
1 and 1.5 millions of pixels. I think its hopeless to try
to show 13 millions of informations on it.
you can try xd3d
http://www.cmap.polytechnique.fr/~jouve/xd3d/
It will do the job :)
--
F.J.
===
Subject: A Least Squares Problem
Im trying to solve a least-squares problem:
min{(Ax+b)(Ax+b)} (1)
whose explicit solution is:
x = (A*A)^-1 * A * b (2)
Matrix A has form, e.g.:
1 0 0 0 1
-1 1 0 1 0
0 -1 1 -1 0
0 0 -1 0 -1
i.e. it has exactly one pair of (1,-1) in each column. There
wont be
any all-zero rows. As you can see, it may contain two or more
linearly
dependent columns (two equal, two mutually negative, or other
combination).
Further, due to nature of the problem (Id skip the physical
background
for now, I can repost if someones interested) condition is
that
size(x) >= size(b)-1
The problem is that equation (2) has explicit solution (A*A
has
inversion) only if
size(x) = size(b)-1 = rank(A)
However, I have to solve the problem in general manner, i.e.
I want
to obtain one of optimal solutions. Rank deficiency of A
should be
preferredly solved by minimizing square sum of x, i.e.
min{|x^2|},
i.e. the task is to find such x that a) minimizes (1) and b)
has
minimal sum(x*x).
I plan to use Matlab as proving environment and
Fortran/Lapack for
deployment (links to Fortran sources which solve that exact
problem are
also welcome). I probably have to find QR decomposition of
AA
(GEQRF),
but Im lost at the moment what to do after that. I admit my
linear
algebra is a little rusty...
Theres an interesting example under Matlabs
help on QR
function,
but its not directly applicable -- in my case,
(AA) is (can
be)
singular.
All suggestions/solutions are welcome. I suspect this is a
well-known
--
Jugoslav
___________
www.geocities.com/jdujic
===
Subject: Re: A Least Squares Problem
An nxp matrix A can be factored into a product of 3 other
matrices :
A=U*L*Vt (where t==transpose) where n for n data and p for p
parameters,
using singular value decomposition SVD. U and V are data and
parameter
space
eigenvectors and L a diagonal matrix containing at most r
non-zero
eigenvalues of G with r<=p. So if
x = (A*A)^-1 * A * b
At*A=V*L^2*Ut*U*L*Vt=V*L^2*Vt ,since Ut*U=I
and (At*A)^-1=V*L^-1*Ut and finally, since Vt*V=I
The least square solution is given by
x=V*L^-1*Ut*b
SVDCMP.f from numerical receips in fortran, and svd(A) from
matlab give
similar results.
SVD is the best way to decompose a matrix since it provides
extra
infomation
about the resolution of the reconstruction (resolution
matrices etc...)
--
---FiLiP Louis---
> please dont redirect as I dont read it
regularly Im
trying to solve a least-squares problem:
> min{(Ax+b)(Ax+b)} (1)
> whose explicit solution is:
> x = (A*A)^-1 * A * b (2)
> Matrix A has form, e.g.:
> 1 0 0 0 1
> -1 1 0 1 0
> 0 -1 1 -1 0
> 0 0 -1 0 -1
> i.e. it has exactly one pair of (1,-1) in each column.
There wont be
> any all-zero rows. As you can see, it may contain two or
more linearly
> dependent columns (two equal, two mutually negative, or
other
combination).
> Further, due to nature of the problem (Id skip the
physical background
> for now, I can repost if someones interested) condition
is
that
> size(x) >= size(b)-1
> The problem is that equation (2) has explicit solution
(A*A has
> inversion) only if
> size(x) = size(b)-1 = rank(A)
> However, I have to solve the problem in general manner,
i.e. I want
> to obtain one of optimal solutions. Rank deficiency of A
should be
> preferredly solved by minimizing square sum of x, i.e.
min{|x^2|},
> i.e. the task is to find such x that a) minimizes (1) and b)
has
> minimal sum(x*x).
> I plan to use Matlab as proving environment and
Fortran/Lapack for
> deployment (links to Fortran sources which solve that exact
problem are
> also welcome). I probably have to find QR decomposition of
AA (GEQRF),
> but Im lost at the moment what to do after that. I admit
my linear
> algebra is a little rusty...
> Theres an interesting example under
Matlabs help on QR
function,
> but its not directly applicable -- in my case,
(AA) is
(can be)
> singular.
> All suggestions/solutions are welcome. I suspect this is a
well-known
> --
> Jugoslav
> ___________
> www.geocities.com/jdujic
===
Subject: Re: A Least Squares Problem
|| Im trying to solve a least-squares problem:
||
|| min{(Ax+b)(Ax+b)} (1)
||
|| whose explicit solution is:
||
|| x = (A*A)^-1 * A * b (2)
||
|| The problem is that equation (2) has explicit solution
(A*A has
|| inversion) only if
|| size(x) = size(b)-1 = rank(A)
||
|| However, I have to solve the problem in general manner,
| An nxp matrix A can be factored into a product of 3 other
matrices :
| A=U*L*Vt (where t==transpose) where n for n data and p for
p parameters,
| using singular value decomposition SVD. U and V are data
and parameter
space
| eigenvectors and L a diagonal matrix containing at most r
non-zero
| eigenvalues of G with r<=p. So if
| x = (A*A)^-1 * A * b
| At*A=V*L^2*Ut*U*L*Vt=V*L^2*Vt ,since Ut*U=I
| and (At*A)^-1=V*L^-1*Ut and finally, since Vt*V=I
| The least square solution is given by
| x=V*L^-1*Ut*b
...
| SVD is the best way to decompose a matrix since it provides
extra
infomation
| about the resolution of the reconstruction (resolution
matrices etc...)
if r
|
| Im trying to solve a least-squares problem:
|
| min{(Ax+b)(Ax+b)} (1)
|
| whose explicit solution is:
|
| x = (A*A)^-1 * A * b (2)
| The problem is that equation (2) has explicit solution
(A*A has
| inversion) only if
| size(x) = size(b)-1 = rank(A)
|
| However, I have to solve the problem in general manner,
i.e. I want
| to obtain one of optimal solutions. Rank deficiency of A
should be
| preferredly solved by minimizing square sum of x, i.e.
min{|x^2|},
| i.e. the task is to find such x that a) minimizes (1) and b)
has
| minimal sum(x*x).
I appear to have found a practical solution. First, sorry,
equation (2)
should have read:
x = A * (A*A)^-1 * b
however, I get excellent results using instead
x = A * (A*A + e*I)^-1 * b
where e is a small constant (e.g. 0.001).
I should have mentioned that a very accurate result is not
required
(entire stuff is part of a larger algorithm which is
semi-discrete
and heuristic, so something like this is perfectly useful).
Id like
to know, however, what I am doing with this -- am I solving
something
like
min{(Ax+b)(Ax+b) + e||x||^2}
?
--
Jugoslav
___________
www.geocities.com/jdujic
===
Subject: Sinc Interpolation Question
I have a vector y with 8 elements and a corresponding 8 x 8
(diagonal)
covariance matrix Sy. I want to perform sinc interpolation on
y and
compute the corresponding covariance.
To sinc interpolate I perform:
y = [224.2944 224.4289 230.3499 239.9113 251.0483 259.5003
260.7612 250.3568];
F = fft(eye(8));
F2 = fft(eye(16));
zf = [0 0 0 0];
y_interp = F2 * circshift([zf; circshift(1/8 *
F * y, 4);
zf],8);
In the above I can replace the circshift function fftshift to
make it
easier to understand.
How do I go about computing the corrected covariance
Se_interp?
Should it be positive-definite?
===
Subject: Re: A Least Squares Problem
>| | please dont redirect as I dont read it
regularly>
Actually, sci.stat.math or sci.stat.consultwould be just as
good.
>Id like
>to know, however, what I am doing with this -- am I solving
something
>like
>min{(Ax+b)(Ax+b) + e||x||^2}
Find some references on ridge regression (dont have any off
the top of my
head, and my reference books are several miles away).
A look at Raos book (Linear Statistical Inference, Wiley)
might help you
in understanding the rank deficient case.
Other possibilities:
1) Use the SVD to determine the inverse. The usual deal there
is to assume
that the singular values near zero are equal to zero. In
exact arithmetic,
this procedure leads to the minimum L2 norm solution you are
looking for.
2) If A is large, try conjugate gradients. It can do
amazingly well on
problems of this sort.
--
My real email address is
mcintosh ##at## research ##dot## telcordia ##dot## com
===
Subject: Multivariate Adaptive Regression Splines (MARS)
Hi All,
I am using MARS 3.5 for scattered data points approximation
2D, 3D, 5D.
What
minimum number of data points is required to start build
approximation in
MARS?
--
rm -r /
Virginia Tech
===
Subject: Re: A Least Squares Problem
>please dont redirect as I dont read it
regularly>
>
>Im trying to solve a least-squares problem:
>
>min{(Ax+b)(Ax+b)} (1)
so Ax+b is a row for you? then: what is Ax ?
>
>whose explicit solution is:
>
>x = (A*A)^-1 * A * b (2)
>
>Matrix A has form, e.g.:
> 1 0 0 0 1
> -1 1 0 1 0
> 0 -1 1 -1 0
> 0 0 -1 0 -1
if this is A then x should have 5 and b 4 components
or otherwise your notation is faulty
for this A A*A is 5 times 5 and of rank 4 hence no inverse
exists.
snip
obviously you are interested in the minimum norm least
squares solution
of
norm(Ax+b,2) = min_x
the best way to compute this is via the svd (available in
matlab and
also in lapack)
it first computes (in the full -non economy version )
U , S and V which yield
A=USV
here U is unitary and has as much rows as A,
UU=UU=I
V is also unitary and has as much columns as A, abd
VV=VV=I
(so U and V are quadratic) and S has the same shape as A, but
has only
some
nonnegative elements along the diagonal, all other being zero
.
these elements are the singular values of A and their number
is the
rank of A.
now you compute the minimum norm least squares solution x as
x=V*S# *U*b
where
S# has the same shape as S (transpose of S) and along the
diagonal
the reciprocals of the nonzero singular values of A at the
positions
they
occur in S:
formally
(S#)(i,i) = 1/S(i,i) if S(i,i) not= 0, and =0 otherwise.
if you compute under roundoff, you must replace the decision
S(i,i) not=0 by S(i,i)>alpha>0 for some reasonable alpha
e.g.
alpha=100*eps*norm(A) in matlab notation.
you can get all his also from LAPACK and its derivatives (in
other
languages).
hth
peter
===
Subject: help needed for runaway ODE?
Hi All,
I am solving a robotic control problem, envolving coupled
nonlinear ODEs. I
found the solution to the ODEs using relaxation method from
the Numeric
Receipe in C for the two-point boundary value problem. But
after I changed
the relaxation method to Runge-Kutta method, using the
initial values
solved
by relaxation method, I got a runaway type of solution with
all solution
grow to very large numbers. I also tried the stiff solver in
the NRC, the
results were similar. I know my derivatives are correct in
Range-Kutta
because I computed derivative from the solution of relaxation
method and
they matched the derivatives from my Range-Kutta routine. But
the
derivatives were wrong when I computed it from the solution
of the
Range-Kutta method. Any suggestions from numeric gurus here?
I have been
struggling on this simple ODE problem for three weeks.
Everett
send your comments to everteq AT sbcglobal DOT net or post it
here. Many
===
Subject: Re: A Least Squares Problem
> | An nxp matrix A can be factored into a product of 3 other
matrices :
> | A=U*L*Vt (where t==transpose) where n for n data and p
for p
parameters,
> | using singular value decomposition SVD. U and V are data
and parameter
space
> | eigenvectors and L a diagonal matrix containing at most r
non-zero
> | eigenvalues of G with r<=p. So if
> | x = (A*A)^-1 * A * b
> | At*A=V*L^2*Ut*U*L*Vt=V*L^2*Vt ,since Ut*U=I
> | and (At*A)^-1=V*L^-1*Ut and finally, since Vt*V=I
> | The least square solution is given by
> | x=V*L^-1*Ut*b
> ...
> | SVD is the best way to decompose a matrix since it
provides extra
infomation
> | about the resolution of the reconstruction (resolution
matrices etc...)
> if r Remind you, the additional criterion for resolving rank
deficiency should
> be min{||x||}.
There is a little mistake in the solution above: If you don
not invert
L itself, which is not possible in the rank deficient case,
but invert only
all nonzero elements of L, you get the uniqe x with minimum
2-norm.
--
Dr. rer. nat. Uwe Schmitt http://www.procoders.net
schmitt@procoders.net A service to open source is a service to
mankind.
===
Subject: Re: help needed for runaway ODE?
> Hi All,
> I am solving a robotic control problem, envolving coupled
nonlinear ODEs.
I
> found the solution to the ODEs using relaxation method from
the Numeric
> Receipe in C for the two-point boundary value problem. But
after I
changed
> the relaxation method to Runge-Kutta method, using the
initial values
solved
> by relaxation method, I got a runaway type of solution with
all solution
> grow to very large numbers. I also tried the stiff solver
in the NRC, the
> results were similar. I know my derivatives are correct in
Range-Kutta
> because I computed derivative from the solution of
relaxation method and
> they matched the derivatives from my Range-Kutta routine.
But the
> derivatives were wrong when I computed it from the solution
of the
> Range-Kutta method. Any suggestions from numeric gurus
here? I have been
> struggling on this simple ODE problem for three weeks.
Its hard to make a substantive suggestion without seeing
the
exact
equations, but it sounds like a job for implicit or
semi-implicit
methods.
===
Subject: Re: =?iso-8859-1?Q?Bernstein=2DB=E9zier?=
coefficients?
boundary=------------2CF45486D96217C5826AF13D
-------------------------------------------------------------
--------
> I have to implement a specific algorithm.
> Just one question: if one talks about a Bernstein-B.8ezier
coefficient
> of a tensor-product surface one means a control point with
x,y,z
> -coordinates, right?
I dont know what specific classification of a
tensor product
surface you
refer to in your mail : a Generic B-Spline or a NURBS or a
Bezier
Surface (since you seem to
refer to Bernstein polynomials and Bezier Basis Functions in
your mail)
, and therefore,
Starting with a generic definition, of a Tensor-Product
Surface,
S(u,v) = B i* Pi(x,y,z)
where Bi is a generic Coefficient (a product of two curve
coefficients
for a surface)
and therefore:
[a] The Coefficient represents a term (of the polynomial
function)
that maps a point in R(3) space to
R(2) (or parametric) space. The point is the control point of
the
Tensor Product Surface for the specific u,v
value in the parametric space where u -> [a,b] and v->[c,d]
[b] A Point on the tensor-product surface and a Control Point
of the
tensor-product surface
are separate concepts and dont mean the same thing.
hope this solves the problem.
Anup R. Joshi
www.me.mtu.edu/~arjoshi
===
Subject: Re: A Least Squares Problem
Google that reference.
[38] S. M. Tan and C. Fox, Inverse Problems 453.707
classnotes, The
University of Auckland, Aucland, New Zealand., 2001.
What you are doing is referred to as regularized solution, or
damped
least
squares.
Read the first 3 chapters of the above reference, excellent
and very brief
reading on what you need to know when performing least
squares, or damped
(regularized) least squares.
Alien+
> An nxp matrix A can be factored into a product of 3 other
matrices :
> A=U*L*Vt (where t==transpose) where n for n data and p for
p parameters,
> using singular value decomposition SVD. U and V are data
and parameter
space
> eigenvectors and L a diagonal matrix containing at most r
non-zero
> eigenvalues of G with r<=p. So if
> x = (A*A)^-1 * A * b
> At*A=V*L^2*Ut*U*L*Vt=V*L^2*Vt ,since Ut*U=I
> and (At*A)^-1=V*L^-1*Ut and finally, since Vt*V=I
> The least square solution is given by
> x=V*L^-1*Ut*b
> SVDCMP.f from numerical receips in fortran, and svd(A) from
matlab give
> similar results.
> SVD is the best way to decompose a matrix since it provides
extra
infomation
> about the resolution of the reconstruction (resolution
matrices etc...)
> --
> ---FiLiP Louis---
> please dont redirect as I dont read it
regularly
> Im trying to solve a least-squares problem:
min{(Ax+b)(Ax+b)} (1)
whose explicit solution is:
x = (A*A)^-1 * A * b (2)
Matrix A has form, e.g.:
> 1 0 0 0 1
> -1 1 0 1 0
> 0 -1 1 -1 0
> 0 0 -1 0 -1
> i.e. it has exactly one pair of (1,-1) in each column.
There wont be
> any all-zero rows. As you can see, it may contain two or
more linearly
> dependent columns (two equal, two mutually negative, or
other
> combination).
> Further, due to nature of the problem (Id skip the
physical background
> for now, I can repost if someones interested) condition
is
that
> size(x) >= size(b)-1
The problem is that equation (2) has explicit solution (A*A
has
> inversion) only if
> size(x) = size(b)-1 = rank(A)
However, I have to solve the problem in general manner, i.e.
I want
> to obtain one of optimal solutions. Rank deficiency of A
should be
> preferredly solved by minimizing square sum of x, i.e.
min{|x^2|},
> i.e. the task is to find such x that a) minimizes (1) and b)
has
> minimal sum(x*x).
I plan to use Matlab as proving environment and
Fortran/Lapack for
> deployment (links to Fortran sources which solve that exact
problem are
> also welcome). I probably have to find QR decomposition of
AA (GEQRF),
> but Im lost at the moment what to do after that. I admit
my linear
> algebra is a little rusty...
Theres an interesting example under Matlabs
help on QR
function,
> but its not directly applicable -- in my case,
(AA) is
(can be)
> singular.
All suggestions/solutions are welcome. I suspect this is a
well-known
--
> Jugoslav
> ___________
> www.geocities.com/jdujic
===
Subject: Re: help needed for runaway ODE?
Hi Everett,
try my ODEs from http://www.delphipages.com/result.cfm?ID=3482
> Hi All,
>
> I am solving a robotic control problem, envolving coupled
nonlinear ODEs.
I
> found the solution to the ODEs using relaxation method from
the Numeric
> Receipe in C for the two-point boundary value problem. But
after I
changed
> the relaxation method to Runge-Kutta method, using the
initial values
solved
> by relaxation method, I got a runaway type of solution with
all solution
> grow to very large numbers. I also tried the stiff solver
in the NRC, the
> results were similar. I know my derivatives are correct in
Range-Kutta
> because I computed derivative from the solution of
relaxation method and
> they matched the derivatives from my Range-Kutta routine.
But the
> derivatives were wrong when I computed it from the solution
of the
> Range-Kutta method. Any suggestions from numeric gurus
here? I have been
> struggling on this simple ODE problem for three weeks.
>
>
> Everett
>
> send your comments to everteq AT sbcglobal DOT net or post
it here. Many
===
Subject: Discrete Chebyshev (or other orthogonal) polynomials.
For my research I want to find out more about discete
Chebyshev
polynomials of the first kind and their relation to Chebyshev
polynomials of the second kind.
In literature a lot of properties of the Chebyshev
polynomials and the
relations to other polynomials are known. But, afaik, all
these nice
propertie are derived under the assumption that we want to
approach a
continue function. I cannot find a list of properties in case
we want to
approach a discrete set of equidistant spaced points with
these
polynomials.
The description in Numerical Recipes for example also assumes
the the
function is known at all locations.
interesting words but nothing came out of it.
===
Subject: Help with adjustment to quadratic regression curve
algorithm
Im trying to understand/get an algorithm to work which
transforms a
quadratic regression curve to be unimodal using Bezier
control points. The
algorithm works out the transform.
Please understand that I am inexperienced in this field.
I can post the whole algorithm if its needed to put things in
context (its
not very long), but what Im particularly stuck on is this
equation:
m
Y0:=1/m SIGMA (Pi-(Y1)Si = (Y2)(Si)^2 )
i=1
where: m=number of points (Si,Pi),
Y1,Y2 are already transformed Ôouter points.
What does the = mean in (Pi-(Y1)Si = (Y2)(Si)^2 ) ?
The text states that:
this statement evaluates the best least squares estimate of
Y0 for the
parabola with new values of Y2 and Y1 resulting from the first
adjustment.
The second adjustment, a vertical translation of the
B.8ezier-reshaped
parabola, occurs only when and Y1 and Y2 have been altered.
In fact, Im not sure that the steps in the algorithm
Ive
done before this
equation are right either:(
So ANY help appreciated.
PS I have access to matlab
TIA
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
HI,
did u take a look at - Conte/DeBoor: Elementary Numerical
Analysis, An
Algorithmic Approach- ?
I found it very useful, many times...
Ivano
Maurice ha scritto nel
messaggio
> For my research I want to find out more about discete
Chebyshev
> polynomials of the first kind and their relation to
Chebyshev
> polynomials of the second kind.
> In literature a lot of properties of the Chebyshev
polynomials and the
> relations to other polynomials are known. But, afaik, all
these nice
> propertie are derived under the assumption that we want to
approach a
> continue function. I cannot find a list of properties in
case we want to
> approach a discrete set of equidistant spaced points with
these
polynomials.
> The description in Numerical Recipes for example also
assumes the the
> function is known at all locations.
> interesting words but nothing came out of it.
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
>
>
> For my research I want to find out more about discete
Chebyshev
> polynomials of the first kind and their relation to
Chebyshev
> polynomials of the second kind.
> In literature a lot of properties of the Chebyshev
polynomials and the
> relations to other polynomials are known. But, afaik, all
these nice
> propertie are derived under the assumption that we want to
approach a
> continue function. I cannot find a list of properties in
case we want to
> approach a discrete set of equidistant spaced points with
these
polynomials.
> The description in Numerical Recipes for example also
assumes the the
> function is known at all locations.
> interesting words but nothing came out of it.
>
Let me give you three references:
1. My class notes for a numerical analysis course: look on
http://www.phys.virginia.edu/classes/551.jvn.fall01/
You want Chapter 2, Representation of Functions.
2. Abramowitz & Stegun, Handbook of Mathematical Functions
(Dover PB)
3. Ralston, Introduction to Numerical Analysis
--
Julian V. Noble
Professor Emeritus of Physics
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
Science knows only one commandment: contribute to science.
-- Bertolt Brecht, Galileo.
===
Subject: Books On Generalized Reduced Gradient Algorithm?
Can any one recommend books or papers I can read that
describe how to
implement a Generalized Reduced Gradient Algorithm (such as
that used in
Excels Solver tool)? Source code implementing this
algorithm
would be
very useful as well.
-Vik
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
>
>For my research I want to find out more about discete
Chebyshev
>polynomials of the first kind and their relation to Chebyshev
>polynomials of the second kind.
>In literature a lot of properties of the Chebyshev
polynomials and the
>relations to other polynomials are known. But, afaik, all
these nice
>propertie are derived under the assumption that we want to
approach a
>continue function. I cannot find a list of properties in case
we want to
>approach a discrete set of equidistant spaced points with
these
polynomials.
>The description in Numerical Recipes for example also
assumes the the
>function is known at all locations.
>interesting words but nothing came out of it.
>
>
discrete orthogonal for what scalar product? the chebyshev
polynomials of
degree <= m of the first kind are orthogonal
on the discrete set {x_k} = zeros of T_{m+1} with the
ordinary scalar
product.
for more information see
Chebyshev Polynomials
by J.C. Mason (University of Huddersfield, UK)
and D.C. Handscomb (lately of Oxford University)
Chapman & Hall / CRC Press, 2002, xiii+341 pp.
Hardback ISBN 0-8493-0355-9, $99.95
or
Pure and Applied Mathematics. New York: John Wiley & Sons,
Inc. xvi, 249
p. (1990)
hth
peter
===
Subject: Re: Books On Generalized Reduced Gradient Algorithm?
I believe this is actually Leon Lasdons GRG2 algorithm.
Some references are:
@misc{ fylstra98design,
author = D. Fylstra and L. Lasdon and A. Warren and J. Watson,
title = Design and use of the microsoft excel solvers,
text = D. Fylstra, L. Lasdon, A. Warren, and J. Watson,
Design and use
of the
microsoft excel solvers. To appear in Interfaces, 1998.,
year = 1998,
url = citeseer.nj.nec.com/fylstra98design.html }
L. S. Lasdon and A. D. Waren. Generalized reduced gradient
software for
linearly and nonlinearly constrained problems. In H. J.
Greenberg, editor,
Design and Implemetation of Optimization Software, pages
335---362. Sijthoff
and Noordhoff, Netherlands, 1978.
Lasdon L. S., Warren A. D., Jain A., and Ratner M. 1978.
Design and Testing
of a GRG Code for Nonlinear Optimization, ACM Trans.
Mathematical Software,
4: 3450.
-------------------------------------------------------------
---
Erwin Kalvelagen
erwin@gams.com, http://www.gams.com/~erwin
-------------------------------------------------------------
---
> Can any one recommend books or papers I can read that
describe how to
> implement a Generalized Reduced Gradient Algorithm (such as
that used in
> Excels Solver tool)? Source code implementing this
algorithm would be
> very useful as well.
> -Vik
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
> 2. Abramowitz & Stegun, Handbook of Mathematical Functions
> (Dover PB)
Or download it for free.
http://jove.prohosting.com/~skripty/
===
Subject: skyline matrix solver
Please,
I need a skyline matrix solver algorithm. Where can I find
one?
Cpplayer
===
Subject: Re: Books On Generalized Reduced Gradient Algorithm?
>Can any one recommend books or papers I can read that
describe how to
>implement a Generalized Reduced Gradient Algorithm (such as
that used in
>Excels Solver tool)? Source code implementing this
algorithm would be
>very useful as well.
>
>
>-Vik
the textbook by
Avriel & Golany : mathematical programming for industrial
engineers
(marcel Decker publisher) contains a chapter on nonlinear
programming
written
by Lasdon and coworkers which gives very much details about
this.
source code is not on the net. this is in commercial (and
rather
expensive)
products.
hth
peter
===
Subject: Eignevector differentiation
Is it possible to differentiate the eigenvalues and
eigenvectors of a
square
symmetric matrix A with respect to any component of A?
I.e. if A has ith eigenvalue lambda_i and
ith eigenvector
v_i, what is :
d(lambda_i)/d(A_mn)
and what is :
d(v_i)/d(A_mn)
Mike Fairbank.
===
Subject: Re: Eignevector differentiation
> Is it possible to differentiate the eigenvalues and
eigenvectors of a
square
> symmetric matrix A with respect to any component of A?
>
> I.e. if A has ith eigenvalue lambda_i and
ith eigenvector
v_i, what is
:
>
> d(lambda_i)/d(A_mn)
>
> and what is :
>
> d(v_i)/d(A_mn)
Basically I think you can differentiate
the eigenvalues, provided they are not multiple. If they are
multiple,
things get
tricky.
The Implicit Value Theorem gives you the derivative for
non-multiple
eigenvalues,
if you plug the characteristic polynomial of the matrix.
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
Distribution: inet
>
>>2. Abramowitz & Stegun, Handbook of Mathematical Functions
>> (Dover PB)
>
> Or download it for free.
> http://jove.prohosting.com/~skripty/
Any evidence that this book is out of copyright, and that
its now okay
to post it on the web?
Andrew
--
Andrew Ross
Industrial and Systems Engineering
Lehigh University, www.lehigh.edu/~inime/
Bethlehem, Pennsylvania, USA
remove all digits <=4 from my e-mail address to reply
===
Subject: Re: Eignevector differentiation
>
> Is it possible to differentiate the eigenvalues and
eigenvectors of a
square
> symmetric matrix A with respect to any component of A?
Yes, if (and only if) the eigenvalue is simple. But the
normalization
condition for the eigenvector must be specified.
Let A=A(t) depend on some parameter t (e.g. t=A_ik).
Then locally (assuming the eigenvalue is simple) there is a
solution of
A(t)x(t)=lam(t)x(t), x(t)^Tx(t)=1
by the implicit function theorem. Differentiate to get
A(t)x(t)+A(t)x(t)=lam(t)x(t)
+lam(t)x(t), x(t)^Tx(t)=0.
Rewrite this as a linear system for x(t) and
lam(t) an
solve by
Gaussian elimination. It is not difficult to prove that the
system
is nonsingular for a simple eigenvalue. The same approach
works for
other normalizations such as e^Tx(t)=1, where e is a unit
vector.
Arnold Neumaier
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
>
>2. Abramowitz & Stegun, Handbook of Mathematical Functions
>> (Dover PB)
Or download it for free.
> http://jove.prohosting.com/~skripty/
>
> Any evidence that this book is out of copyright, and that
its now okay
> to post it on the web?
The Dover version might be copyright, but the hardback is a US
government publication, and hence (IIRC) not
copyright-protected.
Phil Hobbs
===
Subject: Re: Eignevector differentiation
Would that give me the derivatives of the eigenvectors as
well?
BTW I think I might have missed some relevant background to
my problem.
This may or may not make the solution to the above problem
clearer. Here
it
is:
I am investigating the least squares solution to the matrix
equation Ax=B.
This least squares solution is x=(AA)^-1AB.
I want to
differentiate x
with respect to A, i.e. find dx/dA_mn, but when
AA is close
to singular.
For this I need to know how (AA)^-1 changes with respect to
A. I already
have this answer when AA is not close to singular: then you
can just use
d((AA)^-1)/dA_mn
=-(AA)^-1.(d(AA)/d(A_mn)).(AA
)^-1
However as AA it is close to singular I should use the
Singular Value
Decompostion algorithm to invert AA. Instead of SVD, for
efficiency, I am
using an eigenvector finding algorithm since this gives the
same thing, I
think (i.e. when SVD algorithm is applied to
AA=UWV then W
always turns
out to be the diagonal matrix of eigenvalues of AA, and U=V
is the square
matrix of all the eigenvectors of AA). Thats
why I want to
differentiate
the eigenvectors and eigenvalues.
So once you have the SVD, you can invert AA easily as
VW^-1U while
neglecting any terms in W that were close to zero. So
shouldnt it be
possible to differentiate this SVD inversion, since it was
possible when
AA
wasnt nearly singular, and the SVD inverion algorithm is
clearly defined.
Hope that makes things more clear!
Mike Fairbank
> Is it possible to differentiate the eigenvalues and
eigenvectors of a
square
> symmetric matrix A with respect to any component of A?
I.e. if A has ith eigenvalue lambda_i and
ith eigenvector
v_i, what
is
:
d(lambda_i)/d(A_mn)
and what is :
d(v_i)/d(A_mn)
> Basically I think you can differentiate
> the eigenvalues, provided they are not multiple. If they
are multiple,
things get
> tricky.
> The Implicit Value Theorem gives you the derivative for
non-multiple
eigenvalues,
> if you plug the characteristic polynomial of the matrix.
===
Subject: Re: Eignevector differentiation
> Is it possible to differentiate the eigenvalues and
eigenvectors of a
square
> symmetric matrix A with respect to any component of A?
I.e. if A has ith eigenvalue lambda_i and
ith eigenvector
v_i, what is
:
d(lambda_i)/d(A_mn)
and what is :
d(v_i)/d(A_mn)
> Basically I think you can differentiate
> the eigenvalues, provided they are not multiple. If they are
multiple, things get tricky.
> The Implicit Value Theorem gives you the derivative for
non-multiple eigenvalues,
> if you plug the characteristic polynomial of the matrix.
The main trouble is that you do not get an eigenvector but the
whole straight line in its direction: every non-zero multiple
of
an eigenvector is again an eigenvector.
Normalization does narrow down the choices to two (mutually
opposite, with no clear decision procedure), which leads to
the
following implicit equation between A and v (where v is the
transpose of v):
A * v - v * v * A * v = 0 , v non-zero
with the side effect that v * A * v is the eigenvalue.
Implicit differentiation of v with respect to A is a mess.
If you are interested only in pure existence, or are not
afraid
of inverting and integrating to obtain derivatives, here is a
different philosophy:
The projector onto the eigenspace corresponding to an
eigenvalue
lambda is uniquely determined and has a formula due to Riesz:
E = (1/(2*pi*i)) * int[along C] (z*I - A)^(-1) dz
where the smooth curve C runs once around lambda and misses
all
other eigenvalues.
The rank of E equals the multiplicity of lambda.
Differentiation: To get the total differential of E applied to
an increment H in A, you differentiate the resolvent
(z*I-A)^(-1) with respect to H, and integrate. The result is
if H = dA then (abbreviating to avoid long formulas)
dE = (1/(2*pi*i)) * int[along C] W(z) dz
and W(z) = (z*I-A)^(-1) * H * (z*I-A)^(-1)
If the curve C is a circle, polar coordinates centered at
lambda
lead to the integration of a smooth periodic function over
full
period, so Trapezoidal Rule converges very fast.
The eigenvectors are recaptured as a basis of the range of E.
Another welcome property of Riesz formula is that it applies
to
a cluster of eigenvalues, yielding a basis of the invariant
subspace corresponding to that cluster, and also symmetry is
not
required. C can be a sum of simple curves, if convenient
(e.g.in
case of pairs of complex conjugate eigenvelues).
More detailed discussion is found in
H. Baumgaertel (or Baumgartel?)
Analytic perturbation theory for matrices and operators.
Birkhauser Verlag, Basel, 1985.
===
Subject: Percentile Estimation in a Monte Carlo Simulation
I am trying to estimate a fairly extreme percentile (99.7)
for a
distribution obtained from Monte Carlo simulation. I am
facing a
problem because of the instability of data at the extremes of
the
distribution, as in, if the percentile values are calculated
for
different simulation runs, they seem to jump around a lot. I
am
currently using the MATLAB prctile function to estimate the
percentile
and am limited to estimate it with 1000 iterations.
I would really appreciate if anybody can suggest me a better
algorithm
than what is being used in the prctile function, which can
provide a
closer estimate to the true percentile.
Uttam
===
Subject: SPSS.v12.0, SPSS.Web.Deployment.Framework.v2.4, -
new !
SPSS.v12.0, SPSS.Web.Deployment.Framework.v2.4, - new !
for more info please send email
===
Subject: poisson equation
I am trying to solve poisson equation:
diff(diff(u(y,z),y),y) + diff(diff(u(y,z),z),z) = 1
using separation of variables method.
The boundary conditions for the problem are:
u(1,z)=0
u(-1,z)=0
u(y,1)=0
u(y,-1)=0
can you please guide me to a good reference which will help
me in
solving this problem using separation of variables.
Do you know how to solve this equation using mathematica
===
Subject: Fourier transforms (coefficient calculation)...
For a few days now, I have tried a number of methods that are
supposed
to provide the a(k) and b(k) coefficients for a Fourier
series. These
methods are listed at the end of this post (in Java).
However, not
one of these methods seems to provide the correct coefficients
for the
following function;
f(x) = 2 - 2 * cos(x)
...I never see a Ô-2 or
Ô2 in the resulting generated
arrays. I get
a(0) = 4 which is correct, but there is no sign of a(1)
anywhere (when
the expression is changed to f(x) = 2 - 2 * sin(x), things
seem to
work).
Does someone have an algorithm that works for the above? I
have tried
all the Google searches (the source of the methods below),
but to no
avail. Please post some example code either in Java or C++ as
all the
Ôtry searching for ... suggestions seen in
other postings
have been
fruitless. In the meanwhile, perhaps the stuff below will be
of use to
someone.
Matthew
public static void computeFFT(double ar[], double ai[])
{
if (ar.length != ai.length)
{
System.err.println(array dimensions must
match);
}
else if (!isPowerOfTwo(ar.length))
{
System.err.println(array dimensions must be
multiple of 2);
}
else
{
computeFFT(1, ar.length, ar, ai);
if (ai[0] > EPSI)
{
System.err.println(imaginary part of
constant not zero);
}
}
}
public static void computeFFT(int sign, int n, double ar[],
double
ai[])
{
double scale = 2.0 / (double) n;
int i, j;
for (i = j = 0; i < n; ++i)
{
if (j >= i)
{
double tempr = ar[j] * scale;
double tempi = ai[j] * scale;
ar[j] = ar[i] * scale;
ai[j] = ai[i] * scale;
ar[i] = tempr;
ai[i] = tempi;
}
int m = n / 2;
while ((m >= 1) && (j >= m))
{
j -= m;
m /= 2;
}
j += m;
}
int mmax, istep;
for (mmax = 1, istep = 2 * mmax;
mmax < n;
mmax = istep, istep = 2 * mmax)
{
double delta = sign * Math.PI / (double) mmax;
for (int m = 0; m < mmax; ++m)
{
double w = m * delta;
double wr = Math.cos(w);
double wi = Math.sin(w);
for (i = m; i < n; i += istep)
{
j = i + mmax;
double tr = wr * ar[j] - wi *
ai[j];
double ti = wr * ai[j] + wi *
ar[j];
ar[j] = ar[i] - tr;
ai[j] = ai[i] - ti;
ar[i] += tr;
ai[i] += ti;
}
}
mmax = istep;
}
}
public static double[][] fft_1d(double[][] array)
{
double u_r, u_i, w_r, w_i, t_r, t_i;
int ln, nv2, k, l, le, le1, j, ip, i, n;
n = array.length;
ln = (int) (Math.log((double) n) / Math.log(2) + 0.5);
nv2 = n / 2;
j = 1;
for (i = 1; i < n; i++)
{
if (i < j)
{
t_r = array[i - 1][0];
t_i = array[i - 1][1];
array[i - 1][0] = array[j - 1][0];
array[i - 1][1] = array[j - 1][1];
array[j - 1][0] = t_r;
array[j - 1][1] = t_i;
}
k = nv2;
while (k < j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}
for (l = 1; l <= ln; l++) /* loops thru stages */
{
le = (int) (Math.exp((double) l * Math.log(2)) +
0.5);
le1 = le / 2;
u_r = 1.0;
u_i = 0.0;
w_r = Math.cos(Math.PI / (double) le1);
w_i = -Math.sin(Math.PI / (double) le1);
for (j = 1;
j <= le1;
j++) /* loops thru 1/2 twiddle values per
stage */
{
for (i = j;
i <= n;
i += le) /* loops thru points per
1/2 twiddle */
{
ip = i + le1;
t_r = array[ip - 1][0] * u_r - u_i *
array[ip - 1][1];
t_i = array[ip - 1][1] * u_r + u_i *
array[ip - 1][0];
array[ip - 1][0] = array[i - 1][0] -
t_r;
array[ip - 1][1] = array[i - 1][1] -
t_i;
array[i - 1][0] = array[i - 1][0] +
t_r;
array[i - 1][1] = array[i - 1][1] +
t_i;
}
t_r = u_r * w_r - w_i * u_i;
u_i = w_r * u_i + w_i * u_r;
u_r = t_r;
}
}
return array;
}
public static void fft(double[][] array)
{
double u_r, u_i, w_r, w_i, t_r, t_i;
int ln, nv2, k, l, le, le1, j, ip, i, n;
n = array.length;
ln = (int) (Math.log((double) n) / Math.log(2) + 0.5);
nv2 = n / 2;
j = 1;
for (i = 1; i < n; i++)
{
if (i < j)
{
t_r = array[i - 1][0];
t_i = array[i - 1][1];
array[i - 1][0] = array[j - 1][0];
array[i - 1][1] = array[j - 1][1];
array[j - 1][0] = t_r;
array[j - 1][1] = t_i;
}
k = nv2;
while (k < j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}
/* loops thru stages */
for (l = 1; l <= ln; l++)
{
le = (int) (Math.exp((double) l * Math.log(2)) +
0.5);
le1 = le / 2;
u_r = 1.0;
u_i = 0.0;
w_r = Math.cos(Math.PI / (double) le1);
w_i = -Math.sin(Math.PI / (double) le1);
/* loops thru 1/2 twiddle values per stage */
for (j = 1; j <= le1; j++)
{
/* loops thru points per 1/2 twiddle */
for (i = j; i <= n; i += le)
{
ip = i + le1;
t_r = array[ip - 1][0] * u_r - u_i *
array[ip - 1][1];
t_i = array[ip - 1][1] * u_r + u_i *
array[ip - 1][0];
array[ip - 1][0] = array[i - 1][0] -
t_r;
array[ip - 1][1] = array[i - 1][1] -
t_i;
array[i - 1][0] = array[i - 1][0] +
t_r;
array[i - 1][1] = array[i - 1][1] +
t_i;
}
t_r = u_r * w_r - w_i * u_i;
u_i = w_r * u_i + w_i * u_r;
u_r = t_r;
}
}
}
public static void newCompute(int sign, int n, double ar[],
double
ai[])
{
double scale = 2.0 / (double) n;
int i, j;
for (i = j = 0; i < n; ++i)
{
if (j >= i)
{
double tempr = ar[j] * scale;
double tempi = ai[j] * scale;
ar[j] = ar[i] * scale;
ai[j] = ai[i] * scale;
ar[i] = tempr;
ai[i] = tempi;
}
int m = n / 2;
while ((m >= 1) && (j >= m))
{
j -= m;
m /= 2;
}
j += m;
}
int mmax, istep;
for (mmax = 1, istep = 2 * mmax;
mmax < n;
mmax = istep, istep = 2 * mmax)
{
double delta = sign * Math.PI / (double) mmax;
for (int m = 0; m < mmax; ++m)
{
double w = m * delta;
double wr = Math.cos(w);
double wi = Math.sin(w);
for (i = m; i < n; i += istep)
{
j = i + mmax;
double tr = wr * ar[j] - wi *
ai[j];
double ti = wr * ai[j] + wi *
ar[j];
ar[j] = ar[i] - tr;
ai[j] = ai[i] - ti;
ar[i] += tr;
ai[i] += ti;
}
}
mmax = istep;
}
}
public static void otherCompute(int sign, int n, double ar[],
double
ai[])
{
double scale = 2.0 / (double) n;
int i, j;
for (i = j = 0; i < n; ++i)
{
if (j >= i)
{
double tempr = ar[j] * scale;
double tempi = ai[j] * scale;
ar[j] = ar[i] * scale;
ai[j] = ai[i] * scale;
ar[i] = tempr;
ai[i] = tempi;
}
int m = n / 2;
while ((m >= 1) && (j >= m))
{
j -= m;
m /= 2;
}
j += m;
}
int mmax, istep;
for (mmax = 1, istep = 2 * mmax;
mmax < n;
mmax = istep, istep = 2 * mmax)
{
double delta = sign * Math.PI / (double) mmax;
for (int m = 0; m < mmax; ++m)
{
double w = m * delta;
double wr = Math.cos(w);
double wi = Math.sin(w);
for (i = m; i < n; i += istep)
{
j = i + mmax;
double tr = wr * ar[j] - wi *
ai[j];
double ti = wr * ai[j] + wi *
ar[j];
ar[j] = ar[i] - tr;
ai[j] = ai[i] - ti;
ar[i] += tr;
ai[i] += ti;
}
}
mmax = istep;
}
}
public static Ret fourn(
double x_re[],
double x_im[],
int nn[],
int ndim,
int isign)
{
Ret ret = new Ret();
ßoat[] data = new ßoat[2 * nn[0]];
for (int i = 0; i < 2 * nn[0]; i = i + 2)
{
if (i < 2 * x_re.length && i < 2 * x_im.length)
{
data[i] = (ßoat) x_re[i / 2];
data[i + 1] = (ßoat) x_im[i / 2];
}
else
{
data[i] = 0;
data[i + 1] = 0;
}
}
//Replaces data by its ndim-dimensional discreate Fourier
transform,
//if isign is input as 1. nn[0..ndim-1] is an integer
array
containing
//the lengths of each dimension (number of complex values),
which
//MUST all be power of 2. data is a real array of
length twice
the
//product of these lengths, in which the data are stored as
in a
//multidimensional complex array: real and imaginary parts
of each
//element are in consecutive locations, and the right most
index of
//the array inreases most rapidly as one proceeds along
data. For a
//two-dimensionalbe array, this is equivalent to storing the
array
by
//rows. If isign is input as -1, data is replaced by its
inverse
//transform times the product of the lengths of all
dimensions.
int idim;
int i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
int ibit, k1, k2, n, nprev, nrem, ntot;
ßoat tempr, tempi;
//Double precision for the trigonometric recurrences
double wtemp, wr, wpr, wpi, wi, theta;
//Compute total number of complex values
for (ntot = 1, idim = 0; idim < ndim; idim++)
ntot *= nn[idim];
nprev = 1;
for (idim = ndim - 1; idim >= 0; idim--)
{
n = nn[idim];
nrem = ntot / (n * nprev);
ip1 = nprev << 1;
ip2 = ip1 * n;
ip3 = ip2 * nrem;
i2rev = 1;
//This is the bit-reversal section of the routine.
for (i2 = 1; i2 <= ip2; i2 += ip1)
{
if (i2 < i2rev)
{
for (i1 = i2; i1 <= i2 + ip1 - 2; i1
+= 2)
{
for (i3 = i1; i3 <= ip3; i3
+= ip2)
{
i3rev = i2rev + i3 -
i2;
ßoat temp = data[i3
- 1];
data[i3 - 1] =
data[i3rev - 1];
data[i3rev - 1] =
temp;
temp = data[i3];
data[i3] =
data[i3rev];
data[i3rev] = temp;
}
}
}
ibit = ip2 > 1;
while (ibit >= ip1 && i2rev > ibit)
{
i2rev -= ibit;
ibit >= 1;
}
i2rev += ibit;
}
//Here begins the Danielson-Lanczos section of the
routine.
ifp1 = ip1;
while (ifp1 < ip2)
{
ifp2 = ifp1 << 1;
// Initialized for the trig. recurrence
theta = isign * 6.28318530717959 / (ifp2 /
ip1);
wtemp = Math.sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = Math.sin(theta);
wr = 1.0;
wi = 0.0;
for (i3 = 1; i3 <= ifp1; i3 += ip1)
{
for (i1 = i3; i1 <= i3 + ip1 - 2; i1
+= 2)
{
for (i2 = i1; i2 <= ip3; i2
+= ifp2)
{
//Danielson-Lanczos
formula:
k1 = i2;
k2 = k1 + ifp1;
tempr =
(ßoat) wr *
data[k2
-
1]
-
(ßoat) wi * data[k2];
tempi =
(ßoat) wr *
data[k2]
+
(ßoat) wi * data[k2
-
1];
data[k2 - 1] =
data[k1 - 1] - tempr;
data[k2] = data[k1]
- tempi;
data[k1 - 1] +=
tempr;
data[k1] += tempi;
}
}
//Trigonometric recurrence
wr = (wtemp = wr) * wpr - wi * wpi +
wr;
wi = wi * wpr + wtemp * wpi + wi;
}
ifp1 = ifp2;
}
nprev *= n;
}
ret.x_re = new double[data.length / 2];
ret.x_im = new double[data.length / 2];
for (int i = 0; i < data.length; i = i + 2)
{
ret.x_re[i / 2] = data[i];
ret.x_im[i / 2] = data[i + 1];
}
return ret;
}
public class Ret
{
public double[] x_re;
public double[] x_im;
}
===
Subject: Re: Fourier transforms (coefficient calculation)...
>Does someone have an algorithm that works for the above? I
have tried
>all the Google searches (the source of the methods below),
but to no
>avail. Please post some example code either in Java or C++
as all the
>try searching for ... suggestions seen in
other postings
have been
>fruitless. In the meanwhile, perhaps the stuff below will be
of use to
>someone.
This is off-topic in at least two ways here, but I feel nice
today, so
I tried (after making it compilable) your method number I
with the
following main function.
public static void main(String argv[])
{ final int size=8;
double ar[]=new double[size];
double ai[]=new double[size];
for(int i=0;i
>
> For a few days now, I have tried a number of methods that
are supposed
> to provide the a(k) and b(k) coefficients for a Fourier
series. These
> methods are listed at the end of this post (in Java).
However, not
> one of these methods seems to provide the correct
coefficients for the
> following function;
>
> f(x) = 2 - 2 * cos(x)
>
> ...I never see a Ô-2 or
Ô2 in the resulting generated
arrays. I get
> a(0) = 4 which is correct, but there is no sign of a(1)
anywhere (when
> the expression is changed to f(x) = 2 - 2 * sin(x), things
seem to
> work).
>
> Does someone have an algorithm that works for the above? I
have tried
> all the Google searches (the source of the methods below),
but to no
> avail. Please post some example code either in Java or C++
as all the
> Ôtry searching for ... suggestions seen in
other postings
have been
> fruitless. In the meanwhile, perhaps the stuff below will
be of use to
> someone.
>
This is Off Topic
in comp.lang.c comp.lang.c++ and probably comp.lang.java.
Try The Object-Oriented Numerics Page
http://www.oonumerics.org/oon/
GSL -- The GNU Scientific Library
http://sources.redhat.com/gsl/
Available C++ Libraries FAQ
http://www.faqs.org/faqs/C++-faq/libraries/part1/
===
Subject: Re: Computational Evolution without velocities
>* Gordon D. Pusch:
> Newtons laws (i.e. non-relativistic) one can neglect the
fact that the
> update formula
>
> x_new = x_old - Force DeltaTime / 2 Mass
>
> holds for any simulation and we do not need to take care of
velocities.
>> This equation is quite obviously wrong, [...]
>He may be talking about Stoermer-Verlet algorithms,
>which have update formulae that are indeed relatively
>simple and dont contain the velocity. Also, they are
>used for n-body problems and hard sphere ßuids and have
>other interesting properties (look up symplectic algorithms).
>x_new = 2 x_old - x_older + DeltaTime^2 Force / Mass +
O(DeltaTime^4)
>would be a correct formula. That is not far off from his
>version, and given that he heard it as folklore, it may
>well be it.
>(Of course, the velocity is not neglected here. It is
implicitly
>contained in the ODE of 2nd order that is discretized
directly.
>It could also be rewritten into two ODEs of 1st order and
discretized
>afterwards, which would lead to better known schemes like
Runge-Kutta.)
I suspect that you are right. And, after reading the thread
here I suspect that some folks who answered were a bit too
quick on the trigger.
First, the absence of a velocity (first derivative) term in
what
is really a second order differential equation is no problem.
There are several well-known methods for doing such
integrations
when the velocity is not needed.
There are several different Stoermer-Verlet algorithms. One
is given above by Mazzucco. The original and slightly
incorrect
formula given by the original poster is half of what is
usually
called velocity Verlet. There are a pair of coupled equations
x_{n+1} = x_n + hv_n + h^2 F(x_n) / 2 m
v_{n+1} = v_n + h [ F(x_{n+1} + F(x_n) ] / 2 m
where x is the coordinate, v the velocity, and F the force
evaluated at the given argument. Starting with x_0 and v_0
the first equation can be evaluated which allows for the
second to be evaluated since the needed coordinate is now
known. So this is NOT a predictor-corrector formula.
The set of equations is symplectic and while not of high
order, if the system is conservative and describle by a
Hamiltonian in which the kinetic energy depends only on
the velocities and the potential eneryg only on the
coordinates,
then the total energy will not deviate much from its
initial value. That is the deviation in the total energy
is bounded.
Since there is only once function evaluation per cycle
through the equations (after the first trip) this is a
fairly cheap integration method. It is, as the original
poster said, quite popular among those who do modelling
of galaxies and molecules.
There is a large literature on these equations and their
several cousins. A search on Verlet and symplectic
on Google turned up 528 responses, many of which are
quite good.
---- Paul J. Gans
===
Subject: Re: Eignevector differentiation
>
> Would that give me the derivatives of the eigenvectors as
well?
the procedure has already been explained by Arnold Neumaier;
but be
aware that you have to carry out the procedure for each
eigenvector. Since you might be only interested in the
dependency of
the solution vector x on A, you could apply the same method
to your
equation.
In order to find the dependency of all eigenvalues and
eigenvectors on
A, you could apply V.I. Arnolds Normal Form
for matrices,
which has
been outlined in his book on Geometric methods.
> However as AA it is close to singular I should use the
Singular Value
> Decompostion algorithm to invert AA. Instead of SVD, for
> efficiency, I am
I seriously doubt that that would improve performance! As far
as I
know SVD is much cheaper than finding the eigenvalues and
eigenvectors!
Another hint: In almost any textbook you will find the
recommendation
not to use AA, but perform the SVD directly to A.
Maybe you should also apply bifurcation theory to your
problem, since
at to the singularity you just dont have any continuous
dependence.
Also the SVD doesnt behave smoothly, when the singular
values cross
the threshold.
Alois
--
Alois Steindl, Tel.: +43 (1) 58801 / 32558
Inst. for Mechanics II, Fax.: +43 (1) 58801 / 32598
Vienna University of Technology,
A-1040 Wiedner Hauptstr. 8-10
===
Subject: Re: Percentile Estimation in a Monte Carlo Simulation
> I am trying to estimate a fairly extreme percentile (99.7)
for a
> distribution obtained from Monte Carlo simulation. I am
facing a
> problem because of the instability of data at the extremes
of the
> distribution, as in, if the percentile values are
calculated for
> different simulation runs, they seem to jump around a lot.
I am
> currently using the MATLAB prctile function to estimate the
percentile
> and am limited to estimate it with 1000 iterations.
> I would really appreciate if anybody can suggest me a better
algorithm
> than what is being used in the prctile function, which can
provide a
> closer estimate to the true percentile.
> Uttam
(i) You can work out roughly how many samples you would need
to get
reasonable accuracy. Treat the quantile as a fixed value and
the
number of samples exceeding this as a random variable ...
essentially
Binominal. The estimate of the probability is a scaled
version of
this, and hence you can get a formula for the standard
deviation of
the estimated probability as a function of sample size. To
get a
reasonable estimate of the probability when it is near 0.003,
you
might decide that you should estimate it to within plus or
minus
0.0001 .... following through with this will show you need a
very
large number of simulations.
(ii) To improve matters you will need to move away from
straightfoward
direct simulations and make use of some knowledge of the
structure of
whatever you are simulating. A number of general approaches
are
available ...
(a) Variance reduction techniques (qv). ( These are more often
suggested for problems of estimating mean values, but you
might fins
something useful under this heading). Examples are the use of
antithetic pseudo-random numbers and using control variables.
(b) Importance sampling (qv). You fix up the simulations to
generate
more of the extreme events and then downscale the estimated
frequency.
This can result in much improved estimates of the tails of a
distribution for the same number of simulations.
(c) Remove some of the random variation from the simulations.
For
example, if your final variable is generated from a
distribution
conditional on another random variable, you might remove this
step
(the generation of the final variable) and replace it with a
theoretical calculation based on the conditional distribution.
David Jones
===
Subject: Discrete Orthogonal Polynomials - II
A few days i ago i asked a question regarding discrete
Chebyshev
polynomials and the answers were usefull in the sense that it
triggered
me to rethink the problem.
The problem is as follows.
Ive a set of measurements sampled on a uniform (1D or 2D
-grid).
Through these values I want to estimate a function that
describes the
values I have.
Using the estimated function I can obtain a set of equations
that
contain certain parameters Im interested in.
The problem lies in the parameterization of the function. I
Know it is a
fairly smooth function so I can parameterize the function in
terms of a
Taylor series. Instead of a Taylor-series I can also
parameterize it as
a orthogonal polynomial,for example a Chebyshev polyn. But in
order to
do this correctly and to be able to use al nice properties of
the
Chebyshev polynomials I need sampled points on a
Chebyshev-grid. This
is non-uniform and for each order of the Chebychev polynomial
different.
But in literature (Mostly Image Processing) I see examples in
which they
apply Chebychev polynomials on a uniform grid, an image for
example. And
they use the scalar-product rules but now using the uniform
sampled points.
Now are my questions
1 Is this allowed? (I think the answer is no).
2 Can I overcome the grid-problem for the Chebyshev
polynomials?
3 There are a lot of other discrete orthogonal polynomials, is
there one in which I can directly use the points sampled on
the
uniform grid?
Maurice
===
Subject: quadrature
Could somebody help me to answer the following question: Let
$fcolon[0,infty)to[0,infty)$. What do I have to assume in
order to
obtain the convergence of the square quadrature for $f$? Id
be
grateful for any help. Peter
===
Subject: Re: poisson equation
> I am trying to solve poisson equation:
> diff(diff(u(y,z),y),y) + diff(diff(u(y,z),z),z) = 1
> using separation of variables method.
> The boundary conditions for the problem are:
> u(1,z)=0
> u(-1,z)=0
> u(y,1)=0
> u(y,-1)=0
> can you please guide me to a good reference which will help
me in
> solving this problem using separation of variables.
> Do you know how to solve this equation using mathematica
You cannot solve nonhomogeneous equations with separation of
variables.
Try
instead an eigenvector expansion, e.g.,
u(y,z)=Sum_{mn}A_{mn}cos[(2m-1)pi y/2]cos[(2n-1)pi z/2].
Substitute this equation into your PDE, and solve for
coefficients A_{mn}.
See any standard book on Fourier series solution of PDEs.
===
Subject: Re: quadrature
>
> Could somebody help me to answer the following question: Let
> $fcolon[0,infty)to[0,infty)$. What do I have to assume in
order
to
> obtain the convergence of the square quadrature for $f$?
Id be
> grateful for any help. Peter
I presume you are asking what conditions must f(x) satisfy in
order that
the
(improper) integral --LaTeX format, sorry--
int_0^infty {dxleft| {ßeft( x right)} right|^2 }
be finite.
Obviously, assuming f(x) has no singularities other than 0
and infty,
you want |f|^2 to be integrable at x = 0, and also it should
fall off
rapidly enough for large x as to be integrable there also.
Since this sounds like a homework problem, I will content
myself with
hints:
What is the condition on c < 0 for x^c to give a finite
integral
near x = 0?
What is the condition on c < 0 for x^c to yield a finite
integral
as x -> infty ?
How can I translate these bounds to conditions on |f(x)| ?
--
Julian V. Noble
Professor Emeritus of Physics
^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
Science knows only one commandment: contribute to science.
-- Bertolt Brecht, Galileo.
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
>
>>2. Abramowitz & Stegun, Handbook of Mathematical Functions
>> (Dover PB)
Or download it for free.
> http://jove.prohosting.com/~skripty/
Any evidence that this book is out of copyright, and that
its now okay
> to post it on the web?
>
> The Dover version might be copyright, but the hardback is a
US
> government publication, and hence (IIRC) not
copyright-protected.
>
>
> Phil Hobbs
Quite right. However, keep in mind that by the time you print
out this
immense
book, you will have spent a lot more than the Dover PB will
run you. I
dont
know if the US Printing Office still puts out a hardcover
version.
--
Julian V. Noble
Professor Emeritus of Physics
^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
Science knows only one commandment: contribute to science.
-- Bertolt Brecht, Galileo.
===
Subject: Re: Discrete Orthogonal Polynomials - II
>
>
> A few days i ago i asked a question regarding discrete
Chebyshev
> polynomials and the answers were usefull in the sense that
it triggered
> me to rethink the problem.
> The problem is as follows.
> Ive a set of measurements sampled on a uniform (1D or 2D
-grid).
> Through these values I want to estimate a function that
describes the
> values I have.
> Using the estimated function I can obtain a set of
equations that
> contain certain parameters Im interested in.
> The problem lies in the parameterization of the function. I
Know it is a
> fairly smooth function so I can parameterize the function
in terms of a
> Taylor series. Instead of a Taylor-series I can also
parameterize it as
> a orthogonal polynomial,for example a Chebyshev polyn. But
in order to
> do this correctly and to be able to use al nice properties
of the
> Chebyshev polynomials I need sampled points on a
Chebyshev-grid. This
> is non-uniform and for each order of the Chebychev
polynomial different.
> But in literature (Mostly Image Processing) I see examples
in which they
> apply Chebychev polynomials on a uniform grid, an image for
example. And
> they use the scalar-product rules but now using the uniform
sampled
points.
> Now are my questions
> 1 Is this allowed? (I think the answer is no).
> 2 Can I overcome the grid-problem for the Chebyshev
polynomials?
> 3 There are a lot of other discrete orthogonal polynomials,
is
> there one in which I can directly use the points sampled on
the
> uniform grid?
>
> Maurice
You evidently havent looked at my class notes. Let me
STRONGLY suggest you
do--they should answer all your questions, including what to
do if the
grid or the weights (in the inner product) are non-uniform.
(In that case
you
can still define and use orthogonal polynomials. They are
called Gram poly-
nomials and are by far the best way to do high-order linear
least-squares
fits.
--
Julian V. Noble
Professor Emeritus of Physics
^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
Science knows only one commandment: contribute to science.
-- Bertolt Brecht, Galileo.
===
Subject: Re: Information about Wavelet Analysis
> I am looking for information (book name internet source)
about Wavelet
analysis.
Try http://www.ondelette.com/
HTH
===
Subject: Re: poisson equation
>
> I am trying to solve poisson equation:
>
> diff(diff(u(y,z),y),y) + diff(diff(u(y,z),z),z) = 1
>
> using separation of variables method.
>
> The boundary conditions for the problem are:
>
> u(1,z)=0
> u(-1,z)=0
> u(y,1)=0
> u(y,-1)=0
>
> can you please guide me to a good reference which will help
me in
> solving this problem using separation of variables.
>
> Do you know how to solve this equation using mathematica
Try Mathews & Walker, Mathematical Methods of Physics for an
accessible introduction.
--
Julian V. Noble
Professor Emeritus of Physics
^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
Science knows only one commandment: contribute to science.
-- Bertolt Brecht, Galileo.
===
Subject: Re: Books On Generalized Reduced Gradient Algorithm?
-Vik
===
Subject: [newbie] detecting singularity in matrices
The criterion SAS-IML uses for dectecting singular matrices
(in
functions meant for solving linear equations) is that, while
reducing
the orginal matrix to a triangular/diagonal form, all
elements less
than 100 x machine epsilon x A
(where A is largest element of the input matrix in absolute
value)
are zeroed.
What is the reason for this criterion? And why is the
constant 100
chosen?
gc
===
Subject: Re: Discrete Chebyshev (or other orthogonal)
polynomials.
>
> For my research I want to find out more about discete
Chebyshev
> polynomials of the first kind and their relation to
Chebyshev
> polynomials of the second kind.
> In literature a lot of properties of the Chebyshev
polynomials and the
> relations to other polynomials are known. But, afaik, all
these nice
> propertie are derived under the assumption that we want to
approach a
> continue function. I cannot find a list of properties in
case we want to
> approach a discrete set of equidistant spaced points with
these
polynomials.
> The description in Numerical Recipes for example also
assumes the the
> function is known at all locations.
> interesting words but nothing came out of it.
>
Regarding discrete orthogonal polynomials , see the monograph
[1] A.F. NIKIFOROV , S.K. SUSLOV , V.B. UVAROV ,
,, Classical Polynomials of a Discrete Variable ,
Springer Series in Computational Physics , Springer-Verlag
,1991.
In this book there are a lot of interesting facts about the
case when
the ,,mesh is uniform, that is when we suppose that points are
equidistant.
There are also many references regarding applications, for
instance in
Physics.
For general properties of (continuous) Chebychev polynomials,
see
[2] Theodore J. RIVLIN , ,,The Chebyshev Polynomials , Wiley -
Interscience , 1974 .
===
Subject: Czy ktos wie jak sie tworzy diagramy krat???
Pozdrawiam
Justyna
===
Subject: positive definite ?=> diagonally dominant?
Can positive definite results in diagonally dominant? Are
there any
===
Subject: Re: positive definite ?=> diagonally dominant?
> Can positive definite results in diagonally dominant? Are
there any
Linear finite elements on a Poisson problem gives diagonally
dominant pd
matrices; higher order elements are still pd but not
diagonally
dominant.
V.
===
Subject: positive definite matrix
Is the diagonal element of a positive definite matrix the
bigest one
comparing with others in the same roe? that is, is a_{ii} >=
|a_{ij}|?
===
Subject: Re: positive definite matrix
>
> Is the diagonal element of a positive definite matrix the
bigest one
> comparing with others in the same roe? that is, is a_{ii}
>= |a_{ij}|?
Not in general, A=[1 2;2 4].
Yes for Hermitian matrices with constant diagonal, since the
determinant of every principal 2x2 minor is nonnegative.
Arnold Neumaier
===
Subject: Re: positive definite matrix
> Is the diagonal element of a positive definite matrix the
bigest one
> comparing with others in the same roe? that is, is a_{ii}
>= |a_{ij}|?
Take the matrix A=
(2 3)
(0 2)
This is positive definite, since
(x,y)*A*(x,y)^T=2*(x+3y/4)^2+7y^2/8>0