mm-290
==
Subject: Re: ill conditioned non-lin system again
===
>There
are between 1 and 6, say k=1..m differential
equations:>dy1(x)/dx=-yk(x)*exp(a1-a2/x) with initial value
x=x0 and y1(x0)=a0>dy2(x)/dx=-yk(x)*exp(a4-a5/x) with initial
value x=x0 and y2(x0)=a3. etc (x is actually (a function of)
temperature in time during >heating of a sample)Do you really
mean dy/dx=y*exp(a-b/x)or dy/dx=y*(a-b/x)?The former equation
has a solution which looks nothing at all like what
youoriginally pos. The second equation looks much more like
theequations that I've seen in chemistry. Are you solving the
differential equations numerically, or usingan exact analytic
solution? Parameter estimation for differential equations is a
problem whichhas been studied extensively. >Each diffeqn
results in a skew bell-shaped function resulting in a sum
>that has a bumpy bell shape.>In practice the a1 is of the
same size as a2/x and the final x is about >2*x0. This means
that increasing both a1 and a2 makes no dramatic change >and
can easily be compensa by small changes in the parameters of
an >other of the diffeqns.This shows in an intuitive way that
you've got a badly conditionedproblem. >The suggestion of
adding a fraction of the sum of squares of the >parameters as
a kind of penalty to the ssq clearly improves the system. >In
addition it is even physically relevant because in particular
the a1, >a4 .. should not be too large to have a physical
meaning.This is known as ridge regression It's also rela to
Tikhonovregularization. Unfortunately, the estimates that you
get from thisprocedure are biased, and there's no way to get a
bound on the sizeof the bias. Thus you can't get usable
confidence intervals for the estima parameters. i.e. you get
solutions like a1=10, plus or minus50, plus or minus another
term that might be 1000 or more, but we don'tknow how large it
is.-- Brian Borchers borchers@nmt.eduDepartment of Mathematics
http://www.nmt.edu/~borchers/New Mexico Tech Phone:
505-835-5813Socorro, NM 87801 FAX: 505-835-5366 -----= Pos via
Newsfeeds.Com, Uncensored Usenet News
=-----http://www.newsfeeds.com - The #1 Newsgroup Service in
the World!-----== Over 100,000 Newsgroups - 19 Different
Servers! =-----
====

Subject: Re: ill conditioned non-lin system
againThanks for paying attention to my problem,I did indeed
mean dy/dx=y*exp(a-b/x). In this case the solution contains
the exponential integral and in many cases it should actually
bedy/dx=y*exp(a-b/z(x)) where z(x) is of the form
p(1-exp(-qx)). So there is no way around but solving
numerically. I have tried many procedures with Burlish-Stoer
and trapezoidal-Newton as best. In case z(x) has an
exponential shape only trap-Newton works; all others tend to
oscillate.Yes you would intuitively think that the system is
not too well behaved but that lambda_max/lambda_min would be
something like >1E10 (in the one diffeqn case e.g.
lambda={250,1E-7,1E-9})A plot of ssq(a1,a2) shows a narrow
valley with an only slightly curved bottom in the
length-direction.>>There are between 1 and 6, say k=1..m
differential equations:>>dy1(x)/dx=-yk(x)*exp(a1-a2/x) with
initial value x=x0 and
y1(x0)=a0>>dy2(x)/dx=-yk(x)*exp(a4-a5/x) with initial value
x=x0 and y2(x0)=a3>. etc (x is actually (a function of)
temperature in time during >>heating of a sample)> Do you
really mean > dy/dx=y*exp(a-b/x)> or> dy/dx=y*(a-b/x)> ?>
The former equation has a solution which looks nothing at all
like what you> originally pos. The second equation looks much
more like the> equations that I've seen in chemistry. > Are
you solving the differential equations numerically, or using>
an exact analytic solution? > Parameter estimation for
differential equations is a problem which> has been studied
extensively. >>Each diffeqn results in a skew bell-shaped
function resulting in a sum >>that has a bumpy bell
shape.>>In practice the a1 is of the same size as a2/x and
the final x is about >>2*x0. This means that increasing both
a1 and a2 makes no dramatic change >>and can easily be
compensa by small changes in the parameters of an >>other of
the diffeqns.> This shows in an intuitive way that you've got
a badly conditioned> problem. >>The suggestion of adding a
fraction of the sum of squares of the >>parameters as a kind
of penalty to the ssq clearly improves the system. >>In
addition it is even physically relevant because in particular
the a1, >>a4 .. should not be too large to have a physical
meaning.> This is known as ridge regression It's also rela to
Tikhonov> regularization. Unfortunately, the estimates that
you get from this> procedure are biased, and there's no way to
get a bound on the size> of the bias. Thus you can't get usable
confidence intervals for the > estima parameters. i.e. you get
solutions like a1=10, plus or minus> 50, plus or minus another
term that might be 1000 or more, but we don't> know how large
it is.
====

