mm-306 === Subject: : Need help with steepest descent Would someone care to explain how to use the steepest descentalgorithm 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 minimizeRosenbrock's function :-). (I, uh, fell aspeep during the class :-(and now I can't find any references on the net.Thanks in advance,Andrey === 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> Rosenbrock's function :-). (I, uh, fell aspeep during the class :-(> and now I can't find any references on the net.> Thanks in advance,> Andreycheck out mathworld. http://mathworld.wolfram.com/MethodofSteepestDescent.htmlhope 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.( ratherarbitrary amount, you get to decide) by perturb, I mean add orsubtract. it doesn't matter whether you add first or subtract first.let's add incrememnt i to the x1 and x2, then look at what happens tothe 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])/ ithese are numerical derivatives. what happens if the derivative isnegative? this means you are moving downhill towards the local minima.you want to keepp this move.if this is positive than you're moving uphill then you want to makethe opposite move.flank the routine with while loop with conditions to stop it when itreaches a certain minima criteria.good lucksean good luck === Subject: : Re: Need help with steepest descentContent-transfer-encoding: 8bit[[ 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> Rosenbrock's function :-). (I, uh, fell aspeep during the class :-(> and now I can't find any references on the net.For a picture of a descent path down the Rosenbrock banana, see thepicture in the upper left corner of 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 you'llprobably come up with better hits. The basic idea is that from somestarting point near your minimum (or from some arbitrary starting point ifyou have a nice convex function like the above) you repeatedly find thefunction's gradient and choose a new point in the opposite direction.---Roy Stogner === 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 needa rather easy-to-read book).>thanks> 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> peterThank you for your help. Logit model is a model where the dependentvariable takes only one of two possible values - zero or one. Youestimate 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 likelihoodmethod) and just ignore t-test? Can I such an estimation put into myph.d.? === Subject: : Help solving equation.Dear Group,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). zWhere 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 atandrewl at global.net.au I would appreciate it.Best regardsAndrew Lindsay === Subject: : Re: Help solving equation.> Dear Group,> 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).> zI'm 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^3By 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] == 0Cheers,Paul-- Paul Abbott Phone: +61 8 9380 2734School of Physics, M013 Fax: +61 8 9380 1014The University of Western Australia (CRICOS Provider No 00126G) 35 Stirling HighwayCrawley WA 6009 mailto:paul@physics.uwa.edu.au AUSTRALIA http://physics.uwa.edu.au/~paul === Subject: : Re: Help solving equation.> Dear Group,> 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> I'm 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^3That 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]==0Mathematica 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] == 0And that should have read e''[l] + (h/l^2 - b) e[l] == 0Cheers,Paul-- Paul Abbott Phone: +61 8 9380 2734School of Physics, M013 Fax: +61 8 9380 1014The University of Western Australia (CRICOS Provider No 00126G) 35 Stirling HighwayCrawley WA 6009 mailto:paul@physics.uwa.edu.au AUSTRALIA http://physics.uwa.edu.au/~paul === Subject: : Grouping, non-repeatingLet's say you're going to have a dinner party group of 24 people (12couples). 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 Januarycouple A hosts couples B, C, DE hosts F, G, HI hosts J, K, LIn February (and there after) couple A should never be in the samehouse 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 sizesthat must exist in order for the rules to stand up? === Subject: : Re: Grouping, non-repeatingRelated problems are the progressive party problem and the social golfer problem. Google will help you find paperson these problems. These are actually very difficultproblems to solve.-------------------------------------------------------- --------Erwin Kalvelagenerwin@gams.com, http://www.gams.com/~erwin------------------------------------ ----------------------------> Let's say you're 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 velocitiesSorry, there was a typo! I should have squared DeltaTime.This is by no means cartoon physics. People are using this in researchsimulations and I'm wondering if there is a sound basis or not.Thanks!!Free your mind. There is no spoon.************************************************Dr. Patrick Bangerthttp://www.knot-theory.orgResearch Instructor for MathematicsInternational University Bremen> Dear All,> I have heard the folklore that when one wants to simulate some structure> Newton's 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 ofresearch> papers that deal with this issue and perhaps give a proof of this claim asa> theorem. The strange thing is that apparently this is very well known and> widely practised but no one I've spoken to knew where to find the details.> Thank you very much! Best,> Pat> Free your mind. There is no spoon.> ************************************************> Dr. Patrick Bangert> http://www.knot-theory.org> Research Instructor for Mathematics> International University Bremen === Subject: : Re: Computational Evolution without velocitiesX-Enigmail-Version: 0.76.1.0X-Enigmail-Supports: pgp-inline, pgp-mimeThe underlying assumptions allowing to neglect the velocityterm is that it always remains negligible with respect to the forceterm. This situation occurs effectively in very viscous fluids,to random fluctuations. 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 I'm wondering if there is a sound basis or not.> Thanks!!> Free your mind. There is no spoon.> ************************************************> Dr. Patrick Bangert> http://www.knot-theory.org> Research Instructor for Mathematics> International University Bremen>>Dear All,>>I have heard the folklore that when one wants to simulate some structure>>Newton's 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 I've spoken to knew where to find the details.>>Thank you very much! Best,>>Pat>>Free your mind. There is no spoon.>>************************************************>>Dr. Patrick Bangert>>http://www.knot-theory.org>>Research Instructor for Mathematics>>International University Bremen === Subject: : Re: Computational Evolution without velocities> Dear All,> I have heard the folklore that when one wants to simulate some structure> Newton's 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 I've spoken to knew where to find the details.> Thank you very much! Best,> Pat> Free your mind. There is no spoon.> ************************************************> Dr. Patrick Bangert> http://www.knot-theory.org> Research Instructor for Mathematics> International University BremenIf F is constant, the above formula will give a constant incremental distance, i.e a constant velocity. In the absence of drag, that ain't newtonian. Perhaps there was an assumption of a viscous retarding force.-- Joe Legris === Subject: : Re: Computational Evolution without velocities> Dear All,> I have heard the folklore that when one wants to simulate some structure> Newton's laws (i.e. non-relativistic) one can neglect the fact that the> update formula> x_new = x_old - Force DeltaTime / 2 MassYou've got to stop reading those textbooks on cartoon physics.Tom DavidsonRichmond, VA === Subject: : Re: Computational Evolution without velocitiesCc: p.bangert@iu-bremen.de> Dear All,> I have heard the folklore that when one wants to simulate some structure> Newton's 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 can'tsubtract a velocity from a length and get anything meaningful.-- Gordon D. Pusch === Subject: : Re: Computational Evolution without velocities> Dear All,> I have heard the folklore that when one wants to simulate some structure> Newton's 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 can't> subtract a velocity from a length and get anything meaningful.> -- Gordon D. PuschHe obviously meant x_new = x_old - Force * (DeltaTime)^2 / (2 * Mass)However, that does not seem right. Newton's second law is x = F/mIf you wish to eliminate velocities, then in an obvious notation, x(+) = 2x(0) - x(-) + (F/m) * (delta_t)^2will step forward. But you then have to have the values x(0) and x(-).I don't know what algorithm the OP could have been talking about sincedoesn't make sense physically.-- Julian V. NobleProfessor Emeritus of Physicsjvn@lessspamformother.virginia.edu ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. === Subject: : Re: Computational Evolution without velocities* Gordon D. Pusch:>> Newton's 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 relativelysimple and don't contain the velocity. Also, they areused for n-body problems and hard sphere fluids and haveother 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 hisversion, and given that he heard it as folklore, it maywell be it.(Of course, the velocity is not neglected here. It is implicitlycontained in the ODE of 2nd order that is discretized directly.It could also be rewritten into two ODEs of 1st order and discretizedafterwards, which would lead to better known schemes like Runge-Kutta.)-- Rupert === Subject: : Re: Computational Evolution without velocities>* Gordon D. Pusch:> Newton's 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 don't contain the velocity. Also, they are>used for n-body problems and hard sphere fluids 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 threadhere I suspect that some folks who answered were a bit tooquick 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 integrationswhen the velocity is not needed.There are several different Stoermer-Verlet algorithms. Oneis given above by Mazzucco. The original and slightly incorrectformula given by the original poster is half of what is usuallycalled 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 mwhere x is the coordinate, v the velocity, and F the forceevaluated at the given argument. Starting with x_0 and v_0the first equation can be evaluated which allows for thesecond to be evaluated since the needed coordinate is nowknown. So this is NOT a predictor-corrector formula.The set of equations is symplectic and while not of highorder, if the system is conservative and describle by aHamiltonian in which the kinetic energy depends only onthe 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 energyis bounded.Since there is only once function evaluation per cyclethrough the equations (after the first trip) this is afairly cheap integration method. It is, as the originalposter said, quite popular among those who do modellingof galaxies and molecules.There is a large literature on these equations and theirseveral cousins. A search on Verlet and symplecticon Google turned up 528 responses, many of which arequite good. ---- Paul J. Gans === Subject: : Wavelets package in JavaHas anyone ever used a wavelet (de)composition package under Java? Ora program that can be used with Java (looked at Jython, but I didn'tfind a wavelet package for that)?Please mail replies also to gez_75@hotmail.comCheers,Gez === 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 didn't> find a wavelet package for that)?> Please mail replies also to gez_75@hotmail.com> Cheers,> GezThis looks promising: http://www.bearcave.com/software/java/wavelets/. Alsohttp://www.bearcave.com/misl/misl_tech/wavelets/index.html (same site). === Subject: : Coupled system of non-linear PDEHi everybody,I'm trying to solve a non-linear PDE system in u(x,t) e v(x,t) ofthe 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 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 yourefer to in your mail : a Generic B-Spline or a NURBS or a BezierSurface (since you seem torefer 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 coefficientsfor a surface)and therefore:[a] The Coefficient represents a term (of the polynomial function)that maps a point in R(3) space toR(2) (or parametric) space. The point is the control point of theTensor Product Surface for the specific u,vvalue 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 thetensor-product surface are separate concepts and dont mean the same thing.hope this solves the problem.Anup R. Joshiwww.me.mtu.edu/~arjoshi === 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: : WORTH CHECKING OUT!!!!! by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h8MGMjk15740;Has anybody looked at the posting - prime number distribution, solved? Please give me some feedback, of any kind, even criticism could be helpful. If nobody is interested, I'll leave you all alone, but if you are interested, I have more to add. I'm 34 and I've worked on this problem most of my life, so I obviously have a lot invested in it. This is not a crackpot theory, so please take aclose look at it. Thank you.Aaron === Subject: : Revision by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h8MGMo515777;In the previous post, I should have said - any four ODD primes andindeed, any four odd numbers since they all share the samecharacteristic moduli's stated.Please read & respond to prime number distribution - solvedThanks.Synergy === Subject: : new book on Numerical Library in Java by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h8NFuOB11791;A new book with a numerical library in Java was published in August Title: A Numerical Library in Java for Scientists and EngineersAuthor: Hang T. LauISBN: 1584884304publisher: CRC Press email: orders@crcpress.com http://www.crcpress.com/The book contains the source code of a comprehensive numerical libraryin Java. The library can be used to solve a wide range of numericalproblems in linear algebra, the numerical solution of ordinary andpartial differential equations, Fast Fourier Transforms, Time SeriesAnalysis, optimization, parameter estimation and special functions ofmathematical physics. === Subject: : Solving Nonlinear Equation by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h8NDBC900577;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 thesesolvers calculate numerical Jacobian, it perturbs each variable andcalculates jacobian. This makes it very time consuming. Thats why Iused approximate jacobian (It is standard technique in many chemicalproblems) and used Newton raphson method to solve the equations. Now my problem is : Many time I don't have very good initial guess and it is sometimesfar away from the solution. So when i use Jacobian (either approximateor 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..NEVERconverges. 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 o0.2), and in later sateg increase it towarss 0.8, 0.9 when it is closeto convergence criteria. This works and i get solution. But value of w is manual. I am notusing any standard techniqus to get it. So is there any justificationfor this factor? is there any way (fast and reliable) to find thisfactor? Thanks in advance for your helpRegardsPrashant === 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 don't 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? > Thanks in advance for your help don't reenvent the wheel. what you are doing is the damped Newton method with a handcrafted damping. there is a complete and sound theory for thisdeveloped in the seventies of the last century and even downloadable code,the nleq-series of codes from the codelib/elib.(deuflhards 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)))^2where 0 < factor < 1/2 for example factor =0.01 and norm(.)^2 is a scalar product, for example norm=euclidean lengthor 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 solutionif on the set of x's with F-values norm(F(x)) <= norm(F(x(0))) the Jacobianbecomes never singular. for codes see http://plato.la.asu.edu/topics/problems/ zero.htmlhthpeter === Subject: : Re: Solving Nonlinear Equation by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h91Dpvp25153;peter thank you very much for valuable advice. I was looking around forsome information for this. Will apply these theories and see whathappens, thanks againPrashant === 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: : equations (not single equation) by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h8NFuPU11799;hi actually in my previous posting I am using several equations and notsingle equation. Prashant === Subject: : analysis-conneted by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h8NFugr11872;Let p : [a,b] ---> R^3 be a continuous path and a < c < d < b. LetC={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 subsetof the xy-plane in R^3, show that A is still path-conneted. Can youmake 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 don't know for sure, but the way I represented primes in theposting prime number distribution - solved may help you - I've been planning on trying to use it to attack the zeta fcn myself. I've 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 isn't really that hard. Here's a taste of the way it represents primes:105-2^x=y ,(integer x, y<121)y=103,101,97,89,73,41,-23105+2^x=yy=107,109,113for 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 9850GBmodel. 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;Dear AllI am analyzing inter-observer variabilty of a Continuous variable byusing intra-class correlation. Can somebody tell me how to calculatethe 95% confidence interval for this coefficient.Thanksbasu === Subject: : Finding angles using programHow do we find angles in (-180 deg, +180 deg] domain?(If we use dot product, then +45 deg and -45 deg can not bedistinguished)Are there any particular functions in Matlab using which I can findthe angles in the (-180 deg, +180 deg] domain?Thanks,Chetan === 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?>Thanks,>ChetanIt depends. In 2d there is a simple method: for vectors a=(a1,a2) andb=(b1,b2), do the dot product first, then calculate the determinanta1b2 - 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 butthen the angle will be 0 anyway). You then have to adjust your angleaccordingly.(you will notice that this determinant is really the z component ofthe 3d cross-product).In 3d I think you would have to arbitrarily define an orientationrelative to the plane that contains the 2 vectors. === Subject: : Information about Wavelet AnalysisI am looking for information (book name internet source) about Wavelet analysis. === 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: Information about Wavelet AnalysisS. Mallat, A Wavelet Tour of Signal Processing, Academic Pressis said being the bible of wavelet analysis.Axel === Subject: : Re: analysis of a large ODE system. as per Dr. Sspellucci's suggestion.to someone in the forum and perhaps in my school. I went to the library and checked out the book Dr Spellucci hassuggested.( hahn's stability of motion in conjunction with thejacobians and eigenvalues) and guess what?it's been recalled by another person in my school. I wonder who elsefrom my school is reading these forums... I think I'll recall it whenI return it.what are the odds of two people from same school looking for the samebook on the same week right after the book was suggested in anewsgroup.specially when the suggestions for it came after eigenvalues comments.perhaps the book contains some information not found anywhere in theliterature. ( can't imagine what that will be) Though, if you need thebook like me, then you can't be all that great. LOL-----To Dr. Spellucci, that's a great book btw, Dr. Spelluci. (it's very mathematicallyformalized.) perhaps that's why you have suggested it.I have acces,now to some nice introductory materials into nonlineardynamics. I'm sure I would learn a lot based on your comments andthose of others. It really gives me what to focus on when reading thematerials.I really appreciated your comments as well as others. You have alwaysbeen so helpful. I wish people in academics were more like you. (Isound bitter. LOL)respectfully yours,sean from UCIrvine === Subject: : Re: question about 3d plottin software data table limit> Hy, I have a quite noise problem, by a wavelet analysis I've 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 don't 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.> thanks a lot> KioJust remark that a standart computer display has between1 and 1.5 millions of pixels. I think it's hopeless to tryto show 13 millions of informations on it.you can try xd3dhttp://www.cmap.polytechnique.fr/~jouve/xd3d/It will do the job :)-- F.J. === Subject: : A Least Squares Problem charset=iso-8859-1= size(b)-1The problem is that equation (2) has explicit solution (A'*A has inversion) only ifsize(x) = size(b)-1 = rank(A)However, I have to solve the problem in general manner, i.e. I wantto obtain one of optimal solutions. Rank deficiency of A should bepreferredly 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) hasminimal sum(x*x).I plan to use Matlab as proving environment and Fortran/Lapack fordeployment (links to Fortran sources which solve that exact problem are also welcome). I probably have to find QR decomposition of A'A (GEQRF), but I'm lost at the moment what to do after that. I admit my linear algebra is a little rusty...There's an interesting example under Matlab's help on QR function,but it's not directly applicable -- in my case, (A'A) is (can be) singular. All suggestions/solutions are welcome. I suspect this is a well-knownproblem, but it's a little hard to find keywords to google it.-- Jugoslav___________www.geocities.com/jdujic === Subject: : Re: A Least Squares Problem >please don't redirect as I don't read it regularly >I'm 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, U'U=UU'=I V is also unitary and has as much columns as A, abd V'V=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: : Re: A Least Squares Problem charset=iso-8859-1| | | please don't redirect as I don't read it regularlyActually, sci.stat.math or sci.stat.consultwould be just as good.>I'd 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 (don't have any off the top of myhead, and my reference books are several miles away).A look at Rao's book (Linear Statistical Inference, Wiley) might help youin understanding the rank deficient case.Other possibilities:1) Use the SVD to determine the inverse. The usual deal there is to assumethat 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 onproblems of this sort.-- My real email address ismcintosh ##at## research ##dot## telcordia ##dot## com === Subject: : Re: A Least Squares ProblemAn 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 spaceeigenvectors and L a diagonal matrix containing at most r non-zeroeigenvalues of G with r<=p. So ifx = (A'*A)^-1 * A * bAt*A=V*L^2*Ut*U*L*Vt=V*L^2*Vt ,since Ut*U=Iand (At*A)^-1=V*L^-1*Ut and finally, since Vt*V=IThe least square solution is given byx=V*L^-1*Ut*bSVDCMP.f from numerical receips in fortran, and svd(A) from matlab givesimilar results.SVD is the best way to decompose a matrix since it provides extra infomationabout the resolution of the reconstruction (resolution matrices etc...)-----FiLiP Louis---> please don't redirect as I don't read it regularly> I'm 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 won't be> any all-zero rows. As you can see, it may contain two or more linearly> dependent columns (two equal, two mutually negative, or othercombination).> Further, due to nature of the problem (I'd skip the physical background> for now, I can repost if someone's 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 A'A (GEQRF),> but I'm lost at the moment what to do after that. I admit my linear> algebra is a little rusty...> There's an interesting example under Matlab's help on QR function,> but it's not directly applicable -- in my case, (A'A) is (can be)> singular.> All suggestions/solutions are welcome. I suspect this is a well-known> problem, but it's a little hard to find keywords to google it.> Jugoslav> ___________> www.geocities.com/jdujic === Subject: : Re: A Least Squares ProblemGoogle that reference.[38] S. M. Tan and C. Fox, Inverse Problems 453.707 classnotes, TheUniversity of Auckland, Aucland, New Zealand., 2001.What you are doing is referred to as regularized solution, or damped leastsquares.Read the first 3 chapters of the above reference, excellent and very briefreading 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 parameterspace> 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 extrainfomation> about the resolution of the reconstruction (resolution matrices etc...)> ---FiLiP Louis---> please don't redirect as I don't read it regularly> I'm 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 won't 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 (I'd skip the physical background> for now, I can repost if someone's 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 A'A (GEQRF),> but I'm lost at the moment what to do after that. I admit my linear> algebra is a little rusty...> There's an interesting example under Matlab's help on QR function,> but it's not directly applicable -- in my case, (A'A) is (can be)> singular.> All suggestions/solutions are welcome. I suspect this is a well-known> problem, but it's a little hard to find keywords to google it.> --> Jugoslav> ___________> www.geocities.com/jdujic === Subject: : Re: A Least Squares Problem charset=windows-1250|| I'm 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...)Thanks. Please pardon my ignorance, but I still don't see what to doif r | 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...)> Thanks. Please pardon my ignorance, but I still don't see what to do> 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 invertL itself, which is not possible in the rank deficient case, but invert onlyall nonzero elements of L, you get the uniqe x with minimum 2-norm.Greetings Uwe.-- Dr. rer. nat. Uwe Schmitt http://www.procoders.net schmitt@procoders.net A service to open source is a service to mankind. === Subject: : Sinc Interpolation QuestionI 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 andcompute 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 iteasier to understand.How do I go about computing the corrected covariance Se_interp? Should it be positive-definite?Thanks === Subject: : Multivariate Adaptive Regression Splines (MARS)Hi All,I am using MARS 3.5 for scattered data points approximation 2D, 3D, 5D. Whatminimum number of data points is required to start build approximation inMARS?Thank you!-- Best regards,rm -r /Virginia Tech === Subject: : help needed for runaway ODE?Hi All,I am solving a robotic control problem, envolving coupled nonlinear ODEs. Ifound the solution to the ODEs using relaxation method from the NumericReceipe in C for the two-point boundary value problem. But after I changedthe relaxation method to Runge-Kutta method, using the initial values solvedby relaxation method, I got a runaway type of solution with all solutiongrow to very large numbers. I also tried the stiff solver in the NRC, theresults were similar. I know my derivatives are correct in Range-Kuttabecause I computed derivative from the solution of relaxation method andthey matched the derivatives from my Range-Kutta routine. But thederivatives were wrong when I computed it from the solution of theRange-Kutta method. Any suggestions from numeric gurus here? I have beenstruggling on this simple ODE problem for three weeks.Thanks a lot for help.Everettsend your comments to everteq AT sbcglobal DOT net or post it here. Manythanks. === 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.> Thanks a lot for help.> Everett> send your comments to everteq AT sbcglobal DOT net or post it here. Many> thanks. === 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 valuessolved> 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.It's hard to make a substantive suggestion without seeing the exactequations, but it sounds like a job for implicit or semi-implicitmethods.regards, chip === Subject: : Discrete Chebyshev (or other orthogonal) polynomials.Dear all,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.Can anyone help me further, I've already googled on a lot of possible interesting words but nothing came out of it.Thanks === Subject: : Re: Discrete Chebyshev (or other orthogonal) polynomials.> Dear all,> 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.> Can anyone help me further, I've already googled on a lot of possible > interesting words but nothing came out of it.> ThanksRegarding 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 areequidistant.There are also many references regarding applications, for instance inPhysics.For general properties of (continuous) Chebychev polynomials, see[2] Theodore J. RIVLIN , ,,The Chebyshev Polynomials , Wiley -Interscience , 1974 .Regards , Alex === Subject: : Re: Discrete Chebyshev (or other orthogonal) polynomials. >Dear all, >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. >Can anyone help me further, I've already googled on a lot of possible >interesting words but nothing came out of it. >Thanks discrete orthogonal for what scalar product? the chebyshev polynomials ofdegree <= m of the first kind are orthogonal on the discrete set {x_k} = zeros of T_{m+1} with the ordinary scalarproduct. for more information seeChebyshev Polynomialsby 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.95or Pure and Applied Mathematics. New York: John Wiley & Sons, Inc. xvi, 249 p. (1990)hthpeter === Subject: : Re: Discrete Chebyshev (or other orthogonal) polynomials.> Dear all,> 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.> Can anyone help me further, I've already googled on a lot of possible> interesting words but nothing came out of it.> ThanksLet 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. NobleProfessor Emeritus of Physicsjvn@lessspamformother.virginia.edu ^^^^^^^^^^^^^^^^^^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/ === 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 it's now okay to post it on the web?Andrew-- Andrew RossIndustrial and Systems EngineeringLehigh University, www.lehigh.edu/~inime/Bethlehem, Pennsylvania, USAremove all digits <=4 from my e-mail address to reply === 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 it's now okay> to post it on the web?The Dover version might be copyright, but the hardback is a USgovernment publication, and hence (IIRC) not copyright-protected.Cheers,Phil Hobbs === 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 it's 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.> Cheers,> Phil HobbsQuite right. However, keep in mind that by the time you print out this immensebook, you will have spent a lot more than the Dover PB will run you. I don'tknow if the US Printing Office still puts out a hardcover version.-- Julian V. NobleProfessor Emeritus of Physicsjvn@spamfree.virginia.edu ^^^^^^^^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.HI,did u take a look at - Conte/DeBoor: Elementary Numerical Analysis, AnAlgorithmic Approach- ?I found it very useful, many times... IvanoMaurice ha scritto nel messaggio> Dear all,> 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 thesepolynomials.> The description in Numerical Recipes for example also assumes the the> function is known at all locations.> Can anyone help me further, I've already googled on a lot of possible> interesting words but nothing came out of it.> Thanks === Subject: : Help with adjustment to quadratic regression curve algorithmI'm trying to understand/get an algorithm to work which transforms aquadratic regression curve to be unimodal using Bezier control points. Thealgorithm 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 (itsnot very long), but what I'm particularly stuck on is this equation: mY0:=1/m SIGMA (Pi-(Y1)Si = (Y2)(Si)^2 ) i=1where: 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 theparabola with new values of Y2 and Y1 resulting from the first adjustment.The second adjustment, a vertical translation of the B.8ezier-reshapedparabola, occurs only when and Y1 and Y2 have been altered.In fact, I'm not sure that the steps in the algorithm I've done before thisequation are right either:(So ANY help appreciated.PS I have access to matlabTIA === 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 Excel's Solver tool)? Source code implementing this algorithm would be very useful as well. Thanks in advance for any info.-Vik === Subject: : Re: Books On Generalized Reduced Gradient Algorithm?Peter and Erwin, thanks very much for your feedback on this. -Vik === 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 >Excel's Solver tool)? Source code implementing this algorithm would be >very useful as well. >Thanks in advance for any info. >-Vikthe 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.hthpeter === Subject: : Re: Books On Generalized Reduced Gradient Algorithm?I believe this is actually Leon Lasdon's 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 Kalvelagenerwin@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 > Excel's Solver tool)? Source code implementing this algorithm would be > very useful as well. > Thanks in advance for any info.> -Vik === Subject: : skyline matrix solverPlease, I need a skyline matrix solver algorithm. Where can I find one? Thank you,Cpplayer === Subject: : Eignevector differentiationIs it possible to differentiate the eigenvalues and eigenvectors of a squaresymmetric matrix A with respect to any component of A?I.e. if A has i'th eigenvalue lambda_i and i'th eigenvector v_i, what is :d(lambda_i)/d(A_mn)and what is :d(v_i)/d(A_mn)Thanks for any help you can give.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?Yes, if (and only if) the eigenvalue is simple. But the normalizationcondition 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)=1by 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 systemis nonsingular for a simple eigenvalue. The same approach works forother normalizations such as e^Tx(t)=1, where e is a unit vector.Arnold Neumaier === Subject: : Re: Eignevector differentiationX-Enigmail-Version: 0.76.1.0X-Enigmail-Supports: pgp-inline, pgp-mime> 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 i'th eigenvalue lambda_i and i'th 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 differentiatethe eigenvalues, provided they are not multiple. If they are multiple, things gettricky.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 i'th eigenvalue lambda_i and i'th 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 thewhole straight line in its direction: every non-zero multiple ofan eigenvector is again an eigenvector. Normalization does narrow down the choices to two (mutuallyopposite, with no clear decision procedure), which leads to thefollowing implicit equation between A and v (where v' is thetranspose of v): A * v - v * v' * A * v = 0 , v non-zerowith 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 afraidof inverting and integrating to obtain derivatives, here is adifferent philosophy:The projector onto the eigenspace corresponding to an eigenvaluelambda is uniquely determined and has a formula due to Riesz: E = (1/(2*pi*i)) * int[along C] (z*I - A)^(-1) dzwhere the smooth curve C runs once around lambda and misses allother eigenvalues.The rank of E equals the multiplicity of lambda.Differentiation: To get the total differential of E applied toan 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) dzand W(z) = (z*I-A)^(-1) * H * (z*I-A)^(-1)If the curve C is a circle, polar coordinates centered at lambdalead to the integration of a smooth periodic function over fullperiod, 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 toa cluster of eigenvalues, yielding a basis of the invariantsubspace corresponding to that cluster, and also symmetry is notrequired. C can be a sum of simple curves, if convenient (e.g.incase of pairs of complex conjugate eigenvelues).More detailed discussion is found inH. Baumgaertel (or Baumgartel?)Analytic perturbation theory for matrices and operators.Birkhauser Verlag, Basel, 1985.Cheers, ZVK(Slavek). === Subject: : Re: Eignevector differentiationThanks.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 itis:I am investigating the least squares solution to the matrix equation Ax=B.This least squares solution is x=(A'A)^-1A'B. I want to differentiate xwith respect to A, i.e. find dx/dA_mn, but when A'A is close to singular.For this I need to know how (A'A)^-1 changes with respect to A. I alreadyhave this answer when A'A is not close to singular: then you can just used((A'A)^-1)/dA_mn =-(A'A)^-1.(d(A'A)/d(A_mn)).(A'A)^-1However as A'A it is close to singular I should use the Singular ValueDecompostion algorithm to invert A'A. Instead of SVD, for efficiency, I amusing an eigenvector finding algorithm since this gives the same thing, Ithink (i.e. when SVD algorithm is applied to A'A=UWV' then W always turnsout to be the diagonal matrix of eigenvalues of A'A, and U=V is the squarematrix of all the eigenvectors of A'A). That's why I want to differentiatethe eigenvectors and eigenvalues.So once you have the SVD, you can invert A'A easily as VW^-1U' whileneglecting any terms in W that were close to zero. So shouldn't it bepossible to differentiate this SVD inversion, since it was possible when A'Awasn't nearly singular, and the SVD inverion algorithm is clearly defined.Hope that makes things more clear!Thanks again,Mike Fairbank> Is it possible to differentiate the eigenvalues and eigenvectors of asquare> symmetric matrix A with respect to any component of A?> I.e. if A has i'th eigenvalue lambda_i and i'th 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-multipleeigenvalues,> if you plug the characteristic polynomial of the matrix. === Subject: : Re: Eignevector differentiation> Thanks.> Would that give me the derivatives of the eigenvectors as well?the procedure has already been explained by Arnold Neumaier; but beaware that you have to carry out the procedure for eacheigenvector. Since you might be only interested in the dependency ofthe solution vector x on A, you could apply the same method to yourequation. In order to find the dependency of all eigenvalues and eigenvectors onA, you could apply V.I. Arnol'd's Normal Form for matrices, which hasbeen outlined in his book on Geometric methods.> However as A'A it is close to singular I should use the Singular Value> Decompostion algorithm to invert A'A. Instead of SVD, for> efficiency, I amI seriously doubt that that would improve performance! As far as Iknow SVD is much cheaper than finding the eigenvalues andeigenvectors!Another hint: In almost any textbook you will find the recommendationnot to use A'A, but perform the SVD directly to A. Maybe you should also apply bifurcation theory to your problem, sinceat to the singularity you just don't have any continuous dependence. Also the SVD doesn't behave smoothly, when the singular values crossthe threshold.Alois-- Alois Steindl, Tel.: +43 (1) 58801 / 32558 Inst. for Mechanics II, Fax.: +43 (1) 58801 / 32598Vienna University of Technology,A-1040 Wiedner Hauptstr. 8-10 === Subject: : Percentile Estimation in a Monte Carlo SimulationI am trying to estimate a fairly extreme percentile (99.7) for adistribution obtained from Monte Carlo simulation. I am facing aproblem because of the instability of data at the extremes of thedistribution, as in, if the percentile values are calculated fordifferent simulation runs, they seem to jump around a lot. I amcurrently using the MATLAB prctile function to estimate the percentileand am limited to estimate it with 1000 iterations.I would really appreciate if anybody can suggest me a better algorithmthan what is being used in the prctile function, which can provide acloser estimate to the true percentile.Thanks and RegardsUttam === Subject: : Re: Percentile estimation using Monte Carlo by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h93CRHM29151;looks like you are looking for something like a 3-sigma bound. if theends of your distribution are jumping around, you dont have enufpoints and that is it..try more than 1000 points. if you cant do morethan 1000 iterations, do the following:try to predict from the theory of your experiment what kind of pdf(prob density function) you expect. then integrate that so that youcan estimate the 99.97% coverage in terms of the standarddeviation/mean etc (like you do for a gaussuan distribution to get the3-sigma/6-sigma numbers).hope this helps,ankit === 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 thepercentile> and am limited to estimate it with 1000 iterations.> I would really appreciate if anybody can suggest me a betteralgorithm> than what is being used in the prctile function, which can provide a> closer estimate to the true percentile.> Thanks and Regards> Uttam(i) You can work out roughly how many samples you would need to getreasonable accuracy. Treat the quantile as a fixed value and thenumber of samples exceeding this as a random variable ... essentiallyBinominal. The estimate of the probability is a scaled version ofthis, and hence you can get a formula for the standard deviation ofthe estimated probability as a function of sample size. To get areasonable estimate of the probability when it is near 0.003, youmight decide that you should estimate it to within plus or minus0.0001 .... following through with this will show you need a verylarge number of simulations.(ii) To improve matters you will need to move away from straightfowarddirect simulations and make use of some knowledge of the structure ofwhatever you are simulating. A number of general approaches areavailable ... (a) Variance reduction techniques (qv). ( These are more oftensuggested for problems of estimating mean values, but you might finssomething useful under this heading). Examples are the use ofantithetic pseudo-random numbers and using control variables. (b) Importance sampling (qv). You fix up the simulations to generatemore of the extreme events and then downscale the estimated frequency.This can result in much improved estimates of the tails of adistribution for the same number of simulations. (c) Remove some of the random variation from the simulations. Forexample, if your final variable is generated from a distributionconditional on another random variable, you might remove this step(the generation of the final variable) and replace it with atheoretical calculation based on the conditional distribution.David Jones 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 mathematicaTry Mathews & Walker, Mathematical Methods of Physics for anaccessible introduction.-- Julian V. NobleProfessor Emeritus of Physicsjvn@spamfree.virginia.edu ^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. === 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 mathematicaYou cannot solve nonhomogeneous equations with separation of variables. Tryinstead 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: : Fourier transforms (coefficient calculation)...For a few days now, I have tried a number of methods that are supposedto provide the a(k) and b(k) coefficients for a Fourier series. Thesemethods are listed at the end of this post (in Java). However, notone of these methods seems to provide the correct coefficients for thefollowing function; f(x) = 2 - 2 * cos(x)...I never see a '-2' or '2' in the resulting generated arrays. I geta(0) = 4 which is correct, but there is no sign of a(1) anywhere (whenthe expression is changed to f(x) = 2 - 2 * sin(x), things seem towork).Does someone have an algorithm that works for the above? I have triedall the Google searches (the source of the methods below), but to noavail. Please post some example code either in Java or C++ as all the'try searching for ...' suggestions seen in other postings have beenfruitless. In the meanwhile, perhaps the stuff below will be of use tosomeone.Thanks,Matthew=== Subject:=== Subject:=== Subject:=== Subject: Method 1 === Subject:=== Subject:=== Subject:=== Subject:=== Subject:=== Subject:= 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[], doubleai[]) { 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; } }=== Subject:=== Subject:=== Subject:=== Subject:= Methd 2 === Subject:=== Subject:=== Subject:=== Subject:=== Subject:== 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; }=== Subject:=== Subject:=== Subject:=== Subject: Method 3 === Subject:=== Subject:=== Subject:=== Subject:=== Subject:== 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; } } }=== Subject:=== Subject:=== Subject:=== Subject:=== Subject: Method 4 === Subject:=== Subject:=== Subject:=== Subject:=== Subject:= public static void newCompute(int sign, int n, double ar[], doubleai[]) { 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; } }=== Subject:=== Subject:=== Subject:=== Subject:= Method 5 === Subject:=== Subject:=== Subject:=== Subject:== public static void otherCompute(int sign, int n, double ar[], doubleai[]) { 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; } }=== Subject:=== Subject:=== Subject:=== Subject:=== Subject: Method 6 === Subject:=== Subject:=== Subject:=== Subject:=== Subject:=== Subject: public static Ret fourn( double x_re[], double x_im[], int nn[], int ndim, int isign) { Ret ret = new Ret(); float[] data = new float[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] = (float) x_re[i / 2]; data[i + 1] = (float) 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 arraycontaining //the lengths of each dimension (number of complex values), which //MUST all be power of 2. data is a real array of length twicethe //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 arrayby //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; float 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; float 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 = (float) wr * data[k2 - 1] - (float) wi * data[k2]; tempi = (float) wr * data[k2] + (float) 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; }=== Subject:== Class Ret (used above) === Subject:=== Subject:public class Ret{ public double[] x_re; public double[] x_im;} === Subject: : Off Topic in comp.lang.c and comp.lang.c++: 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.This is Off Topicin 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: 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, soI tried (after making it compilable) your method number I with thefollowing 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 Dear all,> 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.> I've 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 I'm 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?> Thanks,> MauriceYou evidently haven't looked at my class notes. Let me STRONGLY suggest youdo--they should answer all your questions, including what to do if thegrid or the weights (in the inner product) are non-uniform. (In that case youcan 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. NobleProfessor Emeritus of Physicsjvn@spamfree.virginia.edu ^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. === Subject: : quadratureCould somebody help me to answer the following question: Let$fcolon[0,infty)to[0,infty)$. What do I have to assume in order toobtain the convergence of the square quadrature for $f$? I'd begrateful for any help. Peter === 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$? I'd be> grateful for any help. PeterI presume you are asking what conditions must f(x) satisfy in order that the(improper) integral --LaTeX format, sorry-- int_0^infty {dxleft| {fleft( 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 offrapidly 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. NobleProfessor Emeritus of Physicsjvn@spamfree.virginia.edu ^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. === Subject: : [newbie] detecting singularity in matricesThe criterion SAS-IML uses for dectecting singular matrices (infunctions meant for solving linear equations) is that, while reducingthe orginal matrix to a triangular/diagonal form, all elements lessthan 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 100chosen?thanks in advance,gc === Subject: : Czy ktos wie jak sie tworzy diagramy krat???PozdrawiamJustyna === Subject: : positive definite ?=> diagonally dominant?Can positive definite results in diagonally dominant? Are there anyexamples if the statement is wrong? Thanks! === Subject: : Re: positive definite ?=> diagonally dominant?> Can positive definite results in diagonally dominant? Are there any> examples if the statement is wrong? Thanks!These are different concepts. For example, is A is diagonally dominant, then (A+b*I) is also diagonally dominant for any scalar b. However, for many values of b, that matrix will not be positive definite. So yes, there are many matrices that are diagonally dominant that are not positive definite.Another way of looking at it is that diagonal dominance is related to quantities like A(i,j)/(A(i,i)-A(j,j)). Shifting the diagonal elements does not change either the numerator or the denominator of that ratio.$.02 -Ron Shepard === Subject: : Re: positive definite ?=> diagonally dominant?> Can positive definite results in diagonally dominant? Are there any> examples if the statement is wrong? Thanks!> These are different concepts. For example, is A is diagonally> dominant, then (A+b*I) is also diagonally dominant for any scalarb. No. As one can read in most textbooks on numerical analysis, diagonal dominancemeans |A_ii|>=sum{k neq i} |A_ik|. Thus your statement is valid only for diagonal matrices.Artnold Neumaier === Subject: : Re: positive definite ?=> diagonally dominant?> No. As one can read in most textbooks on numerical analysis, diagonal > dominance> means |A_ii|>=sum{k neq i} |A_ik|. Thus your statement is > valid only for diagonal matrices.This is relevant to linear equation solutions, not to eigenvalue problems. My post was about the eigenvalue problem, which was the original question. I should have said this explicitly to avoid the confusion.$.02 -Ron Shepard === Subject: : Re: positive definite ?=> diagonally dominant?> No. As one can read in most textbooks on numerical analysis, diagonal> dominance> means |A_ii|>=sum{k neq i} |A_ik|. Thus your statement is> valid only for diagonal matrices.> This is relevant to linear equation solutions, not to eigenvalue> problems. My post was about the eigenvalue problem, which was the> original question. How could you possibly know???The original question by Tong had no context.Apart from that, even books spending a third of their pages toeigenvalues still use diagonally dominant in the sense defined above.See, e.g. Golub and van Loan, Matrix Computations, the bible of Numerical LinearAlgebra. In the eigenvalue context (for the Jacobi method Section 8.5 where thisis relevant) they use careful but different language to denote the way thematrix approaches a diagonal matrix - it does not become more dominant but more diagonal.It pays to preserve the precise usage of established terminology.Arnold Neumaier === Subject: : Re: positive definite ?=> diagonally dominant?> Can positive definite results in diagonally dominant? Are there any> examples if the statement is wrong? Thanks!Linear finite elements on a Poisson problem gives diagonally dominant pdmatrices; higher order elements are still pd but not diagonallydominant.V. === Subject: : positive definite matrixIs 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}|?Thanks. === 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}|?> Thanks.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}|?> Thanks.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 === Subject: : Open source Conjugate GradientFor open-source non-linear conjugate gradient in C,please seehttp://www.willnaylor.com/wnlib.html === Subject: : Re: Closed form determinant of Toeplitz-Hessenberg matrix>I would like to get the determinant of a matrix >typified by (n=4)> 1/2! 1/4! 1/6! 1/8!> -c 1/2! 1/4! 1/6!> -c 1/2! 1/4!> -c 1/2!>where missing entries are zero and c is>a symbolic expression, for arbitrary n.>This is both Toeplitz and Hessenberg (TH).>There is a closed form for dets of Vandermonde >matrices. Is there a similar result for >TH matrices? Thanks.> For n = 4, 5, 6, 7, 8, I tried evaluating the> determinant as a polynomial of degree n-1 in c.> In each case the galois group of this polynomial> was the full symmatric group on n-1 elements,> so there appears not to be a simple formula for its roots.Agree. So far no progress on a closed form. As n->oo the ratiodet(A11)/det(A), with A11 the (1,1) principal minor, appears to approacha Pade form as yet unidentified. === Subject: : 2D FFT for fast polynomial multiplication?I am currently dealing with a problem (details below) where I need to repeatedly multiply twodimensional polynomials (polynomials in two variables). The polynomials are dense (most of the coefficients are != 0) and have approx. the same degree in both variables.I am looking for the fastest way to perform the multiplications and, more importantly, its runtime. In one dimension, a polynomial multiplication can be done in O(n log n) time, where n is the degree (number of coefficients) of the resulting polynomial.My questions are: a) Can two multivariate polynomials (two-variate is sufficient for me) also be multiplied in O(n log n) time, using a multidimensional FFT? Here, n is the number of coefficients of the resulting polynomial, i.e. the product of the degrees in all dimensions. http://www.eptools.com/tn/T0001/PT06.HTM suggests that the 2D FFT can be calculated in O(n log n), so my question probably reduces to whether the multiplication of the transformed polynomials is possible in that time. I guess it is, but I couldn't find any explicit statements on that, and haven't fully understood the theory behind the fast-multiplication-via-FFT. b) As I understand, the FFT and inverse FFT consume O(n log n) runtime and the actual multiplication of the transformed polynomials takes only O(n) time. If I want to perform a lot of multiplications and only need the coefficients of the resulting polynomial at the very end, can this be exploited in order to achieve only O(n) multiplication time? In Cormen's Introduction to Algorithm and various other literature, the fast multiplication is explained via a transformation of the polynomails into a point-value representation. However, this seems to be different from the FFT approach, and I guess that one couldn't be extended to achive O(n) amortised time, because the inverse FFT is necessary to extend the representation to a higher number of (point,value) pairs prior to the multiplication. c) If I have not only two, but k polynomials to be multiplied, what is the time complexity of determining the product? Do I have to subsequently multiply pairs of polynomials, or is there a faster way to work on a bunch of polynomials at once? Can the product of k polynomials in two variables be determined in O(n log n) time, where n is the number of coefficients (product of degrees in the dimensions) in the result?For those who have continued reading till here and are interested in the details, here is the exact function that I need to calculate, just in case it is relevant and the information above doesn't suffice.Given is a (not necessarily binary) tree. For each node, there is a recursive definition of two polynomials p0 and p1: * For a leaf v, p0_v(x,y) = 0 and p1_v(x,y) = 1 or p0_v(x,y) = 1 and p1_v(x,y) = 0 depending on the leaf (there are two classes of leafs, and the input to the algorithm is the tree and the classification of the leaves into the two classes) * For an inner node v, the polynomials p0_v and p1_v are recursively defined throught those for the children w: p0_v(x,y) = product ( p0_w(x,y) + x * p1_w(x,y) ) children w p1_v(x,y) = product ( p1_w(x,y) + y * p0_w(x,y) ) children wNote that the number of edges in the subtree is an upper bound on the degree of the polynomial in both dimensions. As you see, the polynomial is trivial for the leaves and then the operations that need to be performed are: 1.) multiply with x or y 2.) add two polynomials 3.) multiply a set of polynomials (namely those from the >= 2 children)An additional question thus arising might be: d) How could I tackle upper bounds for the total runtime of the algorithm, especially if the tree is unbalanced? Are there any standard papers/estimates on quantising the amount of balancedness (like average depth of a leaf, average size of subtree below a node) in a random tree?However, questions a) to c) are more important at the moment, and d) is probably more of the to solve that is my own business type than the others. Hints to literature are nevertheless appreciated.Any answers to my questions above, pointers to literature or online resources or other hints would be greatly appreciated!Best regards and thanks in advance, Tobias === Subject: : Re: 2D FFT for fast polynomial multiplication?> a) Can two multivariate polynomials (two-variate is sufficient for me)> also be multiplied in O(n log n) time, using a multidimensional FFT?Yes.> b) As I understand, the FFT and inverse FFT consume O(n log n)> runtime and the actual multiplication of the transformed> polynomials takes only O(n) time. If I want to perform a lot of> multiplications and only need the coefficients of the resulting> polynomial at the very end, can this be exploited in order to> achieve only O(n) multiplication time?No.> c) If I have not only two, but k polynomials to be multiplied, what> is the time complexity of determining the product? Do I have to> subsequently multiply pairs of polynomials, or is there a faster> way to work on a bunch of polynomials at once? Can the product> of k polynomials in two variables be determined in O(n log n)> time, where n is the number of coefficients (product of degrees> in the dimensions) in the result?O(k n log n). Yes, there is a faster way. Your third question is redundant. === Subject: : Re: 2D FFT for fast polynomial multiplication?>> c) If I have not only two, but k polynomials to be multiplied, what>> is the time complexity of determining the product? Do I have to>> subsequently multiply pairs of polynomials, or is there a faster>> way to work on a bunch of polynomials at once? Can the product>> of k polynomials in two variables be determined in O(n log n)>> time, where n is the number of coefficients (product of degrees>> in the dimensions) in the result?> O(k n log n). Are you sure about the factor k? Keep in mind that n is the number of coefficients in the *result* - and if there are more factors to start with, the result is also bigger so the k is already contained in n, so to speak.> Yes, there is a faster way.So how do I determine the product of k polynomials faster than by subsequently multiplying pairs of polynomials? Or, which would suffice: which literature could I use as reference to justify a statement of the product of multiple multivariate Polynomials can be determined in O(n log n) time, where n is the size of the result?> Your third question is redundant.Err..the above thing *was* the third question, wasn't it? Why should it be redundant?Best regards, Tobias === Subject: : Re: 2D FFT for fast polynomial multiplication?Message-Id: > b) As I understand, the FFT and inverse FFT consume O(n log n)>> runtime and the actual multiplication of the transformed>> polynomials takes only O(n) time. If I want to perform a lot of>> multiplications and only need the coefficients of the resulting>> polynomial at the very end, can this be exploited in order to>> achieve only O(n) multiplication time?> No.The overall time is still O(n log n) because you still have to do theFFT/IFFT at the beginning and end. However, the time per intermediatemultiplication can be only O(n), if that's what you mean, since you canstay in the Fourier domain for the intermediate operations.Steven === Subject: : Re: 2D FFT for fast polynomial multiplication?> b) As I understand, the FFT and inverse FFT consume O(n log n)> runtime and the actual multiplication of the transformed> polynomials takes only O(n) time. If I want to perform a lot of> multiplications and only need the coefficients of the resulting> polynomial at the very end, can this be exploited in order to> achieve only O(n) multiplication time?>>No.> The overall time is still O(n log n) because you still have to do the> FFT/IFFT at the beginning and end. However, the time per intermediate> multiplication can be only O(n), if that's what you mean, since you can> stay in the Fourier domain for the intermediate operations.Another question about this: how is the fourier-transformed function represented? Is it unnecessary to expand the range of sampling points (which would require an iFFT )? After all, more and more complex functions (polynomials of higher degree) are supposed to require more and more frequency sampling points, aren't they? So if I multiply two fourier-transformed Polynomials of the same m sampling points each, wouldn't I need to evaluate new sampling points, which might require an iFFT?Any hints to literature would be greatly appreciated!Best regards, Tobias === Subject: : Re: 2D FFT for fast polynomial multiplication?> b) As I understand, the FFT and inverse FFT consume O(n log n)> runtime and the actual multiplication of the transformed> polynomials takes only O(n) time. If I want to perform a lot of> multiplications and only need the coefficients of the resulting> polynomial at the very end, can this be exploited in order to> achieve only O(n) multiplication time?>>No.> The overall time is still O(n log n) because you still have to do the> FFT/IFFT at the beginning and end. However, the time per intermediate> multiplication can be only O(n), if that's what you mean,Yes, that's what I meant. Then, the runtime analysis is pretty simple:O(e^2 log e), where e = |E| is the number of Edges in my tree.Best regards, Tobias === Subject: : Counter Intuitive Results: Optimising an FFT routine for SpeedI am confused - I have gone through a two day exercise tuning up anFFT routine, and at the last stage I find that where the biggest gainswere expected, all my conventional understanding of efficient C seemsto be turned on its head.Up until this point I have achieved good speed improvements followingan evolutionary approach, starting from the routine entry point I havemodified code, run, validated and timed a test case and follow thesteepest descent optimisation approach.When I started, with a Cooley-Tukey routine, I was at 2.7ms per 1024point complex FFT (1GHz Celeron - Win2000). I am still using a baby Csystem - Turbo C, but would be surprised if this is the cause of theproblems.I have reduced this time to 0.6ms per FFT call, when I reached theheart of the code, the innermost of 3 loops, where the correlationbetween complex exponentials and data occurs : tempr=wr*data[ixj]-wi*data[ixj+1]; tempi=wr*data[ixj+1]+wi*data[ixj]; data[ixj]=data[ixData]-tempr; data[ixj+1]=data[ixData+1]-tempi; data[ixData] += tempr; data[ixData+1] += tempi;This piece of code is executed 5120 times (ie. N Log_2(N) ) - so a lotof opportunity to improve things. I have removed superfluous arrayindexing and changed to full 64 bit arithmetic (which has proved to befaster than 32bit in previous tweaking) as follows#if defined(FFT_DOUBLE)# define srcPtr0_ doublePtr_1__# define srcPtr1_ doublePtr_2__# define tgtPtr0_ doublePtr_3__# define tgtPtr1_ doublePtr_4__#endif#if defined(FFT_FLOAT)# define srcPtr0_ floatPtr_1__# define srcPtr1_ floatPtr_2__# define tgtPtr0_ floatPtr_3__# define tgtPtr1_ floatPtr_4__#endif srcPtr1_ = srcPtr0_ = &( data[ixj]); ++srcPtr1_; tgtPtr1_ = tgtPtr0_ = &( data[ixData]); ++tgtPtr1_; curReal = wr **srcPtr0_ -wi **srcPtr1_; curImag= wr **srcPtr1_ +wi **srcPtr0_; *srcPtr0_ = *tgtPtr0_ -curReal; *srcPtr1_ = *tgtPtr1_ -curImag; *tgtPtr0_ += curReal; *tgtPtr1_ += curImag;# undef srcPtr0_ # undef srcPtr1_ # undef tgtPtr0_ # undef tgtPtr1_ However now I find the time per call has gone up significantly.Using the FFT_FLOAT option above, the call time increased by 39%compared with the original code.Using the FFT_DOUBLE option, it is a bit quicker, only 20% longer. Soby rights I should just use the original code, with all variableschanged to doubles, but I am loathe to do so just yet until Iunderstand why the current bottleneck exists.Why should accessing vectorised data using pointer references be somuch slower than indexing - surely Turbo-C cannot have pointerde-referencing so badly wrong?Some other ancillary points that I would be grateful for feedback on :I have written a software timer that squeezes better resolution out ofclock() than the CLK_TCK (supposedly 18.2ms), which is normallyavailable.I have timed some runs with a stopwatch, this is what I got_( This is the standard value, ms per each tick in clock())#undef CLK_TCK 18.2 _( This measured over 45s on a Celeron 1GHz, ms per each tick inclock())#define CLK_TCK 56.12 On a 850MHz Athlon (Win-98) system I can get 600ns resolution, howeverwith all the system overheads in Win-2k, even using a 1GHz Celeron, Ionly achieve around 1ms resolution.Is there any easy way to disable or reduce a lot of the Win-2k systemworkload to get more accurate timing? I remember reading about 5years ago of a software stopwatch that could read a far higherresolution CPU timer, anyone have a reference to this?These are typically the sort of results I am gettingFFT_NRC starting to execute 5000 loopsFFT_NRC time = 934456 88(3.5s) nsFFT_NRC_prev starting to execute 5000 loopsFFT_NRC_prev time = 776079 88(3.5s) nsFFT_NRC time change :FFT_NRC_Prev = 158.377 0.124(3.5s) [Micro]s(20.41%)FFT_NRC_prev starting to execute 1000 loopsFFT_NRC_prev time = 777115 451(3.5s) nsFFT_NRC starting to execute 1000 loopsFFT_NRC time = 935104 451(3.5s) nsThanks,Paul. === Subject: : Re: Counter Intuitive Results: Optimising an FFT routine for Speed> When I started, with a Cooley-Tukey routine, I was at 2.7ms per 1024> point complex FFT (1GHz Celeron - Win2000). I am still using a baby C> system - Turbo C, but would be surprised if this is the cause of the> problems.> I have reduced this time to 0.6ms per FFT call, when I reached the> heart of the code, the innermost of 3 loops, where the correlation> between complex exponentials and data occurs :You realize, of course, that your first mistake was starting with the radix-2 Numerical Recipes in C code and spending your time attempting micro-optimizations. The biggest speed gains in an FFT (and most algorithms) come from much higher-level transformations. For example, our FFTW (www.fftw.org) performs the same size-1024 FFT (double-precision) almost 10 times faster than you (0.06-0.07ms on a 1GHz P-III, gcc 3.3).I'm not going to go into the possible transformations that one can do; we have papers on our web site, and better yet you can just go and read the code for some of the faster FFTs (see www.fftw.org/speed for benchmarks) to see their tricks. Such codes aren't as short as Numerical Recipes, I'm afraid, because you trade complexity for performance (as well as generality, e.g. arbitrary sizes). If you just want a fast FFT and are not doing this as a learning experience, you'll be much better off downloading an existing optimized code.Cordially,Steven G. JohnsonPS. With regards to your second question, about timing, you should use a different routine than clock(). The best resolution comes from reading the CPU cycle counter directly (see www.fftw.org/download.html for code to do this on many CPUs).PPS. A decent compiler should transform array references a[j] in a loop into separate incremented pointers for you, if it's advantageous to do so. By being clever you can easily just confuse the compiler's optimizer unless you know what you're doing. === Subject: : Re: Counter Intuitive Results: Optimising an FFT routine for Speed(lots of stuff snipped)> I have reduced this time to 0.6ms per FFT call, when I reached the> heart of the code, the innermost of 3 loops, where the correlation> between complex exponentials and data occurs :> PPS. A decent compiler should transform array references a[j] in a loop> into separate incremented pointers for you, if it's advantageous to do> so. By being clever you can easily just confuse the compiler's> optimizer unless you know what you're doing.I've done the above many a time. The only way I've been able to cope isby learning the instruction set and reading the (equivalent) assembly code.It's pretty easy to tell when you've confused the optimizer.But be aware that modern architectures can execute instructions out of theoriginal sequence and the optimizer is aware of this. So it may refuse toorganize code in the 'best' sequence, knowing that instruction pipeliningwill take care of it. Forcing what appears to be 'best' might goof up therest of the optimization.Most of my experience has been on SPARC, but Pentium architectures sharea lot of thefeatures. === Subject: : Re: Counter Intuitive Results: Optimising an FFT routine for Speed> I am confused - I have gone through a two day exercise tuning up an> FFT routine, and at the last stage I find that where the biggest gains> were expected, all my conventional understanding of efficient C seems> to be turned on its head.> [Replaced array-index operations with explicit pointers, and> the code got slower.]> Why should accessing vectorised data using pointer references be so> much slower than indexing - surely Turbo-C cannot have pointer> de-referencing so badly wrong? You've posted this to quite an array of newsgroups. Here incomp.lang.c (where I've set the followups so as to spare the restof the wide audience), the question is officially off-topic becausethe C Standard makes no promises about the speeds or even relativespeeds of different constructs. On the other hand, this off-topicquestion of efficiency is of great interest, so it winds up gettingwritten about quite a lot ... At the simplest level, your question is number 20.14 in thecomp.lang.c Frequently Asked Questions (FAQ) list http://www.eskimo.com/~scs/C-faq/top.html... but you may be unsatisfied with the answer, which boils downto It depends. For a concrete example of why array indexingmight (*might*) win in your situation, consider this: In codelike your original (much snipped): data[ixj]=data[ixData]-tempr; data[ixj+1]=data[ixData+1]-tempi;... it is fairly easy for the compiler to prove to itself thatdata[ixj] and data[ixj+1] are distinct objects, and data[ixData]and data[ixData+1] are likewise distinct, and that all these aredistinct from tempr and tempi. (A really smart optimizer may evenbe able to prove ixj!=ixData.) If the compiler can determine thatthese data objects do not overlap or otherwise interfere, it canperhaps rearrange the computation for more efficient use of thenumeric pipelines or of the memory access path. Perhaps it canfetch both data[ixj] and data[ixj+1] in one operation, operate onboth of them in registers or other CPU-internal fast locations,and use just one operation to store both results. But when you rewrite with explicit pointers (paraphrased): *out_r = *in_r - tempr; *out_i = *in_i - tempi;... it may not be so easy for the compiler to determine that thepointers are aimed at distinct data objects. And if there's anypossibility that some of these pointers address the same datalocations, the compiler's freedom to rearrange things dwindles.For example, if the compiler thinks it might be possible thatout_r==&tempi, it dares not fetch the value of tempi until afterit's stored the result of the first assignment -- that is, thecompiler may not be able to hold tempr and tempi in registers,and may generate code to fetch them anew every time. Of course, this is speculation: It's only meant to offer aplausible example of the kind of thing that *might* be going on.You'd have to examine the actual generated code to have a hopeof figuring out what's really happening in any particular case.And, of course, the conclusions would only be good for that oneimplementation, and under that one set of circumstances: add aperfectly innocent printf() somewhere and things might change.-- Eric.Sosman@sun.com === Subject: : Re: Counter Intuitive Results: Optimising an FFT routine for SpeedGreetings.[massive post about optimization using Turbo C on Windows snipped]Holy crossposting, Batman! What does this message have to do with Oberon,functional programming, or ISO/ANSI C?For algorithm questions, try an algorithm or programming newsgroup. (Isuppose sci.math.num-analysis would be appropriate for this.)For questions about programming for specific compilers or operating systems,post to a group dedicated to your particular compiler and/or operatingsystem.comp.lang.c is for discussion of standard C only, and not how to tweakspecific compiler or system settings.comp.lang.functional is, I assume, about programming languages such as LISPand Scheme, not procedural languages like C.I have no idea what comp.lang.oberon is doing in the headers.Tristan-- _ _V.-o Tristan Miller [en,(fr,de,ia)] >< Space is limited / |`-' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-= <> In a haiku, so it's hard(7_ http://www.nothingisreal.com/ >< To finish what you === Subject: : 3-d integration routine is neededgoogle but I could not find anything. Do you have any suggestions? Anyreliable routine which involves c/c++/fortran is good enough for me.Cem UZANThanks a lot. === Subject: : Solve Benedict Webb Rubin Equation of state by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h924tx824684;Can any body assist me in solving the BWR equation for all its rootsanalytically.John === Subject: : learning how to use Matlab by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h931Rdg17688;Hello matlab users. I am trying to figure out how to write a Matlabprogram to evaluate the function y=ln(1/1-x)for any user specified value of x, where x is a number <1(note that lnis the natural logarithm to the base e) Use an if structure tovarify that the value passed to the program is legal. If the value ofx is legal, calculate y(x). If not wright a suitable error message andquit.Please give me some pointers on how to write this program or evenbetter if you can solve it email me, it would be better.thanks === Subject: : Code about Testing Self-Similarity by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h94CVFj30133;Hi All!My name is eric i am a student in jerusalem university and i am doing a work about Self-Similarity.can you please send me or give me a reference to source code fortesting self-similarity in Variance-Time Plot,R/S Analysis andPeriodogram Thank all - it will be a real help ! !eric. === Subject: : need the solution by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id h96DQQr28349;question:letwe have sets A and B. A={1,2,3,4,5} (domain) B={7,8,9}(Range)Let a is the member of A and b is the member of B.narate a planthat ensure that the pair (a,b) describes a function. take help of adiagram if necessary.