Subject: Re: ill conditioned non-lin system again
>Thanks for paying attention to my problem, >I did indeed
mean dy/dx=y*exp(a-b/x). In this case the solution contains
>the exponential integral and in many cases it should actually
be >dy/dx=y*exp(a-b/z(x)) where z(x) is of the form
p(1-exp(-qx)). So there >is no way around but solving
numerically. I have tried many procedures >with Burlish-Stoer
and trapezoidal-Newton as best. In case z(x) has an
>exponential shape only trap-Newton works; all others tend to
oscillate. >Yes you would intuitively think that the system
is not too well behaved > but that lambda_max/lambda_min would
be something like >1E10 (in the one >diffeqn case e.g.
lambda={250,1E-7,1E-9}) >A plot of ssq(a1,a2) shows a narrow
valley with an only slightly curved >bottom in the
length-direction. >>There are between 1 and 6, say
k=1..m differential equations: >dy1(x)/dx=-yk(x)*exp(a1-a2/x)
with initial value x=x0 and y1(x0)=a0
>dy2(x)/dx=-yk(x)*exp(a4-a5/x) with initial value x=x0 and
y2(x0)=a3 >>. etc (x is actually (a function of) temperature
in time during >heating of a sample) >> Do you really
mean > dy/dx=y*exp(a-b/x) > or > dy/dx=y*(a-b/x)
> ? > The former equation has a solution which looks
nothing at all like what you >> originally pos. The second
equation looks much more like the >> equations that I've seen
in chemistry. > Are you solving the differential equations
numerically, or using >> an exact analytic solution? >
Parameter estimation for differential equations is a problem
which >> has been studied extensively. >Each diffeqn
results in a skew bell-shaped function resulting in a sum
>that has a bumpy bell shape. >In practice the a1 is of the
same size as a2/x and the final x is about >2*x0. This means
that increasing both a1 and a2 makes no dramatic change >and
can easily be compensa by small changes in the parameters of
an >other of the diffeqns. >> This shows in an intuitive
way that you've got a badly conditioned >> problem. >The
suggestion of adding a fraction of the sum of squares of the
>parameters as a kind of penalty to the ssq clearly improves
the system. >In addition it is even physically relevant
because in particular the a1, >a4 .. should not be too large
to have a physical meaning. >> This is known as ridge
regression It's also rela to Tikhonov >> regularization.
Unfortunately, the estimates that you get from this >>
procedure are biased, and there's no way to get a bound on the
size >> of the bias. Thus you can't get usable confidence
intervals for the >> estima parameters. i.e. you get solutions
like a1=10, plus or minus >> 50, plus or minus another term
that might be 1000 or more, but we don't >> know how large it
is. >> this is parameter identification for ordinary
differential equations andillconditioned per se. did you
already try some well tes software like this one:
http://www.uni-bayreuth.de/departments/math/~kschittkowski/
demo_probs.htmyour observation that only trapezoidal rule is a
suitable integrator indicatesthat your system also exhibits
stiffness, which makes things even worse.if you are forced
doing it yourself, then you should use:a ver good integrator
(say radau5 of hairer and wanner).a very good nonlinear least
squares code using gauss newton with handling of rank
defective problems like elsunc (or enlsip ifyou need
restrictions). see
http://plato.la.asu.edu/topics/problems/nlolsq.htmlhthpeter
===

==
Subject: Re: ill conditioned non-lin system again> Thanks for
paying attention to my problem,> I did indeed mean
dy/dx=y*exp(a-b/x). In this case the solution contains> the
exponential integral and in many cases it should actually be>
dy/dx=y*exp(a-b/z(x)) where z(x) is of the form p(1-exp(-qx)).
So there> is no way around but solving numerically. I have
tried many procedures> with Burlish-Stoer and
trapezoidal-Newton as best. In case z(x) has an> exponential
shape only trap-Newton works; all others tend to oscillate.
You might want to try 2nd-order implicit Runge-Kutta. You have
an ^^^^^^^^essential singularity at x = 0, and rational
approximation schemes arenot likely to give you stable
behavior near that point.I have used this successfully for
equations with singularities.--
^^^^^^^^^^^^^^^^^^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: Indexing sparse
matrix by support1.mathforum.org (8.11.6/8.11.6/The Math
Forum, $Revision: 1.9 primary) id i0MJCxC28831;
===
Hi
Rick,Thank you very much for your feedback. It really worked
elegantly. I did infact change ij into ij-1 before I saw your
reply again. Itaught that was just a typo
mistake.Yacob
====

Subject: Complex numbers proofHi there,I cam
across to bold statements in a math's book but can't see
toprove that they are in fact true. (is the book wrong, no, i
reckonit's just me). Can anyone out there help?First
statement|cos(x+iy)|^2=1/2[cos2x+cosh2y]I can't seem to prove
this/Also, second statementexp ((ib)/2) 1--------------- =
------------exp (-ib) + 1 2 cos (b/2)I can't prove this
either.Any help would be much
apprecia,thanksJonathan
====

Subject: Re: Complex numbers proof>
Hi there,> I cam across to bold statements in a math's book
but can't see to> prove that they are in fact true. (is the
book wrong, no, i reckon> it's just me). Can anyone out there
help?> First statement> |cos(x+iy)|^2=1/2[cos2x+cosh2y] This
one is CORRECT. You prove it by using the following
definitions and identities (== means identical): cos(x+iy) ==
(1/2) * [ e^{ix} * e^{-y} + e^{-ix} * e^y ] | z |^2 == | z^2 |
cos^2 (a) + sin^2 (a) == 1 cosh^2 (b) = 1 + sinh^2 (b) | x+iy |
= sqrt( x^2 + y^2 ) (a + b)^2 = a^2 + 2a*b + b^2 As I don't
want to do your work for you, I leave the rest to you.> I
can't seem to prove this. You should be able to now. > Also,
second statement> exp ((ib)/2) 1> --------------- =
------------> exp (-ib) + 1 2 cos (b/2) ^ > I can't prove this
either. There's a good reason. You have a wrong sign, either in
the downstairs exponent (as I indica with a carat ^ ) or in the
upstairs exponent. > Any help would be much apprecia, Send
money. LOTS of money. Talk is cheap ;-) > thanks You're
welcome. > Jonathan--
^^^^^^^^^^^^^^^^^^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: function interpolation
by using normal distributionsis there an algorithm which
permits to abtain a function interpolationby using a set of
normal distributions?M
====

Subject: Re: first deriviate with
FFT, indexing problems>==
Subject: first deriviate with FFT,
indexing problems>Hey there,>I have n evenly spaced sample
values taken from a bandlimi function.>I assume that sampling
has been done properly (nyquist criteria).>I like to calculate
the first derivative of this function by using>the Derivative
Theorem of the Fourier transform:> df[n]/dt <=> -jnF[n]The
periodic discrete Fourier transform (what you get if you talk
FFT)is defined over a different group than the Fourier
transform you seem to be assuming as the continuous unbounded
Fourier transform is the variant which has the derivative
theorem that you quote. All of this assumes that you are doing
the mathematics at the level you indicate.Operationally you
will find that the DFT has a Finite Difference Theorem, rather
than a Derivative Theorem, that is easy to derive. It will tell
you exactly what form of of modification you need to doto get
the FT of the finite difference and then you will have totreat
the finite difference as an approximation to the
derivative.There are FOUR common Fourier transforms that you
are likely toencounter which are often named Fourier
transform, Fourier series, Fourier sequences and Discrete
Fourier transform depending upon whether the time index is
continuous or discrete and whether it is unbounded or
periodic. Most results will only correspond between the
variants rather than being exactly the same. You have been
tripped up by crossing the variant boundary without noticing
it.>- First I do a Fourier transform. As I only have real
values, the> resulting function is hermitian:> f(x) = f*(-x)
f* means conjugate complex>- Then the multiplication with -j
(-sqrt(-1)) is a rotation of 90deg around> the time axis. Easy
saying:> re_new = im_old> im_new = -re_old>- than I multiply
the result values with their index n> for example size 8:>
index: 0 1 2 3 4 5 6 7 > factor 0 1 2 3 k -3 -2 -1> and there
is the problem! what to do at index 4?> in order to keep the
hermitian property f(4)=f*(4)> I would have to force the
imaginary component to 0, but that would alter > the data and
therefore lead to a less precise results.> The reason why I
want to keep the hermitian property is, that I want> a real
function after the inverse Fourier transform.> Fortunately
this is only a problem when I have an even number of> samples,
so I can always zero pad to a odd number.> But that is not very
satisfying, and I would appreciate a cleaner solution.>Maybe
someone can help me out here!>Thanks a lot>Martin
====

Subject:
Re: first deriviate with FFT, indexing problems> Hey there,> I
have n evenly spaced sample values taken from a bandlimi
function.> I assume that sampling has been done properly
(nyquist criteria).> I like to calculate the first derivative
of this function by using> the Derivative Theorem of the
Fourier transform:> df[n]/dt <=> -jnF[n]> - First I do a
Fourier transform. As I only have real values, the> resulting
function is hermitian:> f(x) = f*(-x) f* means conjugate
complex> - Then the multiplication with -j (-sqrt(-1)) is a
rotation of 90deg around> the time axis. Easy saying:> re_new
= im_old> im_new = -re_old> - than I multiply the result
values with their index n> for example size 8:> index: 0 1 2 3
4 5 6 7> factor 0 1 2 3 k -3 -2 -1> and there is the problem!
what to do at index 4?> in order to keep the hermitian
property f(4)=f*(4)> I would have to force the imaginary
component to 0, but that would alter> the data and therefore
lead to a less precise results.> The reason why I want to keep
the hermitian property is, that I want> a real function after
the inverse Fourier transform.> Fortunately this is only a
problem when I have an even number of> samples, so I can
always zero pad to a odd number.> But that is not very
satisfying, and I would appreciate a cleaner solution.> Maybe
someone can help me out here!> Thanks a lot> MartinYou should
filter your input data so that there is no signal componentfor
the Nyquist frequency. That is in your case f(4)=0.The Nyquist
freq component is a signal of the type:+1, -1, +1, -1, +1, -1,
...How would you define the derivative of such sampled signal ?
Jean-Pierre
====

Subject: Re: Too many iterations in
tqli.cX-AUTHid: minutsil> I agree with you that there does not
seem to be any perfect solutions> for matrix computations in
C++ in the same way that LAPACK exists for> FORTRAN (and C). I
have developed my own class library whose primary> goal is to
interface to the core functionality of LAPACK/ATLAS/Intel> MKL
while hiding the details of interfacing with these libraries.
If> you're interes in seeing what I've got, I'm willing to sht
are.> RyanSure I'm interes. I can contribute too. Also, I've
found it++(google for it), and it seems closest to what I
want. It interfacesblas and lapack, and the syntax resembles
matlab, which is good, but it still implements its own vector
class, so you can only multiply a matrix by their own vector.
With a little hacking, it can be what I need.
====

Subject: Re:
Too many iterations in tqli.cconcerns non- or late termination
in tqli.c of NRC >> we had his already repea time. the code in
NRC is in principle o.k. but they>> use an overly stringent
termination criterion which assumes a small relative error>>
in the elements to be reduced to zero. check this with the
error criteria in>> LAPACK (or the original Wilkinson_Reinsch
volume). adding a small absolute error >> tolerance to the
termination bound removes the problem>> hth>> peter>Exactly
where was this discussed? I'd like to check it out.in the
original source the exit criterion is if abs(e[l])> b then
goto nextitmeans continue iteration . here b is eps (the
matlab eps) times the maximumof
abs(subdiagonal[l])+abs(diagonal element[l]) where the l is
the row numberof the first row in the reduced matrix and the
elements are the transformed onesprior to iteration, hence it
is relative to the value which is to be reducedto zero. the
max iteration count is fixed to thirty. (remark: at that time
typical high precision has been 12 to 14 decimals)in NRC this
reads now dd=fabs(d[l])+fabs(d[l+1]); if { fabs(e[l])+dd == dd
} break; (i.e. accept eigenvalue)d are the diagonal elements !
this is the splitting criterion for thetridiagonal matrix and
the authors of NRC erroneously decided to use this alsoas the
termination criterion. dd is n o t rela to the original
e[l].even worse, an optimizing compiler may translate this
into fabs(e[l]) == 0. n e v e r f o r g e t first to store the
result of theaddition, then recall and compare.maybe the
increased iteration number you observe simply comes from the
higher precision of the modern FPUhthpeter 
====

Subject: New
resourceHi everyone!I have found a very interesting site on
mathematics. It will be usefulfor high school students. You'll
find all theoretical and practicalstuff, and problems solution
variants there.Check it out: www.bymath.com bye!
====

Subject:
Re: Riemann Surfaces - Calculate the integral of 1/(1 -
x^4)^(1/2)> I'm trying to calculate the integral of 1/(1 -
x^4)^(1/2). There are poles at z = +1, -1, +i, -i. Draw a big
circular contouraround the positive half of the z-plane, and
enclose only the pole atz=+i. There is a branch cut on the
positive x-axis, but that is OK. Draw vanishingly small
half-circles around the poles at z = +1, -1. It is easy to
prove the integral vanishes in the far field, lookinglike 1/R
as R->infinity. It is not as easy to show that the
integralsover the half circles around the poles on the real
axis also vanish. I assumed they did. If you evaluate the
2*pi*i times the residue ofthe function at z=+i, then I think
you get pi/sqrt(2) for the answer. I am not too sure of this,
but it's the answer I would put on a HWassignment if it was
due tomorrow.If I remember correctly, the answer can depend on
whether or not youinclude the poles on the real axis in the
contour. Both answerscorrespond to the same problem with
different boundary conditions...orsomething.-John
====

Subject:
Square root of 0,999....Please ignore this if I manage to post
this to the wrong ng! What would the square root of 0,999...
(followed by an infinite numberof nines) be?/Kjell
====

Subject:
Re: Square root of 0,999....> Please ignore this if I manage
to post this to the wrong ng!> What would the square root of
0,999... (followed by an infinite number> of nines) be?Try
sci.math. It often has discussions about why that number is
the sameas 1.00 (one).
====

Subject: gfgfg=====
Subject: seminal
papersI am thinking of going to grad school and study
numerical analysis. I was wondering if one might point me to
seminal (and possibly semi-intelligible to a comp sci senior
in college ;) ) papers written in this area. ALos, any paper
that gives a nice overview of what's going on now would be
nice. Thank you!didier