mm-10 whereA,B are NxN matrices and x is a vector with N elements and^T denotes the transpose.The above equation is true, if A and B are symmetric.Is it possible to relax this symmetry-constraint (at least on B) so that the equation still holds?Frank vector with N elements and >^T denotes the transpose. >The above equation is true, if A and B are symmetric.it is never true. left a column, right a row. you meant (A+B)x=(x^T(A+B))^T >Is it possible to relax this symmetry-constraint (at least >on B) so that the equation still holds? >If dimension is larger than 2, then the nullspace of x has dimesnion >=2. select two linear independent elements y and z in this nullspace and set N=yz^T .then, by choice, Nx=0 and N^Tx=0 .hope this helpspeter Hash: SHA1Online registration form now on our website:http://www.uni-ulm.de/uni/intgruppen/fem/ extension of the FEM Workshop 2002 to an internationalconference turned out to be very successful and also be held in English.The Program of this years workshop is in the web for your informationunder http://www.uni-ulm.de/uni/intgruppen/fem/prog03_ english.htmlDue to the large success of the past workshops, this year we arecelebrating the 10th anniversary of our workshop onThe Finite Element Method in Biomedical Engineering, Biomechanics, andRelated Fields .During the last nine years the workshop has provided an idealpossibility for the exchange of knowledge and experience in differentelds of research. The workshop covers all aspects on computationalmethods (especially the Finite Element Method) applied to Biomechanics,Biomedical Engineering and Biomedicine. A broad spectrum of topics willgive an excellent overview over the present state of the art and shows avariety of application possibilities of the Finite Element Method inBiomechanics. In particular, young scientists will have the opportunityto present and discuss their research results in an open-mindedscientic environment.We cordially invite you to Ulm and would be pleased to welcome you as aparticipant to contribute actively to the discussions of this workshop.Should you be interested in participating in the workshop, wekindly ask you to complete the online registration form on our homepage:http://www.uni-ulm.de/uni/intgruppen/fem/ registration.htmlFor further information regarding workshop topics, registration, traveland accommodation please visit the workshop website:http://www.uni-ulm.de/uni/intgruppen/femWe kindly ask you to pass this invitation on to interested colleagues. What are the classical good methods for 1D interpolation ?I am interested in methods which are not too much expensive to compute, evenif the result of interpolation is of lower quality.Does anybody know if Compactly Supported Radial Basis function [as discussedin papers from Wendland for example] can be used for 1D interpolation.Pierre-Alain. What are the classical good methods for 1D interpolation ?> I am interested in methods which are not too much expensive to compute, even> if the result of interpolation is of lower quality.> Does anybody know if Compactly Supported Radial Basis function [as discussed> in papers from Wendland for example] can be used for 1D interpolation.> Pierre-Alain.Quality of interpolation is not well dened; you can t (1) a LegendrePolynomial up to nth degree and the (polynomial) error will be of order n+1.But in many cases the result will be noisy, especially on the derivatives ofthe t. (2) Spline functions is a classical method for a good stable t.The Spline minimizes the elastic tension of the t and thereby mimic an oldengineering practice of curve tting by bending a wire of some sort. (3) AFourier series (FFT) t gives direct control of the wave content andfrequency cut-off point. Johanneshttp://sizetter.com What are the classical good methods for 1D interpolation ?> I am interested in methods which are not too much expensive to compute, even> if the result of interpolation is of lower quality.> Does anybody know if Compactly Supported Radial Basis function [as discussed> in papers from Wendland for example] can be used for 1D interpolation.> Pierre-Alain.Try Lagrange interpolation. This is discussed in my lecture notes oncomputational methods, downloadable from http://www.phys.virginia.edu/classes/551.jvn.fall01/Look at Chapter 2, Representation of functions.-- Julian V. ^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. Here are two code fragments for calculating the machine epsilonthat give different results on my computer:double temp1, temp2, mchEps;temp1 = 1.0;double mchEps2 = 1.0;// method 1do{ mchEps2 /= 2; cout << mchEps2 = << mchEps2 << endl; }while (1.0 + mchEps2 > 1);// medthod 2do { mchEps = temp1; cout << mchEps = << mchEps << endl ; temp1 /= 2; temp2 = 1.0 + temp1; } while (temp2 > 1.0);These methods give 2 different results and from what I can gather,the second piece of code is doing something to prevent a compileroptimization of some kind from taking place. In fact the 2nd methoddoes yield a result that equals DBL_EPSILON on my machine.Is there some intuitive way of explaining what is going on here? =>> Here are two code fragments for calculating the machine epsilon> that give different results on my computer:>> double temp1, temp2, mchEps;> temp1 = 1.0;> double mchEps2 = 1.0;>> // method 1> do> {> mchEps2 /= 2;> cout << mchEps2 = << mchEps2 << endl;> }> while (1.0 + mchEps2 > 1);>> // medthod 2> do {> mchEps = temp1;> cout << mchEps = << mchEps << endl ;> temp1 /= 2;> temp2 = 1.0 + temp1;> } while (temp2 > 1.0);>> These methods give 2 different results and from what I can gather,> the second piece of code is doing something to prevent a compiler> optimization of some kind from taking place. In fact the 2nd method> does yield a result that equals DBL_EPSILON on my machine.>> Is there some intuitive way of explaining what is going on here?>Optimization stymies the calculation of fp characteristics. So includeoat.h and forget about calculating them.HTH,Gerry T.ps in method 1, change while (1.0 + mchEps2 > 1) to while ((1.0 +mchEps2) > 1) and turn optimization off. =of the optimization that is causing this to happen. I aminterested in this from a theoretical point of view (numericalmethods) more than as a practical matter.Allen>> Here are two code fragments for calculating the machine epsilon> that give different results on my computer:>> double temp1, temp2, mchEps;> temp1 = 1.0;> double mchEps2 = 1.0;>> // method 1> do {> mchEps2 /= 2;> cout \ << mchEps2 = << mchEps2 << endl; }> while (1.0 + mchEps2 > 1);>> // medthod \ 2> do { mchEps = temp1;> cout << mchEps = << mchEps << endl ;> temp1 /= 2;> temp2 = 1.0 + temp1;> } while (temp2 > 1.0);>> These methods give 2 different results and from what I can gather,> the second piece of code is doing something to prevent a compiler> optimization of some kind from taking place. In fact the 2nd method> does yield a result that equals DBL_EPSILON on my machine.>> Is there some intuitive way of explaining what is going on here?>> Optimization stymies the calculation of fp characteristics. So include> oat.h and forget about calculating them.> HTH,> Gerry T.> ps in method 1, change while (1.0 + mchEps2 > 1) to while ((1.0 +> mchEps2) > 1) and turn optimization off. =>> Here are two code fragments for calculating the machine epsilon> that give different results on my computer:>> double temp1, temp2, mchEps;> temp1 = 1.0;> double mchEps2 = 1.0;>> // method 1> do> {> mchEps2 /= 2;> cout << mchEps2 = << mchEps2 << endl;> }> while (1.0 + mchEps2 > 1);>> // medthod 2> do {> mchEps = temp1;> cout << mchEps = << mchEps << endl ;> temp1 /= 2;> temp2 = 1.0 + temp1;> } while (temp2 > 1.0);>> These methods give 2 different results and from what I can gather,> the second piece of code is doing something to prevent a compiler> optimization of some kind from taking place. In fact the 2nd method> does yield a result that equals DBL_EPSILON on my machine.>> Is there some intuitive way of explaining what is going on here?>Optimizations are compiler dependent. But the loops you have constructedterminate under different conditions so it isnt surprising that theresults are different -- even if no optimization is going on.--There are two things you must never attempt to prove: the unprovable --and the obvious.--Democracy: The triumph of popularity over principle.--http://www.crbond.com =>> of the optimization that is causing this to happen. I am> interested in this from a theoretical point of view (numerical> methods) more than as a practical matter.>Take a look athttp://www.stat.unipg.it/stat/statlib/cmlib/ machar.c.netlibThis works with _no_ optimization. Optimization is proprietary to thecompiler.HTH,Gerry T. = Here are two code fragments for calculating the machine epsilon >that give different results on my computer: double temp1, temp2, \ mchEps; >temp1 = 1.0; >double mchEps2 = 1.0; >// method 1 >do >{ > mchEps2 /= 2; > cout << mchEps2 = << mchEps2 << endl; > } >while (1.0 + mchEps2 > 1); >// medthod 2 >do { > mchEps = temp1; > cout << mchEps = << mchEps << endl ; > temp1 /= 2; > temp2 = 1.0 + temp1; >} while (temp2 > 1.0); >These methods give 2 different results and from what I can gather, >the second piece of code is doing something to prevent a compiler >optimization of some kind from taking place. In fact the 2nd method >does yield a result that equals DBL_EPSILON on my machine. > Is there some intuitive way of explaining what is going on here?there are two obvious reasons for this:if oyu work with IEEE754, then internally in the cpu precision is higherthan in stored form in memory. version one of your code allows the compilerto keep 1.0+mchEps2 in register, hence to compute this with increased precision and returning this.(once I experienced an even stronger optimization:the compiler detected the 1 on both sides of the >, eliminated itand the outcome was the smallest positive machine number rather than machineprecision. version 2 is the correct way to do it, since here you force temp2 being stored inmemory, hence rounded to the precision used in memory.hope this helpspeter => Here are two code fragments for calculating the machine epsilon> that give different results on my computer:> ...In addition to the answers you already got, I want to remind you of thevolatile classier for variables. This prevents certain optimizationslike putting them in registers only.Nicolas. =Your explanation was at the level I was looking for, and it>>Here are two code fragments for calculating the machine epsilon>that give different results on my computer:>>double temp1, temp2, mchEps;>temp1 = 1.0;>double mchEps2 = 1.0;>>// method 1>do>{> mchEps2 /= 2;> cout << mchEps2 = << mchEps2 << endl;> }>while (1.0 + mchEps2 > 1);>>// medthod 2>do {> mchEps = temp1;> cout << mchEps = << mchEps << endl ;> temp1 /= 2;> temp2 = 1.0 + temp1; >} while (temp2 > 1.0);>>These methods give 2 different results and from what I can gather,>the second piece of code is doing something to prevent a compiler>optimization of some kind from taking place. In fact the 2nd method>does yield a result that equals DBL_EPSILON on my machine.>> Is there some intuitive way of explaining what is going on here?> there are two obvious reasons for this:> if oyu work with IEEE754, then internally in the cpu precision is higher> than in stored form in memory. version one of your code allows the compiler> to keep 1.0+mchEps2 in register, hence to compute this with increased > precision and returning this.> (once I experienced an even stronger optimization:> the compiler detected the 1 on both sides of the >, eliminated it> and the outcome was the smallest positive machine number rather than machine> precision. > version 2 is the correct way to do it, since here you force temp2 being stored in> memory, hence rounded to the precision used in memory.> hope this helps> peter =I work in the eld of cardiac modelling and my multicellular modelling isdone using a nite difference technique. The stencils that I use are in:1D: |1 -2 1| / h^2 (3-point)2D: |0 1 0| / h^2 (5-point) |1 -4 1| |0 1 0|3D: 7-pointAt the moment, I allow hx, hy and hz, the grid spacings, to be different. Infact, I also allow different grid spacings within a given direction, forinstance (in 2D):+-----+---+| | | hy2+-----+---+| | || | | hy1| | |+-----+---+ hx1 hx2In other words, I have:d2u/dx2 + d2u/dy2 = -2*u(x,y)/(hx1*hx2) -2*u(x,y)/(hy1*hy2) +2*u(x-1,y)/(hx1*(hx1+hx2)) +2*u(x+1,y)/(hx2*(hx1+hx2)) +2*u(x,y-1)/(hy1*(hy1+hy2)) +2*u(x,y+1)/(hy2*(hy1+hy2))Well, this all works ne, but what I would like to do now is to incorporateber orientation. Say that, in 2D, Dx (the diffusion coefcient in the xdirection) is n times bigger than Dy (the diffusion coefcient in the ydirection). Therefore, if the ber orientation is 0 degrees, then diffusionwill occur in an elliptic manner (with the long axis of the ellipse parallelto the x axis). If the ber orientation is at a given angle (different fromzero degrees), the elliptic diffusion should be rotated by that angle andthis is where I am kind of stuck...I believe (please correct me if I am wrong) that there is no way I canimplement ber orientation using a 5- and 7-point nite difference schemein 2D and 3D respectively. I have the feeling that I should use a 9- and27-point stencil instead. I have found a 9-point stencil, which is:|1 4 1||4 -20 4| / (6*h^2)|1 4 1|It works ne, but that assumes hx1 = hx2 = hy1 = hy2, which I would like toavoid if possible... Also, I have been looking for a 27-point stencil, buthave been unlucky so far...So, in summary, I would appreciate if someone could help me with thefollowing:1) Where to nd both a 9- and a 27-point stencil with unequal grid spacingsboth in the different directions and within a same direction, (in the worstcase, I am happy to take anything you may have...)2) How to implement ber orientation in 2D and 3D using a 9- and 27-pointstencil, respectively (or a 5- and 7-point stencil, respectively, in case itis possible)?PS: regarding the implementation of ber orientation, I have given it a tryin 2D using a 9-point stencil, but at the moment it only works ne at 45deg + k*90 deg, with k an integer. At other angles, I still get an ellipse,but the short axis is too long... I guess my formula contains a sillymistake... = >I work in the eld of cardiac modelling and my multicellular modelling is >done using a nite difference technique. The stencils that I use are in: >1D: |1 -2 1| / h^2 (3-point) >2D: |0 1 0| / h^2 (5-point) > |1 -4 1| > |0 1 0| >3D: 7-point >At the moment, I allow hx, hy and hz, the grid spacings, to be different. In >fact, I also allow different grid spacings within a given direction, for >instance (in 2D): +-----+---+ >| | | hy2 \ >+-----+---+ >| | | >| | | hy1 >| | | >+-----+---+ > hx1 hx2 >In other words, I have: >d2u/dx2 + d2u/dy2 = -2*u(x,y)/(hx1*hx2) -2*u(x,y)/(hy1*hy2) > +2*u(x-1,y)/(hx1*(hx1+hx2)) > +2*u(x+1,y)/(hx2*(hx1+hx2)) > +2*u(x,y-1)/(hy1*(hy1+hy2)) > +2*u(x,y+1)/(hy2*(hy1+hy2)) Well, this all \ works ne, but what I would like to do now is to incorporate >ber orientation. Say that, in 2D, Dx (the diffusion coefcient in the x >direction) is n times bigger than Dy (the diffusion coefcient in the y >direction). Therefore, if the ber orientation is 0 degrees, then diffusion >will occur in an elliptic manner (with the long axis of the ellipse parallel >to the x axis). If the ber orientation is at a given angle (different from >zero degrees), the elliptic diffusion should be rotated by that angle and >this is where I am kind of stuck... >I believe (please correct me if I am wrong) that there is no way I can >implement ber orientation using a 5- and 7-point nite difference scheme >in 2D and 3D respectively. I have the feeling that I should use a 9- and >27-point stencil instead. I have found a 9-point stencil, which is: >|1 4 1| >|4 -20 4| / (6*h^2) >|1 4 1| >It works ne, but that assumes hx1 = hx2 = hy1 = hy2, which I would like to >avoid if possible... Also, I have been looking for a 27-point stencil, but >have been unlucky so far... >So, in summary, I would appreciate if someone could help me with the >following: >1) Where to nd both a 9- and a 27-point stencil with unequal grid spacings >both in the different directions and within a same direction, (in the worst >case, I am happy to take anything you may have...) >2) How to implement ber orientation in 2D and 3D using a 9- and 27-point >stencil, respectively (or a 5- and 7-point stencil, respectively, in case it >is possible)? >PS: regarding the implementation of ber orientation, I have given it a try >in 2D using a 9-point stencil, but at the moment it only works ne at 45 >deg + k*90 deg, with k an integer. At other angles, I still get an ellipse, >but the short axis is too long... I guess my formula contains a silly >mistake... > in principle it will be possible, with formulae depending on the orientationand the relative spacings. but not with the reduced number of nodes you imagineand this will be also especially cumbersome with respect to the boundaryconditions. you have to write the taylorexpansion for the grid points with oen grid point as the selected center and now need sufcient equations in orderto eliminate all terms up to one you want to approximate and up to higher ordernegelectable terms:in 2Df(x+c*hx+s*hy,y+c*hy-s*hx)=f(x,y)+f_x(.)*(c*hx+s*hy)+f_y(.) *(-s*hx+c*hy)+ (1/2)*f_xx(.)*(.....) write this for several values of hx and hy, [c s ; -s c] being the rotation matrix for the ber direction. now combine these equtions linearly in orderto eliminate f, f_x, f_y, f_xy, f_yy (in order to get f_xx) etc. why not using nite elements where you can achieve all this quite easily even with standard software?(there exist nite difference formulae for hexagonal grids which may also apply here, but those published use xed orientation)hope this helpspeter spellucci@mathematik.tu-darmstadt.de>>the normal 2d polynomial t model is >> z = model(x,y,a) = a00 + a10*x+a01*y +a20*x^2+a11*x*y + .........>>given data (x(i),y(i),z(i)) you wnat to compute a00, a10, a01, a20,..> such that the >sum of squares distances z(i)-model(x(i),y(i),a) is minimal. >of course, the parameters will not appear as an array with two indices but >simply as a vector and the meaning of the coefcient is connected with >the powers of x and y you store in the corresponding matrix column.>>you can do this (but should not) like in the source you mention,>namely explicitly computing and inverting the normal equations matrix.>better use some decent linear least squares code like dgelss from lapack>or equivalents, see e.g.>http://plato.la.asu.edu/topics/problems/nlolsq.html.>but well. this was not your question. >>scaling is indeed necessary After thinking about it for awhile, my curiosity has gotten the bestof me. Why do you say you can do this (but should not)?I understand the merits of using a high quality linear least squarescode (or any high quality code for that matter).Our inhouse linear least squares solver evolved to C++ from our legacyfortran libraries (that probably evolved from something lapack-ish).It has worked quite well for us and has had 25+ years worth of GISdata from all over the world put through it. The problem we are having only occurs when using a mercator projectionwhich has a very large x and y range. Mercator coordinates are quitelarge (numbers with 6-7 digits to the left of the decimal point).This same problem can be seen for example with the calculation of thepolygon centroid. Given the following 4 2D (mercator projection)points (that form a closed triangle)1 - ( 302083.31500000, 5045133.0490000 )2 - ( 302085.53600000, 5045129.9390000 )3 - ( 302081.31600000, 5045129.9390000 )4 - ( 302083.31500000, 5045133.0490000 )If one calculates the centroid the normal way, the centroid:( 302082.87052229, 5045122.3160023 )is not inside the triangle.The solution here was to normalize (with a shift only) the data pointsand to apply the reverse shift to the centoid point.We knew that this problem existed with the polynomial t for sometime and have been able to get around it because the le format ofthe transformed data stores a transformation origin. We use just ashift for this (and hence my original question about shift).Now we want to write the transformed data out to another similar butnot identical format that does not store the transformation origin(but does store the original control points and the coefcients).By scaling to the mean (of the control points), calculating andunscaling the coefcients, we do not need to store the origin (or themean).Which brings me back to my original question. Why do you say you cando this (but should not)? Is there something I dont know about?Todd spellucci@mathematik.tu-darmstadt.de >>>the normal 2d polynomial t model is >>> z = model(x,y,a) = a00 + a10*x+a01*y +a20*x^2+a11*x*y + ......... >>>given data (x(i),y(i),z(i)) you wnat to compute a00, a10, a01, a20,.. >> such that the >>sum of squares distances z(i)-model(x(i),y(i),a) is minimal. >>of course, the parameters will not appear as an array with two indices but >>simply as a vector and the meaning of the coefcient is connected with >>the powers of x and y you store in the corresponding matrix column. >>>you can do this (but should not) like in the source you mention, >>namely explicitly computing and inverting the normal equations matrix. >>better use some decent linear least squares code like dgelss from lapack >>or equivalents, see e.g. >>http://plato.la.asu.edu/topics/problems/nlolsq.html. >>but well. this was not your question. >>>scaling is indeed necessary >After thinking about it for awhile, my curiosity has gotten the best >of me. Why do you say you can do this (but should not)? >I understand the merits of using a high quality linear least squares >code (or any high quality code for that matter). >Our inhouse linear least squares solver evolved to C++ from our legacy >fortran libraries (that probably evolved from something lapack-ish). >It has worked quite well for us and has had 25+ years worth of GIS >data from all over the world put through it. >The problem we are having only occurs when using a mercator projection >which has a very large x and y range. Mercator coordinates are quite >large (numbers with 6-7 digits to the left of the decimal point). >This same problem can be seen for example with the calculation of the >polygon centroid. Given the following 4 2D (mercator projection) >points (that form a closed triangle) >1 - ( 302083.31500000, 5045133.0490000 ) >2 - ( 302085.53600000, 5045129.9390000 ) > 3 - ( 302081.31600000, 5045129.9390000 ) >4 - ( 302083.31500000, 5045133.0490000 ) >If one calculates the centroid the normal way, the centroid: >( 302082.87052229, 5045122.3160023 ) >is not inside the triangle. >The solution here was to normalize (with a shift only) the data points >and to apply the reverse shift to the centoid point. >We knew that this problem existed with the polynomial t for some >time and have been able to get around it because the le format of >the transformed data stores a transformation origin. We use just a >shift for this (and hence my original question about shift). >Now we want to write the transformed data out to another similar but >not identical format that does not store the transformation origin >(but does store the original control points and the coefcients). >By scaling to the mean (of the control points), calculating and >unscaling the coefcients, we do not need to store the origin (or the >mean). >Which brings me back to my original question. Why do you say you can >do this (but should not)? Is there something I dont know about? >Todd,you gave an URL, I looked there and detected there was used an explicitmatrix inversion for the normal equations. this is clearly not the stateof the art. hence I assumed that you also intend to do this. ifyou use an inhouse lsq-solver which evolved from somewhat lapack like, this iscompletely o.k. I see your point concerning the representation of the dataand the solution I proposed will work, but is not optimal, since you then work on[0,n] times [0,m] for (n,m) grid (in the worst case), not on [-1,1] times [-1,1]which would be nicer. nevertheless the situation is much improved and ifall the x_i resp y_j are in the same order magnitude you are working approximately on [0,1] times [0,1].hence this should be o.k.peter =xd3d is a simple scientic visualization tool designed to be easy to learn.It can plot 2d and 3d meshes, with shadowing, contour plots, vector elds,iso-contour (3d), as well as 3d surfaces z=f(x,y) dened by an algebraicexpression or a cloud of points.It generates high quality vector PostScript les for scientic publicationsand still or animated bitmap images.It includes the graph plotter xgraphic.Have a look athttp://www.cmap.polytechnique.fr/~jouve/xd3d/-- Francois Jouve (Francois.Jouve@polytechnique.fr)Centre de Math. Appliquees, Ecole Polytechnique, 91128 Palaiseau (France)http://www.cmap.polytechnique.fr/~jouvehttp:// www.cmap.polytechnique.fr/~optopo =I am very pleased to inform you that we have just released SansGUI version1.2, an update to the modeling and simulation environment for developingand deploying scientic and engineering simulation programs withoutwriting any graphical user interface (GUI) code. The new version supportsindustry standard OpenGL programming with Compaq Visual Fortran andMicrosoft Visual C++ on Windows platforms for animated 3D graphics anduser controls.As you all know, laying out and coding GUI front and back ends for yourcomputational programs is a major hurdle in your development process. Youmay use the programming language facility to call low level APIs directly,or some intermediate layers of APIs from a 3rd party. No matter what youdo, you still need to handle low level events, such as window creation andmaintenance, data memory and GUI contents synchronization, mouse buttonclicks, printing, undo/redo operations, etc. When you have the underlyingspecications changed, it is very tedious to synchronize such changeswith your existing data les, your data structures in the memory, and theGUI layouts. With the SansGUI environment, you can combine an interactiveschematic editor for hierarchical model building, a graphical userinterface for data entry and validation, a dynamic charting facility andOpenGL animated 3D graphics programming support for data exploration andvisualization; all of these without the need to write any GUI code. Aschema version evolution facility is built-in for synchronizing thechanges of your data specications with the users data les and the GUIautomatically, with both forward and backward compatibility.How can all these be done without low level GUI API calls, you may ask?Weve done it by implementing a meta-schema that lets you dene theschemas of your classes (model building blocks) for data objects. Theseclass schemas provide essential information for data objectsynchronization among different versions. The data objects are packagedin a universal data object format and passed around to your own code viaDLL functions, which has only a single API function prototype. In fact,we also use these data objects in many tasks inside SansGUI. The SansGUIdata object format and the single API function prototype are the onlythings you need to learn in order to program SansGUI, not hundred orthousands of API functions.The OpenGL 3D graphics support is done by introducing a Graphics class,which uses exactly the same data object format. The Graphics classprovides you with an application framework that manages multiple 3Dgraphic windows and interactive user interface features, such as 3D objecttranslation, rotation, zooming, selection, continuous run, single step orfast forward animation controls, printing, exporting image les, and manymore. Your code simply fetches the 3D operational data, along with thedata of your custom elds, from the data object and calls OpenGL routinesdirectly to render the graphical scene.For legacy code integration, SansGUI provides an external process controlstand-alone program to run in a separate process via a customizablecommand script le.For more details, please visit:http://protodesign-inc.com/sansgui.htm - for the main product pagehttp://protodesign-inc.com/sg_exStart.htm#solid - for a new 3D graphicsexample with full source code in CVF and MSVC, implemented independentlyhttp://protodesign-inc.com/les/SGtutSho.pps - a new Getting Startedwith SansGUI slide show in PowerPoint format (60 slides)http://protodesign-inc.com/les/SGtutSho.pdf - same slide show but in PDFformatinformation.-- Greg Chienhttp://protodesign-inc.com note to inform you that you can now order your copy online andmake the payment by credit card.For further information please go to:http://www.gene-expression-programming.com/gep/Books/ index.aspCandida Ferreira----------------------------------------------------- ----------Candida Ferreira, Ph.D.Chief Scientist, Gepsoft73 Elmtree DriveBristol BS13 8NA, UKph: +44 (0) 117 330 9272www.gepsoft.com/ gepsoftwww.gene-expression-programming.com/author.asp-------- ------------------------------------------------------- = 100 miles far from the truth: Iknow James Harris, and he is the most brilliant mathematician ever.To tell the truth... well... Im not sure James will appreciate... butyou can check the following link, where youll nd:1) The social status of JSH.Now you can meet the myth ...http://www.jal.cc.il.us/math/faculty/jim.html mathematicians,> old friend JSH; the fact is youre all 100 miles far from the truth: I> know James Harris, and he is the most brilliant mathematician ever.> To tell the truth... well... Im not sure James will appreciate... but> you can check the following link, where youll nd:> 1) The social status of JSH.> Now you can meet the myth ...> http://www.jal.cc.il.us/math/faculty/jim.htmlwhat about this guy? http://www-ee.stanford.edu/~harris/ mathematicians,>> old friend JSH; the fact is youre all 100 miles far from the truth: I> know James Harris, and he is the most brilliant mathematician ever.> To tell the truth... well... Im not sure James will appreciate... but> you can check the following link, where youll nd:>> 1) The social status of JSH.>> Now you can meet the myth ...>> http://www.jal.cc.il.us/math/faculty/jim.html> what about this guy? http://www-ee.stanford.edu/~harris/Who? The panda, or the one holding the panda?GC =According to my own tests (pentadiagonal A, Ax=B) UMFPACK requiresslightly more time than old Y12M direct solver(Zlatev). The precision is comparable. The CG approach is fasterbut one component of solution (x) is always false.M.>>I need to solve a sparse matrix many times. (A*x=b, b varies)> Then UMFPACK would let you save the factorization and apply that many> times. Since the factorisation is the expensive part, doing this is a> good idea.> Its much harder to reuse information in a CG approach.> V. =Ch NAG Statistics Package is released. Ch is a C/C++ interpreterfor numerical computing and 2D/3D plotting.Ch NAG Statistics package is for rapid statistical applicationdevelopment, it makes a set of the robust NAG statistical routinesavailable to users of Ch.With SoftIntegrations free Ch ODBC and Ch CGI toolkit, it is easy tostatisticalanalysis with graphical presentation.Ch NAG Statistics Package makes the rigorously tested statisticalroutinesin the NAG C library available to Ch users.Key FeaturesThe Ch NAG Statistics Package features routines for:1. Correlation and regression, multivariate methods, and analysis of variance;2. Random number generation;3. Nonparametric statistics, smoothing, contingency table analysis, and survival analysis;4. Time series analysis;5. Interactive execution of routines for rapid application developmentand deployment.More about Ch NAG Statistics Package can be found at http://www.softintegration.com/products/package/statistics/== ==We have just released our program VisuMap 1.4. The standard edition of theprogram is free for non-commercial use. It can be downloaded at:http://www.visumap.net/Products.htmlVisuMap is a visualization tool for high dimensional data. It basicallyturns any numerical table into a 2D map that preserves, as much as possible,the distance information of the original dataset. VisuMap is based on a newdimensionality reducing method called relational perspective map (RPM). RPMworks like traditional multidimensional scaling (MDS) methods, but showsquite different characteristics.Additionally VisuMap also implements other MDS methods like Sammon mapping,curvilinear component analysise, principal component analysis as well asclustering methods like K-mean clustering and self-organizing map (i.e.Kohonen net).The following link is an example that depicts DOW 30 index stocks. Each spoton the map represents a stock. Closely located stocks have similar weeklyprice history in the year 2002. The price histograms of some stocks aredepicted as examples.http://www.visumap.net/Applications/ AppFinancial.htmlNew features and major changes in VisuMap 1.4 optimization algorithm for RPM map that is about10 times faster than previous simulated annealing based algorithm Implemented two new mapping algorithms: Sammon mapping and thecurvilinear component analysis. Many enhancements in user interfaces: congurable toolbars,drag-and-drop supports, customizable glyphs etc. Bug xes. Upgraded the platform from .NET 1.0 to .NET 1.1.James X. Lijaemsxli@visumap.net =I am now using ARPACK to solve a generalized eigenvalue problem. I nevergot the right results. There must be something wrong in my understandingof the code.I just modied the dsdrv3.f under the examples/sym/For test, I just test on two jordan canonic matricesa=[2 2 0 ] m=[3 2 ] [2 3 ] [2 3 ] [ 4 2 ] [ 3 2 ] [ 2 5 ] [ 2 3 ] [ 6 2 ] [ 3 2 ] [ 2 7 ] [ 2 3 ] [ ... ] [ ... ] [ 150 2 ] [ 3 2] [ 2 151 ] [ 2 3]They are both p.d. and symmetric, what I did are1) replacing the 2 call av and the 2 call mv (one is mainloop after dsaupd, the other is in computing the residue error)with theworkd(n+1:2*n)=matmul(a,workd(1:n))workd(n+1:2*n)=matmul(m ,workd(1:n)),ax=matmul(a,v(1:n,j)); mx=matmul(m, v(1:n,j))2) when ido=1 or ido =-1, replace the dgttrs with dgesv ( the Lapacksubroutine for solve generalized matrix), workds on m to getinv[m]*A*xThe results are incorrect, even though no warning, no error info. Ritz values and relative residuals ---------------------------------- Col 1 Col 2 Row 1: 8.58771D+03 8.25552D-01 Row 2: 1.91678D+04 1.44619D+02 Row 3: 2.91610D+06 1.79778D+03 Row 4: 6.24538D+10 2.63107D+05 Row 5: 2.83229D+19 5.60303D+09 Row 6: 5.80098D+36 2.53574D+18 _SDRV3 Size of the matrix is 150 The number of Ritz values requested is 6 The number of Arnoldi vectors generated (NCV) is 15 What portion of the spectrum: LA The number of converged Ritz values is 6 The number of Implicit Arnoldi update iterations taken is 18 The number of OP*x is 131 The convergence criterion is 1.1102230246251565E-16While with the Lapack routine, dsygv, I did get right answers. Cansomebody help me out of this puzzle? You cannot imagine how I appreciate it.I spent 3 weeks on the sparse generalized eigenvalue problem, from LANZ toARPACK. I am really upset that I cannot make them run.Best Rgds. =The problem:I have a relation R wheredomain set is: dom R = {A,B,C}, andrange set is: ran R = {1,2,3}.The relations (maplets) are:{A>1,A>3,B>2,B>3,C>1,C>3}The following sets would be performed, if we use at most once every elementfrom each side of the relationship, within each new set. (E.g. {A>2,B>1} isok but {A>2,A>3,C>3} is not ok because A and 3 is being twice in the set).We have:rst {A>1,B>2,C>3}second {A>1,B>3} (not C3 because 3 is being used once in the set)third {A>1,C>3}fourth {A>3,B>2,C>1}etc.I try to capture the sets with the most number of elements (rst and fourthline, that have 3 elements each) and be able to store these somehow, forfurther manipulation.Here is an example (in case I didnt explain the problem well) :the manager A hasnt got the qualications to manage department 2, thus B>1and C>2 are not applicable as well.Each manager could manage only one department and each department is managedby only one manager.I try to nd a way to occupy as many managers as I can with as manydepartment possible (highest number of matches).I have spent long time with this problem, and I looked into the jobassignment problem in operations research and in maths, but still, I couldnt create any algorithm or code (vb or java) .Any help is appreciated.Jim =Suppose a positive random varialbe X with pdf f_X(x) and 0 Given a value A between x1 and x2,i.e. x1 into two pars, [x1, A] and [A, x2]. I need to nd the average of X in the> range [x1,A], and additionly the average value of X in the range [A, x2].> how to write the expression?Assuming that f is continuous and P(a<=X<=b)>0 b Integral x*f(x) dx a E(X|a<=X<=b) = ------------------ b Integral f(x) dx a Horst =Maybe I have not stated the question clearly.A positive random variable X varies in the range [a,b] with pdf f_X(x).Suppose a real number r greater than a and less than b, i.e. a> Maybe I have not stated the question clearly.> A positive random variable X varies in the range [a,b] with pdf f_X(x).> Suppose a real number r greater than a and less than b, i.e. a of the random variable will fall in the range [a,r] or the range [r,b].> Dene the event> A = value of random variable falls in the range [a,r];> B = value of random variable falls in the range [r,b];> Then, how to compute the expected value of X with the value in the range> [a,r] and in the range [r,b]. namely, to nd> E(X|A) and E(X|B)You stated your question clearly and Im going to repeat my answer A Integral x*f_X(x) dx aE(X|A) = -------------------- A Integral f_X(x) dx a b Integral x*f_X(x) dx AE(X|B) = ------------------ b Integral f_X(x) dx Aassuming that f_X ist continuous and Pr(A) and Pr(B) are not 0.Horst =You must decide what it is you want to average. It is certainly _NOT_number of cars so dividing by zerois not a problem. It _COULD_ be average price of a car. In thatcase, the sign of the cash owdoesnt matter, so the average would be:(140+170+180+130+140+150)/6 = 910/6 = 455/3 = $151.67On the other hand, it could be (and I agree with David Wilson, thiswould make more sense,Net prot per sell/buy exchange, in which case your average would be((140-130)+(170-140)+(180-150))/3or (10+30+30)/3 = 70/3 = $23.33Or you could choose to average the net prot per trade. Since therewere 6 trades, thatresult would be70/6 = $11.67Dividing by 4 is denitely wrong. Whoever it was who right enough, but they led you off down a side road.Its a not to count them in the average. If I play 100 games at the $1 slots and get one $10 payoff, my average prot is (10-100)/100 =-$.90, not +$10.Jack> I need to calculate a net average around 0.> I have sold 3 cars and bought 3 new cars. My net total of cars is 0.> if I sold the cars for:> 140> 170> 180> and bought the new ones for:> -130> -140> -150> what is my net average for these amounts?> I cant divide by 0 so how can I calculate the answer?> I have been be:> (170+180-130-150)/3 = 70/4> but that seems wrong!> Surely the average takes no account of the net no. of sales and the answer is:> (140+170+180-130-140-150)/6 = 70/6> Is there a formula for this relating to the net no. =0??> Michelle =The Kansas State University Laboratory for Knowledge Discovery inDatabases is pleased to announce the initial pre-alpha of BayesianNetwork Tools in Java (BNJ) v.2.BNJ is an open-source suite of Java tools for research andapplications development in probabilistic reasoning, and currently includes the following modules:- Exact inference (junction tree, variable elimination)- Sampling-based inference (logic sampling, importance sampling)- K2 algorithm for structure learning- Experimental packages for structure learning (stochastic score-based structure search, variable ordering)- Format interconversion suiteA GUI-driven editor and front end is currently in redevelopmentand will be available with the full release of BNJ v.2.The latest stable release of BNJ is always available on Sourceforgeand can be downloaded with or without Javadoc.BNJ is compatible with the XML Bayesian Network interchange formatand can convert among many network formats, including commercial(Hugin, Netica, Microsoft BN) and academic ones (U. Pittsburgh Genie,U. Sao Paulo Bayesian Interchange Format, Hebrew University LibB).For more information, please see the BNJ distribution page:http://bndev.sourceforge.netand the BNJ development group:http://groups.yahoo.com/group/bndevBNJ is distributed under the GNU General Public License.-William Hsu-- William H. Hsu, Ph.D. Assistant Professor of CIS, Kansas State University Director, Lab for Knowledge Discovery in Databases bhsu-AT-cis.ksu.edu, bhsu-AT-ncsa.uiuc.edu SeriesPossible Real World System Constructshttp://tonylance.mybravenet.com/cricket.html Access page to 138K ZIP leAstrophysics net ring Access siteNewsgroup Reviews including alt.support.attn-decitcomplete with programs, source code in listing format and documentation.Pastures Software Package.Sub-atomic Mesons, Baryons and Leptons Classication System.(C) Copyright Tony Lance 1997Distribute complete and free of charge to comply.Big Bertha Thing cricket The Rules of Cricketby Tony Lance 1. 1st side goes in till all out. (10 out of 11) 2nd side goes in till all out. 1st side wins.Exceptions;-i Rain stops play is draw.ii 2nd side scores more, game stops they win.iii Equal nal scores are a draw. 2. Every run counts one.Exceptions;-i A no runs boundary hit counts 4.ii A no runs clean over boundary hit counts 6.iii A no hit too wide counts 1.iv A no hit boundary counts 4.v A no hit over boundary counts 6. (missing rule)vi A no hit can still be run. 3. Caught out or bowled out. (owzat!)Exceptions;-i Stumped out.ii Run out.iii Thrown out.iv Trod on wicket.v Retired injured. Exception to all of above.Runner for slightly injured. Street Cricket Forget all the exceptions and rules 1 and 2.Equipment;-i Cut out bat.ii Rubber ball.iii Piece of chalk. Owner of bat goes rst.Owner of ball goes rst.Fielder bowls next.Bowled out, bowler bats.Caught out, catcher bats and swaps bat for ball.The end(C) Copyright Tony Lance 2000Distribute complete and free of charge to comply.Tony Lancepobox47@big-bertha-thing.comBig Bertha Thing debitUK big banks have now made debit cards mandatory.1. Every cheque book has a debit card.2. You do not sign agreement for Insure debit cards for 1000 pounds sterling.5. For every 8 pounds loss charged to the customer, the bank can lend 92 pounds. (Multiplier)6. Cheque guarantee card included.7. Cash Card included.8. Phantom withdrawals included.9. 50 pounds cashback at supermarkets.10. 250 pounds Hole-in-the-wall withdrawals.11. 100 pounds Post Ofce and bank withdrawals.12. How do you prove that you have cut the card up and disposed of it?13. How do you prove that you did not use it before destruction?14. Cut the card up and give to your solicitor or third party for safe-keeping. (Escrow)15. The survival of the banks depends on putting customers to the sword. (Barbaric) Obviously not all customers, just enough.16. Unauthorized access to debit card accounts, includes both criminals and bank insiders, which seem to be synonymous.17. Authorized access to debit card accounts, includes bank insiders, grazing on the customers like milk cows, in their ofcial capacity, as per job description.18. You no longer need to prove this in a court of law, since it is endemic and ubiquitous. (Common knowledge) Case law sufcient for class action. =In the biharmonic problem:- boundary conditions on function and rst derivative are Dirichlet(essential or kinematic)- boundary conditions on second and thid derivative are Von Neumann (naturalor static)is that correct?Antonio => prurposes ?Grace Wahba has published a lot on smoothing splines (i.e., polynomialsplines with a penalty term that encourages smoothness). Her theoreticalthere are many papers of hers on the web, including some introductorypapers that talk about reproducing kernels. Try looking for her home page, or see also http://www.researchindex.com/.For what its worth,Robert Dodier--``If I have not seen as far as others, it is because giants were standing on my shoulders. -- Hal Abelson => prurposes ?> Grace Wahba has published a lot on smoothing splines (i.e., polynomial> splines with a penalty term that encourages smoothness). Her theoretical> there are many papers of hers on the web, including some introductory> papers that talk about reproducing kernels. > Try looking for her home page, or see also http://www.researchindex.com/.> For what its worth,> Robert Dodierwhy splines ? Whar about collcation methods ?Mathematical Methods of Physics and Engineering with MathematicaFerdinand F . Cap, CRC-Press/Chapman and Hall, www.crcpress.com1 Introduction1.1 What is a boundary problem?1.2 Classication of partial differential equations1.3 Types of boundary conditions and the collocation method1.4 Differential equations as models for nature2 Boundary problems of ordinary differential equations2.1 Linear differential equations2.2 Solving linear differential equations2.3 Differential equations of physics and engineering2.4 Boundary value problems and eigenvalues2.5 Boundary value problem as initial value problem2.6 Nonlinear ordinary differential equations2.7 Solutions of nonlinear differential equations3 Partial differential equations3.1 Coordinate systems and separability 3.2 Other methods to reduce partial to ordinary differential equations3.3 The method of characteristics3.4 Nonlinear partial differential equations4 Boundary problems with one closed boundary dened by coordinate lines4.1 Laplace and Poisson equation4.2 Conformal mapping in two and three dimensions4.3 DAlembert wave equation and string vibrations4.4 Helmholtz equation and membrane vibrations4.5 Rods and the plate equation4.6 Approximation methods4.7 Variational calculus4.8 Collocation methods5 Boundary problems with two closed boundaries5.1 Inseparable problems5.2 Holes in the domain.Two boundaries belonging to differentcoordinate systems5.3 Corners in the boundary6 Nonlinear boundary problems6.1 Some denitions and examples6.2 Moving and free boundaries6.3 Waves of large amplitudes. Solitons6.4 Rupture of an embankment-type dam6.5 Gas ow in a combustion engineList of codes (1 - 56)see: www.crcpress.com/downloads(F.Cap, Mathematical Methods in Physics and Engineering withMathematica, ccrcpress and of a parachutistc2 Differentiate and integratec3 The differential equation describing the spread of an epidemic diseasec4 The Wronskian of exp(x), exp(2x)c5 The most general linear differential equation of second orderc6 Inhomogeneous equation of oscillationsc7 An initial value problemc8 Homogeneous boundary value problemc9 Inhomogeneous boundary value problemc10 Pick out values of a numerical solution of a problem with varying boundary valuesc11 Preparing the shooting method for inhomogeneous boundary value problemsc12 Series expansion of a Lie-series solutionc13 Learn a loop for the shooting methodc14 Limit cycle of the Van der Pol equationc16 Phase portrait of the Dufng equationc17 Phase portrait Mathieu equationc18 Phase portrait of the pendulum equationc19 Jacobian matrix of spherical coordinatesc20 Vector analysis: curl in spherical coordinatesc21 Separation setup for the Laplacianc22 Laplace transformationc23 Characteristics of one-dimensional ow (3.3.27), see page 119 of the bookc24 Solution of the heat conduction equation using a similaritytransformationc25 Harmonic polynomials as solution of the Laplacianc26 Biharmonic polynomials solve the static homogeneous plate equationc27 Generalized Bessel and Kummer equationsc28 Whittaker, Gegenbauer and Weber equationc29 Legendre equation with Rodriguez formulac30 Laguerre equation with Rodriguez formulac31 Hermite equation with Rodriguez formulac32 Vector eld - cover of the bookc33 Conformal mappingc34 2D ow around a cylinderc35 Solution of the damped Helmholtz equationc36 Free bending vibrations of a rodc37 Two boundaries for the Laplace equationc38 Inhomogeneous boundary problem for the Laplace equationc39 Nontrivial homogeneous boundary problem of the Laplace eigenvalues of a circular membrane in CARTESIAN COORDINATESc41 Eigenvalue problem for a circular membrane. Calculation in cartesian coordinatesc42 Clamped Cassini membrane in cartesian coordinatesc43 Circular membrane with varying surface mass density. Includes four-fold loopc44 Circular plate with two homogeneous boundary conditionsc45 Laplace equation in polar coordinates with a singularityc46 Laplace equation with two different boundaries, inner inhomogeneous on a square and outer homogeneous on circlec47 Membrane with an homogeneous boundary condition on an inner circle and an inhomogeneous outer boundary conditionc48 Corner in a boundary curvec49 Steepening up of a large-amplitude wavec50 Envelope solitonc51 Laplace equation with outer inhomogeneous boundary values on square and inner homogeneous boundary values on circlec52 Gaussian eliminationc53 The shooting methodc54 Navier solutions of the plate equationc56 Toroidal wave guides =hI! ive this cauchy problem but im not able to solve it.any helP? :Py-(1/x*y)=x^3y(1)=-1 =Na srpskom ne postoji, izuzev opste numerike i uglavnom konacnih razlika kodZ.Petrovica.Milovan Peric pise na engleskom.Mozda braca Hrvati imajunesto?Cuo sam za Muzaferiju i Lileka, mozda su oni nesto pisali na nasemtrojeziku? =Does anyone have a rough knowledge in the performanceof various ways to solve ODE:s on a uniform time grid of about24 points?I am considering the following:* Ordinary ODE-solvers* FEM-methods (e.g. cG(q))* CollocationShould there be more methods added to the list for consideration?Does anyone have a hint of what levels of computational complexityI will encounter when trying the different methods? Are ODE-solversreally that advantageous compared to the other methods???The problem I am working with is stiff.Any thoughts?-- /Stefan =I have been posting this in sci.stat.math but maybe this is the place where it should have gone:Can anyone explain me how to calculate the condence intervals from the jacobian returned by minpacks lmdif?Is the following correct?conf ~ diagonal(inverse(matrixmultiply(transpose(jacobian),jacobian) )) => I have been posting this in sci.stat.math but maybe this is the place > where it should have gone:> Can anyone explain me how to calculate the condence intervals from the > jacobian returned by minpacks lmdif?You can see MINPACK covar.f which is designed to work with the outputof lmdif. http://scicomp.ewha.ac.kr/netlib/minpack/Good luck,Craig-- ------------------------------------------------------------- -------------Craig B. Markwardt, Ph.D. EMAIL: craigmnet@cow.physics.wisc.eduAstrophysics, IDL, Finance, Derivatives | Remove net for better response----------------------------------------------------- --------------------- =>>Can anyone explain me how to calculate the condence intervals from the >>jacobian returned by minpacks lmdif?> You can see MINPACK covar.f which is designed to work with the output> of lmdif.That statistics...Ill get the standard errors from the square root of the diagonal of the covariance matrix. Then I multiply them be e.g. 95% and get the condence intervals, right?Christian = >>Can anyone explain me how to calculate the condence intervals from the >jacobian returned by minpacks lmdif? >>>> You can see MINPACK covar.f which is designed to work with the output >> of terms of statistics... >Ill get the standard errors from the square root of the diagonal of the >covariance matrix. Then I multiply them be e.g. 95% and get the >condence intervals, right? >Christian > no. the diagonal of the inverse of the covariance matrix gives you a weightwhich is needed in the computation of the intervals of condence. the computation also involves the quantiles of the t-distribution and theleast sqaures error. look up an applied statistics text from my annotation: (from Petr Kuzmic)begin{citation}The best reference I can think of is %-------------------@book { Bate88,author = D. M. Bates and D. W. Watts,title = Nonlinear Regression Analysis and Its Applications,publisher = Wiley, address = New York, year = 1988}There are _ve_ different kinds of condence intervals you might beinterested in: (a) The approximate marginal condence interval (whichis the one you are asking about) is equation 1.12 on page 6. (b) Theapproximate condence band is equation 1.13 on page 6. (c) Thelikelihood region is equation 1.14 on page 6. (d) The Bayesianinference region is equation 1.17 on page 7.These four kinds of condence intervals (including type (a) which youasked about) have one very serious disadvantage: they assume that thetting model is sufciently _linear_ in parameters. That is, themodel is assumed to be approximately linear within the condenceinterval for the parameters. [Yes, somewhat circular logic.]A kind of condence interval that is very much preferable to all fourtypes above is (e) likelihood region determined by a systematic _search_of the least-squares surface in the parameter space. It is still basedon a _linear_ statistical theory, but much less so than (a) through(d). For (e) there isnt a formula similar to the one you are askingabout. Instead, there is a search algorithm referred to as theprole-t method. The method is sketched out on page 205, and anactual algorithm is given in pseudo-code on pages 302-303.---- =Can anybody recommend a good textbook onData Structures and Algorithms for Scientists and Engineers? => Can anybody recommend a good textbook on> Data Structures and Algorithms for Scientists and Engineers?MaybeAlgorithms in C, Robert Sedgewick, Addison-Wesley, ReadingI have only looked at parts, but it looks good. I think he has a general book on algorithms not narrowed to C, too.-- Lou Pecora - My views are my own. => Can anybody recommend a good textbook on> Data Structures and Algorithms for Scientists and Engineers?> Maybe> Algorithms in C, Robert Sedgewick, Addison-Wesley, Reading> I have only looked at parts, but it looks good.> I think he has a general book on algorithms not narrowed to C, too.> --> Lou Pecora> - My views are my own.Try Numerical Recipes by Press, et al. in any of your favoritelanguages. (There is a Fortran version, a C version, and--IIRC--a C++ version. Also Pascal, in case anyone is still using that.)-- Julian V. NobleProfessor Emeritus of ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. =Its ages since Ive done any math, can anybody show me how to denean equation from the following data??x|0|1|10--------y|0|2|10 =Heres one way:assume a parabolic of form y=ax^2+bx+cthen1) y=0, x=0 => c=02) y=2,x=1 => 2 = a + b3) y=10,x=10 => 10 = 100a + 10bsolving 2 & 3 for a and b:a= -1/9, b=19/9so egn is y = -x^2/9 +19x/9Jim>> Its ages since Ive done any math, can anybody show me how to dene> an equation from the following data??>> x|0|1|10> --------> y|0|2|10> <3ED72AC6.1E61D610@univie.ac.at> =Yes - I agree - there are many examples where A,Bpos defand AB+BA not pos.sem.defBut is there some criteria to determine when this happens?It appears to have a lot to do with the angles between eigenspacesi.e. if P,Q are the orthogonal matrices that diagonalize A,B, resp. thenit is related to the condition number of PQFor example if PQ=I then pos. def. is preserved, i.e. if A,B commute thepos. def. is preserved.So - conjecture - there is an angle condition dependent on the dimensionthat guarantees pos.def. is preserved.Perhaps what is needed is the measure of nonnormality norm(AB-BA)Univ. of Waterloo |URL http://orion.math.uwaterloo.ca/~hwolkowi>I think that the original denition for positive denite referred to>a --- quadratic form --- and that is why there is confusion with the>denition, i.e. a quadratic form is positive denite if> 0 for all x not equal 0> and so this is clearly true iff A+A* is positive denite.> (I suspect that quadratic forms came before the study of matrices, i.e.> before Sylvestor.) So the denition in Horn-Johnson makes sense.>> I have a related question - or converse question. Suppose you have two> positive denite matrices A,B (real, symmetric for simplicity).> When is K=AB positive denite (i.e. AB+BA positive denite)> or positive semidenite. I think that this relates to the study of so-called> K-pd info>>-->Univ. of Waterloo |URL http://orion.math.uwaterloo.ca/~hwolkowi>> this damned problem puzzled me some years ago at no avail until ...>>> A=diag([1 100]);>> U=[0.8 0.6 ; 0.6 -0.8];>> B=diag([100 1]);>> C=U*A*U*B;>> C>> C =>> 1.0e+03 *>> 3.6640 -0.0475> -4.7520 0.0644>>> D=C+C;>> eig(D)>> ans =>> 1.0e+03 *>> -2.2710> 9.7278> sorry> peter>-- Univ. of Waterloo |URL http://orion.math.uwaterloo.ca/~hwolkowi <3ED72AC6.1E61D610@univie.ac.at> = [ deleted ]It is IMPOSSIBLE to dene an ordering for a non-hermitian matrix. Thereason is basically the same as the reason we cant dene an orderingfor complex numbers.An operator or matrix M is positive denite if (x,Mx) > 0for ANY vector x in the vector space. Then Im (x,Mx) = 0 for all x.(Otherwise you are trying to nd an ordering for compex numbers!)Now, (x,Mx)* = (x,Mx) and since x is arbitrary, this implies that (M^T)^* = M M adjoint = Mi.e. M is hermitian.Now for your case, namely A and B both hermitian and positive, why isnt AB + BA > 0 ?The answer is really simple: write AB + BA = (1/2) [ (A+B)^2 - (A-B)^2 ]that is, as the difference between two positive denite matrices.The problem is this: although A+B > A-B, it is not in generaltrue that (A+B)^2 > (A-B)^2 . It would only be true if A and Bcommute, for then they can be diagonalized simultaneously andthe ordering relation (A+B)^2 > (A-B)^2 would then reduce toa relation on the individual eigenvalues in the diagonalizingbasis, i.e. the relation for real numbers. But of course forx,y real and positive we trivially have xy + yx > 0.-- Julian V. NobleProfessor Emeritus of ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. =Suppose I have a square symmetric similarity matrix where M{ij} is thenumber of joint occurrences of dichotomous attributes i and j in apopulation n. For example, there are 20 attributes and the matrixrepresents the number of cases containing each each pair ofattributes. Is it possible to derive from this the generatingcase-by-attribute matrix? Is there a unique solution? =I have two logarithmic decibel values, one at the start of the soundand one at the end of a sound displayed on an arithmetic graph. Howcould I determine the proper volume level at any point in between thestart and the stop of the sound? I initially thought that a ratiowould to it, but that doesnt take into account the logarithmicproperty of sound. Any help would be appreciated. =I thought that a good pseudo random number generator could generate twouncorrelated sequences of random numbers with two different seeds. Letthe random number generator be rand() and it takes a seed to initializeit. So I thought rand(-1) and rand(-2) should produce two sequences ofuncorrelated random numbers. But my professor told me that it was notguaranteed and the two sequences were possibly correlated or even highlycorrelated. He said only two sequences of random number generated withthe same seed could be uncorrelated. What do you think about this? Idoubt what he said. If a random number is good enough, such as ran1 andran2 in Numerical Recipes, it should produce two uncorrelated sequences ofrandom numbers with two different seeds. David =>>I thought that a good pseudo random number generator could generate two>uncorrelated sequences of random numbers with two different seeds. Let>the random number generator be rand() and it takes a seed to initialize>it. So I thought rand(-1) and rand(-2) should produce two sequences of>uncorrelated random numbers. But my professor told me that it was not>guaranteed and the two sequences were possibly correlated or even highly>correlated. He said only two sequences of random number generated with>the same seed could be uncorrelated. What do you think about this? I>doubt what he said. If a random number is good enough, such as ran1 and>ran2 in Numerical Recipes, it should produce two uncorrelated sequences of>random numbers with two different seeds. >Most generators are of the form x_n = f(x_{n-1}), where f is a carefullydesigned function that garantees that x_n and x_{n-1} are not correlated.The initial value derives from the seed that you provide, sayx_0 = g(s)Now, imagine that you happen to choose s_1 such that x_0 = g(s_1), ands_2 such that x_1 = g(s_2). Do you think that the two sequencesare uncorrelated ?Even with a less extreme case, you dont know of the method to choosetwo seeds so that they do not introduce a correlation.If you want 2 uncorrelated sequences, take the sub-sequences of odd and evenindices in a single sequence that you know is perfectly uncorrelated.Sub-sequences of a perfectly uncorrelated sequence are uncorrelated byconstruction. =Sorry, thats incorrect. Let n be a seed. In a properly-constructedRNG, the sequences S1(n) and S2(n+1) are _NOT_ correlated. Any correlation coefcient thatyou compute shouldnot be higher than one would expect by chance. While its obviouslytrue that two sequenceswill correlate perfectly (in fact be equal), given equal seeds, itreally shouldnt matter whether the two seeds differ by 1 or 1,000,000. Distance is no measureof correlation.Jack>>>> Most generators are of the form x_n = f(x_{n-1}), where f is a carefully>> designed function that garantees that x_n and x_{n-1} are not correlated.>> The initial value derives from the seed that you provide, say>> x_0 = g(s)>>> Now, imagine that you happen to choose s_1 such that x_0 = g(s_1), and>> s_2 such that x_1 = g(s_2). Do you think that the two sequences>> are uncorrelated ?>>It sounds that you are saying the two sequences are correlated to whatever>extend it might be. But I cant quite understand your argument. What is>> Lets take the simple example of a random :> if ( seed provided ) saved_state = seed> saved_state = (saved_state * 843314861 + 453816693) modulo 2^31> random = saved_state / 2^31> If you start a sequence with seed = 1, the next value of state> is 1297131554.> If you start another sequence with seed = 1297131554, that> sequence is the rst one shifted by one position, thus> perfectly correlated to the rst one.> Thats the extreme case.> Now the sequences here are periodic with a period P less than 2^31.> There are a number of possible cycles, depending on the seed.> If you choose a second seed close to the rst one in one cycle,> we just saw the sequences are just shifts one from the other.> If you choose a second seed in another cycle, the two cycles may> run in parallel with high correlation.> Of course, not all generators are as simple as this one, but similar> reasons exist that may create correlation.> Michel => Sorry, thats incorrect. Let n be a seed. In a properly-constructed> RNG, the sequences S1(n) and S2(n+1) are _NOT_ correlated. Any correlation > coefcient that you compute should not be higher than one would expect by > chance. This is surely just a question of denition. I have certainly useda random number generator which was good in other ways, but had unexpectedcoincidences for different seed values early on in the sequence. Looking at myarchives, it turns out it was a Tausworthe generator. =For many random number generators, your professor would be wrong. Different seeds do not produce different sequences, only differentstarting points in a single sequence.Most RNGs are of the form n(i+1) = (A * n(i) + B) mod mThe Numerical Recipes generate ran0 drops the B. Such generatorsdepends _ONLY_ on the current value of n to compute the next value. IfA and B are chosen properly, they will have a long run sequence, meaningthat they will generate a large number of different integers beforereturning to the original value again. For example, in ran0 the runlength is about 2 billion. If you start with a known seed, say 1, theRNG will generate 2 billion integers, one time each, before returningback to 1 again. After that, of course, the sequence repeats. Think of it as a huge bike wheel, in which each spoke is labelled with adifferent address. As you advance the wheel, all you do is switch tothe next spoke in the wheel, with a different number written on it.The starting seed simply decides where the wheel is when the sequence isstarted. The sequence doesnt change, only the place where one startsout on it.ran1 is a bit different. It uses a shufe table that shufes, alsoin random fashion, the order in which values are drawn from thesequence. For an RNG like that, your professor would be _PARTIALLY_right, in the sense that each seed creates a unique initial set of tableentries, and therefore a unique sequence.I think what your prof is saying is that, with such an RNG, one cannotprove, via mathematics, that rand(1) and rand(2) [forget the -signs ...thats an artifact of the implementation] are uncorrelated. Likewise,its impossible to prove, except by direct test, that they _ARE_correlated. In the end, all one can do is run a whole bunch of tests and compute awhole bunch of correlation coefcients. According to Press, et al., noone has ever demonstrated a case where ran1 didnt pass such tests. Your prof would say, Ok, but that doesnt prove that it wont fail thenext test. He would be quite correct in saying so, but the opinion isonly of academic interest.Jack> I thought that a good pseudo random number generator could generate two> uncorrelated sequences of random numbers with two different seeds. Let> the random number generator be rand() and it takes a seed to initialize> it. So I thought rand(-1) and rand(-2) should produce two sequences of> uncorrelated random numbers. But my professor told me that it was not> guaranteed and the two sequences were possibly correlated or even highly> correlated. He said only two sequences of random number generated with> the same seed could be uncorrelated. What do you think about this? I> doubt what he said. If a random number is good enough, such as ran1 and> ran2 in Numerical Recipes, it should produce two uncorrelated sequences of> random numbers with two different seeds.> David =>> Most generators are of the form x_n = f(x_{n-1}), where f is a carefully>> designed function that garantees that x_n and x_{n-1} are not correlated.>> The initial value derives from the seed that you provide, say>> x_0 = g(s)> Now, imagine that you happen to choose s_1 such that x_0 = g(s_1), and>> s_2 such that x_1 = g(s_2). Do you think that the two sequences>> are uncorrelated ?>>It sounds that you are saying the two sequences are correlated to whatever>extend it might be. But I cant quite understand your argument. What is>Lets take the simple example of a random :if ( seed provided ) saved_state = seed saved_state = (saved_state * 843314861 + 453816693) modulo 2^31random = saved_state / 2^31If you start a sequence with seed = 1, the next value of stateis 1297131554.If you start another sequence with seed = 1297131554, thatsequence is the rst one shifted by one position, thusperfectly correlated to the rst one.Thats the extreme case.Now the sequences here are periodic with a period P less than 2^31.There are a number of possible cycles, depending on the seed.If you choose a second seed close to the rst one in one cycle,we just saw the sequences are just shifts one from the other.If you choose a second seed in another cycle, the two cycles mayrun in parallel with high correlation.Of course, not all generators are as simple as this one, but similarreasons exist that may create correlation.Michel =So what are we saying here? If one sets out to nd a random numbergenerator sufciently bad,one can nd it?Jack> Sorry, thats incorrect. Let n be a seed. In a properly-constructed> RNG, the sequences S1(n) and S2(n+1) are _NOT_ correlated. Any correlation> coefcient that you compute should not be higher than one would expect by> chance.> This is surely just a question of denition. I have certainly used> a random number generator which was good in other ways, but had unexpected> coincidences for different seed values early on in the sequence. Looking at my> archives, it turns out it was a Tausworthe generator. = I need some routines that implement directed rounding using onlyround-to-nearest oating point operations for +, -, * and /.Actually, I only need round to -Innity and to +Innity. Also, whengiven a oating point number x, that I know it is the result of anoperation with 1 ulp precision, how can I get the smallest oatingpoint interval that contains the result ? I need to do these in Java,but some pseudocode would be sufcient.Daniel Aioanei =can anyone offer some advice?im writing an application the numerically solves the diffusionequation in 2D. the spatial domain is rectangular, with equallyspaced x and y gridpoints. im using a alternating-direction, fullyimplict nite differencing scheme (backward time, center space). theboundries are neumann on the top an bottom (du/dn = 0) andux-conservative (cyclic) on the left and right.i would like to be able to assign different diffusion coefcients tocertain portions of the 2D space, in order to simulate diffusionthrough multiple materials in the same grid space. however, i am notable to apply discontinuous diffusion coefcients - i get artifactsat the boundry between two different coefcients using this method.can anyone offer some advice on how to solve this problem?-john biafore =>> can anyone offer some advice?>> im writing an application the numerically solves the diffusion> equation in 2D. the spatial domain is rectangular, with equally> spaced x and y gridpoints. im using a alternating-direction, fully> implict nite differencing scheme (backward time, center space).the> boundries are neumann on the top an bottom (du/dn = 0) and> ux-conservative (cyclic) on the left and right.>> i would like to be able to assign different diffusion coefcientsto> certain portions of the 2D space, in order to simulate diffusion> through multiple materials in the same grid space. however, i amnot> able to apply discontinuous diffusion coefcients - i get artifacts> at the boundry between two different coefcients using this method.>Irrespective of the grid and the solution scheme, you havent givensufcient details of the prevailing physics. Can you elaborate?Ciao,Gerry T. =My requirement is as follows.I have datapoints (x1, y1) (x2,y2) (x3,y3) and so on till ( Xn , Yn).Now i have estimate/predicate the next value ie ( Xn+1 , Yn+1) usingleast square method with order 1 (ie using straight line). Could youplease let me know the formula which i should use to do the same...Sankar => I have datapoints (x1, y1) (x2,y2) (x3,y3) and so on till ( Xn , Yn).> Now i have estimate/predicate the next value ie ( Xn+1 , Yn+1) using> least square method with order 1 (ie using straight line). Could you> please let me know the formula which i should use to do the same...> You have to nd the function f(x)=ax+b by minimising the followingexpression sum_{i=1}^n (f(x_i) - y_i)^2Then you take y_{n+1}=f(x_n+1)-- Wit Jakuczun =Im trying to plot a histogram on some data that I receive at runtime. Im trying to measure distances between two points, and these distances can be from the range 1 to 5000000. I want to use a logarithmic scale with a range -2^29 to +2^29. How can I determine at any point in time which strata (i.e which cohort or portion of scale) the value belongs to (without using the square or square-root functions)? Im trying to implement this on a computer.I could have mask-bits of all possible scales, such as 2^2, 2^3, 2^4 upto 2^29 and AND each of these with the value to gure this out but is there a cleaner way of doing this?--Andre => Im trying to plot a histogram on some data that I receive at runtime. > Im trying to measure distances between two points, and these distances > can be from the range 1 to 5000000. I want to use a logarithmic scale > with a range -2^29 to +2^29. How can I determine at any point in time > which strata (i.e which cohort or portion of scale) the value belongs to > (without using the square or square-root functions)? Im trying to > implement this on a computer.> I could have mask-bits of all possible scales, such as 2^2, 2^3, 2^4 > upto 2^29 and AND each of these with the value to gure this out but is > there a cleaner way of doing this?A distance cant be less than zero...Maybe you meant to write: 1/(2^29) to 2^29?In some languages, like C, you can use bit-shift operations to getapproximations of logs in base 2 of integers...Anyway, for integer inputs, maybe this can be of some help:David Bernier-----------------------------------------------# include int main(void){ long input, inputcopy; int sign, power; printf(nplease enter your number:n); scanf(%ld, &input); inputcopy = input; if(input>=0) { sign = 1; } else { sign = -1; inputcopy = -input; } power = 0; while((inputcopy>>power)>0) /*** >> bitwise-shift ***/ { power++; }printf(ninput = %ldtsign = %dtpower = %dnn, input, sign, power); return 0;}---------------------------------------------------------- - =Actually in my implementation I get negative distances the code excerpt :)--Andre>> Im trying to plot a histogram on some data that I receive at runtime. >> Im trying to measure distances between two points, and these >> distances can be from the range 1 to 5000000. I want to use a >> logarithmic scale with a range -2^29 to +2^29. How can I determine at >> any point in time which strata (i.e which cohort or portion of scale) >> the value belongs to (without using the square or square-root >> functions)? Im trying to implement this on a computer.>> I could have mask-bits of all possible scales, such as 2^2, 2^3, 2^4 >> upto 2^29 and AND each of these with the value to gure this out but >> is there a cleaner way of doing this?> A distance cant be less than zero...> Maybe you meant to write: 1/(2^29) to 2^29?> In some languages, like C, you can use bit-shift operations to get> approximations of logs in base 2 of integers...> Anyway, for integer inputs, maybe this can be of some help:> David Bernier> -----------------------------------------------> #include > int main(void)> {> long input, inputcopy;> int sign, power;> printf(nplease enter your number:n);> scanf(%ld, &input);> inputcopy = input;> if(input>=0)> {> sign = 1;> }> else> {> sign = -1;> inputcopy = -input;> }> power = 0;> while((inputcopy>>power)>0) /*** >> bitwise-shift ***/> {> power++;> }> printf(ninput = %ldtsign = %dtpower = %dnn, input, sign, power); > return 0;> }> ----------------------------------------------------------- => A colleague asked me for algorithms performing efcient evaluation of> spherical harmonic functions.> Can someone suggest a reference? I havent worked with harmonics for some> years.> As always,> Ed Wright> [...]The best place to start is probably NCARs SPHEREPACK 3.0: A Model Development Facilityhttp://www.scd.ucar.edu/css/software/spherepack/which has a lot of tools for spherical harmonics.Hugo Pfoertnerhttp://www.pfoertner.org/ = >A colleague asked me for algorithms performing efcient evaluation of >spherical harmonic functions. >Can someone suggest a reference? I havent worked with harmonics for some >years. http://www.netlib.org/toms/629 from Atkinson contains some code for computing these.hope this helpspeter =If you have to perform sum_{k,n} Y^n_k for non-uniform distributed nodes on thesphere you may read Fast spherical Fourier algorithmshttp://www.math.mu-luebeck.de/workers/potts/ publikationen/stefan =Suppose I have a series of numbers in a one-dimensional array.I want to nd the median of 5 or 7 numbers centered about thepoint of interest. I must effectively sort these numbers to nd themedian. Now I want to do this for the whole array. Each time,I roll one number off the group and introduce a new number tothe group. Does this allow an algorithm with less operations tond the new median, or must I sort the group all over? => Suppose I have a series of numbers in a one-dimensional array.> I want to nd the median of 5 or 7 numbers centered about the> point of interest. I must effectively sort these numbers to nd the> median. Now I want to do this for the whole array. Each time,> I roll one number off the group and introduce a new number to> the group. Does this allow an algorithm with less operations to> nd the new median, or must I sort the group all over?You can use a (balanced) tree which will allow you to add and deletein O(log(n)) but if you are talking about 5 to 7 numbers I think youare better off with a brand new linear insertion sort each time.(If I remember correctly, Knuth showed that for a MIX implementationthe cutoff between linear insertion sort and better sorting methodswas at 11 elements. Current processors are much different but Iwould expect that for 5 to 7 elements linear insertion is still best.There are also specic algorithms for order statistics which might bemore efcient. Again, for 5 to 7 elements I dont think it will makemuch difference.)-- Use of tools distinguishes Man from Beast. And UNIX users from WINDOZE lusers. =>Suppose I have a series of numbers in a one-dimensional array.>I want to nd the median of 5 or 7 numbers centered about the>point of interest. I must effectively sort these numbers to nd the>median. Now I want to do this for the whole array. Each time,>I roll one number off the group and introduce a new number to>the group. Does this allow an algorithm with less operations to>nd the new median, or must I sort the group all over?For any order statistic, there is an algorithm to nd thatorder statistic in time O(N). This proceeds by computingmedians, but it is not the case that it can do so byupdating a current median.The algorithm is not complicated, although it may look soat rst glance.If N is sufciently small, updating the full order may bemore efcient.-- This address is for information only. I do not claim that these viewsare those of the Statistics Department or of Purdue University.Herman Rubin, Deptartment of Statistics, Purdue University matrix A (nearly band-diagonal, maybe Hermitian),> and I need to compute trace(1/A) as efciently as possible for large> N.If you can afford factoring your matrix, you can compute the sparse inverse (see the book by Duff, Erisman, Reid) and take its trace.If your sparsity pattern is unfavorable for direct methods,you can probably adapt the methods used in math.asu.edu/~hongbin/biformslide_L.pdf for computing u^*f(A)v approximately.Arnold Neumaier =Im facing the following problem:Given a set of point in R3 Im looking for the ellipsoid of minimalsurface and arbitrary semiaxes and orientation that circumscribes thegiven points. Ive already posted this question in the Matlabnews-group* and got a quite useful reply of how to nd the minimizerof the volume. Finding the minimizer of the surface seems to be moreinvolved ;-)Any threadm=bemfgv%24p7u%241%40ginger.mathworks.com&rnum=1&prev=/ groups%3Fq%3Djohan%2Bfranzen%2Bellipsoid%26hl%3Dde%26lr%3D% 26ie%3DUTF-8%26oe%3DUTF-8%26selm%3Dbemfgv%2524p7u%25241% 2540ginger.mathworks.com%26rnum%3D1 = >Im facing the following problem: >Given a set of point in R3 Im looking for the ellipsoid of minimal >surface and arbitrary semiaxes and orientation that circumscribes the >given points. Ive already posted this question in the Matlab >news-group* and got a quite useful reply of how to nd the minimizer >of the volume. Finding the minimizer of the surface seems to be more >involved ;-) >you write the ellipsoid as (x-x0)*L*L*(x-x0) =1 with x0 and the triangular matrix L as unknowns, with L(i,i)>=eps>0 eps specied (very small, this as a safety in order to avoid the degenerate case. with a proper point cloud these constraints will never be active, hence choice of eps is not critical. then you constraints are (x(i)-x0)*L*L*(x(i)-x0) <=1 for the given points x(i) in R^3. the surface is now given by the eigensystem decomposition of L*L -> eig in MATLAB and the computation of the surface using the ordinary surface integral (an elliptic integral) and depends solely on the three eigenvalues of L*L.you must do this numerically. routines for these integrals are available somewhere (at least in f77 nad c in netlib.)this all gives a nonlinear optimization problem with nonlinear inequality constraints, a job for fmincon (in MATLAB optimization toolbox).hence you have 9 unknowns, a nonlinear objective function anda set of nonlinear inequality constraints, but everything smooth.hope this helpspeter =>Im facing the following problem:>>Given a set of point in R3 Im looking for the ellipsoid of minimal>surface and arbitrary semiaxes and orientation that circumscribes the>given points. Ive already posted this question in the Matlab>news-group* and got a quite useful reply of how to nd the minimizer>of the volume. Finding the minimizer of the surface seems to be more>involved ;-)>> you write the ellipsoid as > (x-x0)*L*L*(x-x0) =1 > with x0 and the triangular matrix L as unknowns, with L(i,i)>=eps>0 > eps specied (very small, this as a safety in order to avoid the degenerate > case. with a proper point cloud these constraints will never be active, hence > choice of eps is not critical. > then you constraints are > (x(i)-x0)*L*L*(x(i)-x0) <=1 > for the given points x(i) in R^3. the surface is now given by the > eigensystem decomposition of L*L -> eig in MATLAB > and the computation of the surface using the ordinary surface integral > (an elliptic integral) and depends solely on the three eigenvalues of L*L.> you must do this numerically. routines for these integrals are available > somewhere (at least in f77 nad c in netlib.)> this all gives a nonlinear optimization problem with nonlinear inequality > constraints, a job for fmincon (in MATLAB optimization toolbox).> hence you have 9 unknowns, a nonlinear objective function and> a set of nonlinear inequality constraints, but everything smooth.> hope this helps> peterIndeed, ... thank you for your clarifying words, Peter. Your help isgreatly appreciated.johan =I still have troubles with polynomials.I have a polynomial in the forma_0 + a_1 x + a_2 x^2 + ... + a_n x^n =p(x) (1)(everything real here)How to plot this polynomial in a reliable way, nearby a root ?That is, how to plot, for example p(x)=(1-x)^20nearby x=1, using the expression (1) ?I have the book Elementary numerical analysis by Conte and de Boor on my desk,so I refer to page 73, where they show my problem, :-))))even if they dont seem to provide a solution to it. :-((((Luciano--Luciano RibichiniMax Planck Insitute for Gravitational PhysicsCallinstr. 38D 30167 Hannover Germany = >I still have troubles with polynomials. >I have a polynomial in the form >a_0 + a_1 x + a_2 x^2 + ... + a_n x^n =p(x) (1) >(everything real here) >How to plot this polynomial in a reliable way, nearby a root ? >That is, how to plot, for example > p(x)=(1-x)^20 >nearby x=1, using the expression (1) ? >I have the book Elementary numerical analysis by Conte and de Boor on my desk, >so I refer to page 73, where they show my problem, :-)))) >even if they dont seem to provide a solution to it. :-(((( > think about the size of the coefcients in (1) if you expand (1-x)^20 ! the largest one will be 184756. now in order to get an adequate approximationsay for x=0.9 (i.e. 1.0e-20), you need about 28 digits precision for 1% (graphics) precision. this is your problem and there is no way around other than:dont use the representation (1) for a polynomial. you will always need to represent it in a transformed scale (and possibly in some other base of the polynomials) hope this helpspeter =>I still have troubles with polynomials.>>I have a polynomial in the form>>a_0 + a_1 x + a_2 x^2 + ... + a_n x^n =p(x) (1)>>(everything real here)>>How to plot this polynomial in a reliable way, nearby a root ?>>That is, how to plot, for example>> p(x)=(1-x)^20>>nearby x=1, using the expression (1) ?> Let me quote Pete Stewart(Afternotes on Numerical Analysis, ISBN 0-89871-362-5),Lecture 8: What has killed us is not the certicate.In this spirit, my condolences. As Peter Spellucci explained below, the software is tortured tocombine huge numbers (up to binomial(20,10)) to obtain a tinynumber. I did mention earlier that the power form (1) is often theclumsiest form of a polynomial; if it was originally in anotherform then even the conversion into power form is laden witherrors, not to mention subsequent use in evaluation.(More from me below.)>I have the book Elementary numerical analysis by Conte and de Boor on my desk,>>so I refer to page 73, where they show my problem, :-))))>>even if they dont seem to provide a solution to it. :-((((> think about the size of the coefcients in (1) if you> expand (1-x)^20 !> the largest one will be 184756. now in order to get an> adequate approximation say for x=0.9 (i.e. 1.0e-20), you> need about 28 digits precision for 1% (graphics)> precision. this is your problem and there is no way around> other than> :dont use the representation (1) for a polynomial. you will> always need to represent it in a transformed scale (and> possibly in some other base of the polynomials)> hope this helps> peterThe original inquiry did not supply enough information:Are the coefcients exact (for example, guaranteed tobe integers of manageable size)?Are the roots guaranteed to be multiple, rather than justclustered? In the yes-yes case, there is a way: the greatest commondivisor d(x) of p(x) and of p(x) will be non-constant, andq(x)=p(x)/d(x) (symbolic division all along) will have the sameroots as p(x), each of them exactly once. Lowering the multiplicity of factors, and then restoring thepolynomial p(x) as the product of powers of irreducible factorswill lead to much more stable evaluation. Factors involving the roots of (temporary) interest can befurther improved by Taylors expansion centered near thesingled-out root. Attention: Power form (1) is Taylors expansion centered at 0, and there is no promise of good performance far away from 0. Thats the lost information on the way to the funeral.The usual scenario: The coefcients of p(x) in form (1) arerounded real numbers; in this case, multiple roots are a myth,and clusters of simple roots are a reality with probability 1.(Find out about Weierstasss Preparation Theorem which tells youhow the perturbations of multiple roots travel in areworks-like fashion from their centre.)In this case, always suspect the form (1) as coming from anotherform, with possibly many roundings along the way.Concerning (1-x)^20 versus its binomial expansion: it brings theold saying to mind if it aint broke, dont x it. =I have implemented a Poisson solver for a simple 3d rectilinear gridwith equal sized voxels. I apply Neumann boundary conditions (as faras I remember :), writing out a linear system of equations, and solvethis with the Conjugate Gradient (CG) method (with optional incompleteCholesky preconditioning).But recently Ive heard of so called Fast Poisson Solvers. As solvingthe Poisson equation is at the moment the most time-consuming part ofmy simulator, my interest was hooked. Still, I found very littleinformation on these methods. I have three questions regarding thesubject, which Ill write in a simple list,1.) What (good) online papers are there out of the subject?2.) Ive heard that the Fast Poisson Solvers include a method based onthe FFT. I also know that its possible to solve the Poisson equationefciently using FFT but only if there are no boundary conditions(ie. the uid wraps around, yes, Im doing uid simulation). CanFast FFT Poisson Solvers overcome this problem of having no boundaryconditions?3.) I also heard that its possible to accelerate CG convergence usinga Fast Poisson Solver as the preconditioner. If a FFT based FastSolver is used in this case, shouldnt it have to be able to handleboundary conditions so that the preconditioner would be meaningfulin any way? - Mikko Kauppila =I am no expert, but I think Multi-Grid should also be on your short list.> I have implemented a Poisson solver for a simple 3d rectilinear grid> with equal sized voxels. I apply Neumann boundary conditions (as far> as I remember :), writing out a linear system of equations, and solve> this with the Conjugate Gradient (CG) method (with optional incomplete> Cholesky preconditioning).> But recently Ive heard of so called Fast Poisson Solvers. As solving> the Poisson equation is at the moment the most time-consuming part of> my simulator, my interest was hooked. Still, I found very little> information on these methods. I have three questions regarding the> subject, which Ill write in a simple list,> 1.) What (good) online papers are there out of the subject?> 2.) Ive heard that the Fast Poisson Solvers include a method based on> the FFT. I also know that its possible to solve the Poisson equation> efciently using FFT but only if there are no boundary conditions> (ie. the uid wraps around, yes, Im doing uid simulation). Can> Fast FFT Poisson Solvers overcome this problem of having no boundary> conditions?> 3.) I also heard that its possible to accelerate CG convergence using> a Fast Poisson Solver as the preconditioner. If a FFT based Fast> Solver is used in this case, shouldnt it have to be able to handle> boundary conditions so that the preconditioner would be meaningful> in any way?> - Mikko Kauppila-- Use of tools distinguishes Man from Beast. And UNIX users from WINDOZE lusers. =>> I have implemented a Poisson solver for a simple 3d rectilinear grid> with equal sized voxels. I apply Neumann boundary conditions (as far> as I remember :), writing out a linear system of equations, and solve> this with the Conjugate Gradient (CG) method (with optional incomplete> Cholesky preconditioning).>> But recently Ive heard of so called Fast Poisson Solvers. As solving> the Poisson equation is at the moment the most time-consuming part of> my simulator, my interest was hooked. Still, I found very little> information on these methods. I have three questions regarding the> subject, which Ill write in a simple list,>> 1.) What (good) online papers are there out of the subject?>> 2.) Ive heard that the Fast Poisson Solvers include a method based on> the FFT. I also know that its possible to solve the Poisson equation> efciently using FFT but only if there are no boundary conditions> (ie. the uid wraps around, yes, Im doing uid simulation). Can> Fast FFT Poisson Solvers overcome this problem of having no boundary> conditions?Try searching multigrid methodsByeAntonio => 1.) What (good) online papers are there out of the subject?First of all, ignore the people who advocate multigrid. There are FastSolver codes out there that will give you, straight out of the box notuning required, the solution in optimal time. If thats possible withMG I havent heard of it.Anyway, there is a package called Fishpack. Search for that, and forpapers by its authors, Swarztrauber and Sweet. Heres an overview paper:author = {P.N. Swarztrauber},title = {The methods of cyclic reduction, {F}ourier analysis and the {FACR} algorithm for the discrete solution of {P}oissons equation on a {1977},pages = {490--501},keywords = {FFT, fast fourier, cyclic reduction, Poisson equation}}Regarding your questions on boundary conditions: sorry, Im no expert onfast solvers.V.-- = have a look at:E.F.F. Botta, K. Dekker, Y. Notay, A. van der Ploeg, C. Vuik,F.W. Wubs, P.M. de Zeeuw,How fast the Laplace equation was solved in 1995,J. Applied Numerical Mathematics 24 (1997) 439--455.Paul de Zeeuwhttp://www.cwi.nl/~pauldz/ =I want to model the heat ow from a hot cylinder into a cold one, where the hot one ist standing on the cold one:XXXXXXXXXX XX HOT XX XX XXXXXXXXXXCCCCCCCCCC CC COLD CC CCCCCCCCCCAs the cylinders are made by different materials and due to other reasons, there is a heat resistance between those two that I want model just like convection:dq/dt=h*(T_hot-T_cold)(T_hot and T_cold are the (r-dependant) temperatures at the cold-top/hot-bottom-surface.I want to avoid using coupling variables if possible as this would slow down my simulation a lot. So I wondered if there is a way to dene some kind of boundary condition for the contact surface of the two bodies or maybe an even simpler solution.Any help is highly appreciated,Nils =Im new to Discrete Fourier Transforms and I have a problem.I have constructed a signal in mathematica from a few different sine wavesand am attempting to use a band-pass lter to select one of the components.The ltering seems to work, but the problem is when I do the Inverse FFT toget back to the time domain, the signal is only half the strength Iexpected.The inverse FFT contained some complex terms, which I had not expected(since the original signal was just real data), so I just plotted the realpart.Now after the Nyquist frequency the FFT repeats itself and my lter (beinga bandpass lter) has totally eliminated this repeated part of thespectrum. Perhaps that is why I only end up with half the strength when I dothe Inverse FFT.Any ideas?P.S. I also get the same effect using Fourier Analysis in Microsoft Excel. =the cause of your troubles is that for a real sequence the FFT returns a symmetric spectrum (with respect to the Nyquist frequency). Your ltering must not destroy this symetry! Zdenen Hurak> Im new to Discrete Fourier Transforms and I have a problem.> I have constructed a signal in mathematica from a few different sine waves> and am attempting to use a band-pass lter to select one of the> components. The ltering seems to work, but the problem is when I do the> Inverse FFT to get back to the time domain, the signal is only half the> strength I expected.> The inverse FFT contained some complex terms, which I had not expected> (since the original signal was just real data), so I just plotted the real> part.> Now after the Nyquist frequency the FFT repeats itself and my lter> (being a bandpass lter) has totally eliminated this repeated part of the> spectrum. Perhaps that is why I only end up with half the strength when I> do the Inverse FFT.> Any ideas?> P.S. I also get the same effect using Fourier Analysis in Microsoft Excel. =The book Im using (Neural, Novel & Hybrid Algorithms for Time SeriesPrediction) describes a simple gaussian lter in the frequency domain.i.e. H(f) = exp-((f-f0)/s)^2Where f0 is the centre frequency and s is the width.I applied this right across the frequency domain, and of course, once itwent past the Nyquist frequency(for this particular example) and into themirrored part, the lter frequency f, from the sample numbern.I was just using fn = n/(N*Delta) where N is the total number of inputpoints and Delta is the sampling interval.So as I went from n=0 to n = N-1 the frequency I was plugging in to thelter just got higher and higher as I wentinto the mirrored part of the spectrum.Steven Martin>> the cause of your troubles is that for a real sequence the FFT returns a> symmetric spectrum (with respect to the Nyquist frequency). Your ltering> must not destroy this symetry!>> Zdenen Hurak>> =just apply your lter to the lower half of the frequency response (rst N/2 samples of frequency response). Then build a new ltered frequency response by mirroring the ltered part. Then use the Inverse FFT to get a time-domain signal.You must understand that the higher half of the frequency-domain data (called frequency response in engineering) is somewhat redundant here. The nonredundant part of your frequency-domain data is only as high as the Nyquist frequency. Note however, that if you needed to lter complex time-domain signals, this claim about redundancy and symmetry would no longer be valid.Zdenek> The book Im using (Neural, Novel & Hybrid Algorithms for Time Series> Prediction) describes a simple gaussian lter in the frequency domain.> i.e. H(f) = exp-((f-f0)/s)^2> Where f0 is the centre frequency and s is the width.> I applied this right across the frequency domain, and of course, once it> went past the Nyquist frequency(for this particular example) and into the> mistake is in calculating the frequency f, from the sample> number n.> I was just using fn = n/(N*Delta) where N is the total number of input> points and Delta is the sampling interval.> So as I went from n=0 to n = N-1 the frequency I was plugging in to the> lter just got higher and higher as I went> into the mirrored part of the spectrum.> Steven Martin>> the cause of your troubles is that for a real sequence the FFT returns a>> symmetric spectrum (with respect to the Nyquist frequency). Your>> ltering must not destroy this symetry!>> Zdenen Hurak>> =I had made a fft c program and it is working perfectly for discretevalues. I have also checked it on MATLAB and quite sure that theprogram is correct.However I have two doubts.I inputted the samples of a sine wave of 1000 Hz with a samplingfrequency 4000Hz. For a simple 16 point FFT I got two peaks one atx[4] and other at x[12]. Now shouldnt a simple sine wave show justone peak at 1000Hz.Is it that the values I get are centered along the imaginary axis. Iam not able to understand where the frequency samples after doing fftactually lie in the spectrum.Please help me sorting this fundamental problem. I hope many of youknow the answer.Aarul Jain for a real sequence is SYMMETRIC! That is why you obtain two peaks for a simple sine signal. Just note that with 16 points, the rst peak is at the (0+4th)position while the second peak is at (16-4)th position. Also note that spectrum of a sampled signal is PERIODIC with period equal to the sampling freqency, i.e., 4000 Hz in your case.To summarize, the frequency interval is 0 - 4000 Hz and moreover is symmetric in a sense the the second half of the spectrum is just a mirror image of the rst half.Zdenek Hurak> I had made a fft c program and it is working perfectly for discrete> values. I have also checked it on MATLAB and quite sure that the> program is correct.> However I have two doubts.> I inputted the samples of a sine wave of 1000 Hz with a sampling> frequency 4000Hz. For a simple 16 point FFT I got two peaks one at> x[4] and other at x[12]. Now shouldnt a simple sine wave show just> one peak at 1000Hz. Is it that the \ values I get are centered along the imaginary axis. I> am not able to understand where the frequency samples after doing fft> actually lie in the spectrum.> Please help me sorting this fundamental problem. I hope many of you> know the answer.> Aarul Jain =Wishing understand that the FFT for a real sequence is SYMMETRIC! > That is why you obtain two peaks for a simple sine signal. Just note that > with 16 points, the rst peak is at the (0+4th)position while the second > peak is at (16-4)th position. > Also note that spectrum of a sampled signal is PERIODIC with period equal to > the sampling freqency, i.e., 4000 Hz in your case.> To summarize, the frequency interval is 0 - 4000 Hz and moreover is > symmetric in a sense the the second half of the spectrum is just a mirror > image of the rst half.> Zdenek Hurak> I had made a fft c program and it is working perfectly for discrete> values. I have also checked it on MATLAB and quite sure that the> program is correct.> However I have two doubts.> I inputted the samples of a sine wave of 1000 Hz with a sampling> frequency 4000Hz. For a simple 16 point FFT I got two peaks one at> x[4] and other at x[12]. Now shouldnt a simple sine wave show just> one peak at 1000Hz.> Is it that the values I get are centered along the imaginary axis. I> am not able to understand where the frequency samples after doing fft> actually lie in the spectrum. Please help me \ sorting this fundamental problem. I hope many of you> know the answer.> Aarul Jain 05:56:09 -0700, alexandru.lupas@ulbsibiu.ro>>>To nd a stictly increasing sequence (x_n)_{n>=0} of real numbers >such that>> x_0=2, \ x_1 = r , r being rational > 2 ,and>> (k+1)*x_{k+2}=(k+2)*x_{k} for all k=0,1,... .>>The question semms to be very simple, but >> I havent found any such sequence. Please give an example . >>>>> X(2k+1)/X(2k) = [(2k+1)/(2k) X(2k-1)] / [ (2k)/(2k-1) X(2k-2)]>> = (1 - 1/ (4k^2)) * X(2k-1) / X(2k-2)>>> = Product {t = 1 to k: [1-1/(4t^2)]} r/2>>> The limit as k --> inf of >> Product {t = 1 to k: [1-1/(4t^2)]} = 2/pi>>> Hence for any r < pi, eventually there is a k s.t. X(2k+1) < X(2k). >>> Similarly,>> X(2k)/X(2k-1) = [(2k)/(2k-1) X(2k-2)] / [(2k-1)/(2k-2) X(2k-3)]>> = 4k(k-1)/(2k-1)^2 X(2k-2)/X(2k-3)>> = Product {t = 2 to k: [4t(t-1)/(2t-1)^2]} X(2)/r>> = Product {t = 2 to k: [4t(t-1)/(2t-1)^2]} 4/r>>> Product {t = 2 to k: [4t(t-1)/(2t-1)^2]} --> 4/pi as k --> infSorry for a typo, the above limit is pi/4. The rest holds.>> so if 4/r < 4/pi, we can eventually nd a k s.t. X(2k) < X(2k-1).>>> Hence we need a rational r which is neither lesser nor greater than>> pi...>>Indeed the parameter r is very important. Its unique determined,>namely as you have observed r=pi. Therefore such sequence does not>exists.Alex.> non-equidistant grids, using nite-differencescan cause problems. I cant really imagine why this should be thecase. Is there an easy explanation why this can be a problem or have Ibeen told something wrong?Bye, Christof =>>Ive been told that on non-equidistant grids, using nite-differences>can cause problems. I cant really imagine why this should be the>case. Is there an easy explanation why this can be a problem or have I>been told something wrong?>>Bye,> Christofdepends on your application. if the grid is very irregular, you will getlarger interpolation errors in replacing the derivative (or higher derivatives)by those of the interpolating polynomial. if you use these formulae for discretizing differential equations, you will havemore trouble in order to achieve the normal consistency order, need a largerdiscretization star and so on. also you may loose the M-matrix property in theresulting system and the related stability property. hope this helpspeter => Ive been told that, on non-equidistant grids,> using nite-differences can cause problems.Show us the nite difference equationthat you are using on non-equidistant grids.> I cant really imagine why this should be the case.And we cant *guess* what problem you might be talking about.> Is there an easy explanation why this can be a problem.> Or have I been told something wrong?It is more likely that you simply missedor failed to report some important detailsabout what you have been told. =Im working on a project that requires sketch a surface knowing several(rather many but not enough) points on that surface. I need to estimateother points...Some one told me that Fourier analysis can handle this kind of problem.Could anyone advice what books or online materials are a good ones to readabout Fourier analysis? => Im working on a project that requires sketch a surface knowing several> (rather many but not enough) points on that surface. I need to estimate> other points...>> Some one told me that Fourier analysis can handle this kind of problem.>> Could anyone advice what books or online materials are a good ones to read> about Fourier analysis?>Radial basis functions (RBF) spring to mind. Here one constructs the surfaceusing the scattered shifts of some basis function. If you know the surfaceat the points x_1,...,x_M then one could use:sum_{j=1}^M a_j phi(x-x_j)as an approximation to the surface given an appropriate basis function phi.The coefcients a_j being determined by the ensuing interpolationequations:sum_{j=1}^M a_j phi(x_i-x_j) = f(x_i), i=1,...,M.There are several popular choices for the basis function. One could usephi(x)=exp(-|x|^2/2) for example in any space dimension. In 2-dimensionsphi(x)=|x|^2ln|x| is popular and the resulting interpolant is often calledthe thin plate spline (TPS). Sometimes (for example with the TPS) apolynomial term and additional constraints are added to the above toguarantee the solvability of the system.RBF Interpolants are a multivariate analogy of the well-loved naturalsplines from 1-d.R =In the case of phi(x) = |x|^2ln|x| isnt there a problem when x = x_jwhen generating the a_j coefcients? I think you end up with ln(0),which is not very useful.How are you supposed to work around this?Russell.> Im working on a project that requires sketch a surface knowing several> (rather many but not enough) points on that surface. I need to estimate> other points...>> Some one told me that Fourier analysis can handle this kind of problem.>> Could anyone advice what books or online materials are a good ones to read> about Fourier analysis?>>>> Radial basis functions (RBF) spring to mind. Here one constructs the surface> using the scattered shifts of some basis function. If you know the surface> at the points x_1,...,x_M then one could use:> sum_{j=1}^M a_j phi(x-x_j)> as an approximation to the surface given an appropriate basis function phi.> The coefcients a_j being determined by the ensuing interpolation> equations:> sum_{j=1}^M a_j phi(x_i-x_j) = f(x_i), i=1,...,M.> There are several popular choices for the basis function. One could use> phi(x)=exp(-|x|^2/2) for example in any space dimension. In 2-dimensions> phi(x)=|x|^2ln|x| is popular and the resulting interpolant is often called> the thin plate spline (TPS). Sometimes (for example with the TPS) a> polynomial term and additional constraints are added to the above to> guarantee the solvability of the system.> RBF Interpolants are a multivariate analogy of the well-loved natural> splines from 1-d.> R-- => Im working on a project that requires sketch a surface knowing several> (rather many but not enough) points on that surface. I need to estimate> other points...>> Some one told me that Fourier analysis can handle this kind of problem.>> Could anyone advice what books or online materials are a good ones to read> about Fourier analysis?>>To get started the most straightforward computational help will be jpegcompression, this is 2D fourier transform then a run length encode (the gifalgorithm). There should be plenty of self help on that with code.Best if you get the data transformed, then visualise to see what approachis best from there. I suspect that several repeating surfaces would workbest with this approach.Herc => In the case of phi(x) = |x|^2ln|x| isnt there a problem when x = x_j> when generating the a_j coefcients? I think you end up with ln(0),> which is not very useful.> How are you supposed to work around this?lim phi(x) = lim phi(x) = 0 as x to 0.Arnold Neumaier =>> In the case of phi(x) = |x|^2ln|x| isnt there a problem when x = x_j> when generating the a_j coefcients? I think you end up with ln(0),> which is not very useful.>> How are you supposed to work around this?>> lim phi(x) = lim phi(x) = 0 as x to 0.>> Arnold NeumaierWhat he said :D =>> In the case of phi(x) = |x|^2ln|x| isnt there a problem when x = x_j> when generating the a_j coefcients? I think you end up with ln(0),> which is not very useful.>> How are you supposed to work around this?>> lim phi(x) = lim phi(x) = 0 as x to 0.>> Arnold Neumaier>> What he said :D>Like Jerry Seineld says, Neumaier!Herc Analysts:I am trying to nd the 3 Dimensional discretization formulae for FTCSmethod applied for diffusion equation.Thus I need formulae for Cartetian, Cylindrical and Spherical coordinatesystem. Can you please suggest me a good reference for these formulae?By Diffusion equation (Heat Equation) I meandu/dt = alpha*Laplacian(u) + source_term. =I retrieved from Netlib the algorithm 661 from TOMS, which is aboutcomputing an interpolant using scattered data, based on an algorithm fromRenka (a modied Shepard method)I compiled the code under windows: - compilation of the fortran code with the G77 fortran compiler fromcygwin. - f2c ed the fortran code and compiled it with the Visual C++ compiler(Visual 7 aka .NET)I have different results in output:G77 compiled fortran code: MAXIMUM ABSOLUTE ERRORS IN THE INTERPOLANT Q AND PARTIAL DERIVATIVES (QX,QY,QZ) RELATIVE TO MACHINE PRECISION EPS FUNCTION MAX ERROR MAX ERROR/EPS Q 0.182E-06 1.53 QX 0.106E-05 QY 0.351E-05 QZ 0.129E-05 ABSOLUTE ERRORS IN THE INTERPOLANT Q AND PARTIAL DERIVATIVES (QX,QY,QZ) RELATIVE TO MACHINE PRECISION EPS FUNCTION MAX ERROR MAX ERROR/EPS Q .165E-06 1.38 QX .163E-05 QY .130E-05 QZ .215E-05Has anyone any idea of what is happening and why ?I dont think it is of big importance, but I am just curious.Pierre-Alain FAYOLLE. => I retrieved from Netlib the algorithm 661 from TOMS, which is about[...]> I compiled the code under windows:> - compilation of the fortran code with the G77 fortran compiler from> cygwin.> - f2c ed the fortran code and compiled it with the Visual C++ compiler> (Visual 7 aka .NET)> I have different results in output:[...]> Has anyone any idea of what is happening and why ?Different optimizations. Pentiums oating point registers have larger accuracies than ordinary double oats stored in memory. Depending on how much of your calculation is done in fp registers (as opposed to from memory) your answer will be slightly more or less precise. Nothing to worry about.Matthias =**** Gambling Dice Machine probability ****A gambling machine:A possiblity table of the possiblity of each side of a dice, from 1 to 6:__________________________________________|side | 1, 2, 3, 4, 5, 6 |------------------------------------------|possiblity | 0.1, 0.2, 0.2,0.15,0.2,0.15 | (total possiblity of one)__________________________________________There are four dices on the screen. __ __ __ __ |__| |__| |__| |__|Each prize draw will roll the dices to show 4 random numbers based on the probability provided in the table.The award prize is the sum of all the numbers appeared in the 4 slots. E.g. 3351 will have a prize of 3+3+5+1=12. So the prize is in the range of 4 to 24.Q1: What is the average prize for this gambling machine?Q2: Generalize the question: n numbers, w_1, w_2, w_3 .... w_n, and their corresponding possiblities: p_1, p_2, p_3 .... p_n. With replacement, 4 boxes, each box lled with one number, what is the average sum of these four numbers? =I have got a formula, but I do think it is the smartest way to do ... I state as following: for a single case: like 3351, for example, it isprobability = 0.2*0.2*0.2*0.1prize = 3+3+5+1= 12weight of 3351 = probability * prize = 12 * (0.2*0.2*0.2*0.1)So the forumula of weight for a single case is: weight = (w_i+w_j+w_k+w_l)*(p_i*p_j*p_k*p_l)Enumerate all the case, from 1111 to 6666, get all the weight, and add them together, then get the average: for i=1~6 for j=1~6 for k=1~6 for l=1~6W_{avg}=sum_{l=1}^{6}sum_{k=1}^{6}sum_{j=1}^{6}sum_{i=1} ^{6}(w_{i}+w_{j}+w_{k}+w_{l}) *(p_{i} * p_{j} * p_{k} * p_{l})But this algorithm will end up with a high complexity, because there is redudency: e.g. 1121 and 1112, in the enumeration, have the same sum, but were counted twice, similarly, 2111, 1211.So I guess that there will be a smarter way than mine... Any suggestion is highly appreciated! => **** Gambling Dice Machine probability ****> A possiblity table of the possiblity of each side of a dice, from 1 to 6:> __________________________________________> |side | 1, 2, 3, 4, 5, 6 |> ------------------------------------------> |possiblity | 0.1, 0.2, 0.2,0.15,0.2,0.15 | (total possiblity of one)> __________________________________________> There are four dices on the screen.> Each prize draw will roll the dices to show 4 random numbers based on > the probability provided in the table.> The award prize is the sum of all the numbers appeared in the 4 slots. > E.g. 3351 will have a prize of 3+3+5+1=12. So the prize is in the range > of 4 to 24.> Q1: What is the average prize for this gambling machine?> Q2: Generalize the question: n numbers, w_1, w_2, w_3 .... w_n, and > their corresponding possiblities: p_1, p_2, p_3 .... p_n. With > replacement, 4 boxes, each box lled with one number, what is the > average sum of these four numbers?Hmm, sounds vaguely like a homework problem. In fact, I might make it a homework problem in my class....Your follow-up post gave a way to compute it that will work, but as you say, will be very slow, and hard to generalize.Have you ever seen the formulaE[X1 + X2 + X3 + X4] = E[X1] + E[X2] + E[X3] + E[X4] ?I think that would help in this situation.Andrew-- Andrew RossIndustrial and Systems EngineeringLehigh University, www.lehigh.edu/~inime/Bethlehem, Pennsylvania, USA =>> **** Gambling Dice Machine probability ****>> A possiblity table of the possiblity of each side of a dice, from 1 to 6:>> __________________________________________>> |side | 1, 2, 3, 4, 5, 6 |>> ------------------------------------------>> |possiblity | 0.1, 0.2, 0.2,0.15,0.2,0.15 | (total possiblity of one)>> __________________________________________>> There are four dices on the screen.>> Each prize draw will roll the dices to show 4 random numbers based on >> the probability provided in the table.>> The award prize is the sum of all the numbers appeared in the 4 >> slots. E.g. 3351 will have a prize of 3+3+5+1=12. So the prize is in >> the range of 4 to 24.>> Q1: What is the average prize for this gambling machine?>> Q2: Generalize the question: n numbers, w_1, w_2, w_3 .... w_n, and >> their corresponding possiblities: p_1, p_2, p_3 .... p_n. With >> replacement, 4 boxes, each box lled with one number, what is the >> average sum of these four numbers? Hmm, sounds \ vaguely like a homework problem. In fact, I might make it a > homework problem in my class....> Your follow-up post gave a way to compute it that will work, but as you > say, will be very slow, and hard to generalize.> Have you ever seen the formula> E[X1 + X2 + X3 + X4] = E[X1] + E[X2] + E[X3] + E[X4] ?> I think that would help in this situation.I am very glad that I could make some contribution to your class.There is a similar problem, I think it might not be the most effective solution as well.... Looks like it is, but not 100% sureProblem statement:A probability table of each side of a dice, from 1 to 6:__________________________________________|side | 1, 2, 3, 4, 5, 6 |------------------------------------------|possiblity | p_1, p_2, p_3, p_4, p_5, p_6 | (total possiblity of one)__________________________________________There are four dices on the screen. __ __ __ __ |__| |__| |__| |__|Each prize draw will roll the dices to show 4 random numbers based on the probability provided in the table.The difference in that, the prize of award of each roll is the maximum number of the four dices:e.gif a roll is 3351, then the prize will be max(3,3,5,1) = 5if a roll is 1121, then the prize will be max(1,1,2,1) = 2Q: What is the average prize of award for this scheme?Solution (maybe not the best):1. caculating the probability of each possible prize, from 1 to 6. P(i),e.g. P(5), is the possiblity that a roll of four dices will have at least one ve, e.g. 1235, 1154, 5554, 1151, thus, prize is 5P(6)P(5)P(4)P(3)P(2)P(1)The caculation of P(5) is:P(5) = (p_5+p_4+p_3+p_2+p_1)4 - (p_4+p_3+p_2+p_1)4Where (p_5+p_4+p_3+p_2+p_1)4 is the probability that for the 4 dices,1,2,3,4,5 could appear in each dice. (e.g the probability of all from 1111 to 5555)Minus the possiblity that there is no ve exists (p_4+p_3+p_2+p_1)4, i.e. the probability of all from 1111 to 4444So it is the case, that at least one ve exist.Generalized formula:P(i)=(sum_(a=1)^(i)p_{a})4 - (sum_{i=b}^{i-1}p_{b})42 Weight each number with its probability, and sum all of them:Average = 6*P(6) + 5*P(5) + 4*P(4) + 3*P(3) + 2*P(2) + 1*P(1)Average Formula = sum_{i=1}^{6} P(i)*iI am not sure is highly appreciated!! =Sorry, the power 4 was not clearly shown in the previous msg. Here is the same post again. Sorry for the redundency...I am very glad that I could make some contribution to your class.There is a similar problem, I think it might not be the most effective solution as well.... Looks like it is, but not 100% sureProblem statement:A probability table of each side of a dice, from 1 to 6:__________________________________________|side | 1, 2, 3, 4, 5, 6 |------------------------------------------|possiblity | p_1, p_2, p_3, p_4, p_5, p_6 | (total possiblity of one)__________________________________________There are four dices on the screen. __ __ __ __ |__| |__| |__| |__|Each prize draw will roll the dices to show 4 random numbers based on the probability provided in the table.The difference in that, the prize of award of each roll is the maximum number of the four dices:e.gif a roll is 3351, then the prize will be max(3,3,5,1) = 5if a roll is 1121, then the prize will be max(1,1,2,1) = 2Q: What is the average prize of award for this scheme?Solution (maybe not the best):1. caculating the probability of each possible prize, from 1 to 6. P(i),e.g. P(5), is the possiblity that a roll of four dices will have at least one ve, e.g. 1235, 1154, 5554, 1151, thus, prize is 5P(6)P(5)P(4)P(3)P(2)P(1)The caculation of P(5) is:P(5) = (p_5+p_4+p_3+p_2+p_1)^4 - (p_4+p_3+p_2+p_1)^4Where (p_5+p_4+p_3+p_2+p_1)^4 is the probability that for the 4 dices,1,2,3,4,5 could appear in each dice. (e.g the probability of all from 1111 to 5555)Minus the possiblity that there is no ve exists (p_4+p_3+p_2+p_1)^4, i.e. the probability of all from 1111 to 4444So it is the case, that at least one ve exist.Generalized formula:P(i)=(sum_(a=1)^(i)p_{a})^4 - (sum_{i=b}^{i-1}p_{b})^42 Weight each number with its probability, and sum all of them:Average = 6*P(6) + 5*P(5) + 4*P(4) + 3*P(3) + 2*P(2) + 1*P(1)Average Formula = sum_{i=1}^{6} P(i)*iI am not sure is it clear to say like appreciated!! =What is the value of this limit =>> What is the value of this limit :>> asymptotic series, eg Handbok of Mathematical Functions (its on line ifyou havent got a copy!)Bill =>>> What is the innity);>> Chris>> 1/2>> Use asymptotic series, eg Handbok of Mathematical Functions (its on lineif> you havent got a copy!)It is interesting that in this case we dont even have to know theasymptotics. It is enough to know that the asymptotics by the followingversion of Laplace method.For unimodal functions f with a sharp maximum at t=t0, fast approaching 0 att->0 and at t-> innity, the integral int(f(t),t=0..innity) can beapproximated by I=int(exp(g(t)),t=-innity..innity) whereg(t)=a-b*(t-t0)^2 is the beginning of the Taylor series of ln f(t) centeredat t=t0. In this case, the integrals int(f(t),t=0..t0) andint(f(t),t=t0..innity) can be well approximated byint(exp(g(t)),t=-innity,t0) and int(exp(g(t)),t=t0..innity). Each ofthem equals the half of I, because exp(g(t)) is an even function of t-t0.Thus, for such functions f, int(f(t),t=t0..innity)/int(f(t),t=0..innity)is close to 1/2, and the better is the approximation ont(f(t),t=0..innity) by I, the closer is that fraction to 1/2.Now, note that f(t)=t^(x-1)*exp(-t) is achieving maximum at t0=x-1 ,t=t0..innity),so the limit in question is the limit ont(f(t),t=t0..innity)/int(f(t),t=0..innity) which is 1/2 because themethod described above leads to the Stirlings Mihailovshttp://webpages.shepherd.edu/amihailo/ =How do you interpolate smoothly (i.e. with continuous rst derivative) when the sample points are scattered randomly in two or more dimensions?-- Anton Sherwood, http://www.ogre.nu/ => How do you interpolate smoothly (i.e. with continuous rst derivative)> when the sample points are scattered randomly in two or more dimensions?What you want is called scattered data interpolation;many links are given in http://www.mat.univie.ac.at/~neum/surface.html[~ is a tilde]Arnold Neumaier =>How do you interpolate smoothly (i.e. with continuous rst derivative) >when the sample points are scattered randomly in two or more dimensions?>you can do this in quite different techniques: global interpolation orpiecewise polynomial. global techniques dene the interpolating functionas a linear combination of functions dened everywhere (e.g. radial basisfunctions). piecewise polynomial functions mimic whta is known from splines in 1D.you need then a triangulation (tedrahedral decomposition in 3d)here are some sourceshope this helps peter http://www.netlib.org/toms toms/684 gams: E2b for: C1 and C2 interpolation on triangles with quintic and nonic bivariate polynomials by: A. Preusser ref: ACM TOMS 16 (1990) 253-257 interpolation, piecewise polynomial interpolation gams: E2b title: IDBVIP and IDSFFT for: bivariate interpolation and smooth surface tting for irregularly distributed data points by: H. Akima ref: ACM TOMS 4 (1978) 160-164 toms/660 gams: E2b title: QSHEP2D for: quadratic Shepard method for bivariate interpolation of scattered data by: R. J. Renka ref: ACM TOMS 14 (1988) 149-150 toms/661 gams: E2b title: QSHEP3D for: quadratic Shepard method for trivariate interpolation of scattered data by: R. J. Renka ref: ACM TOMS 14 (1988) 151-152 TOMS 25,1 (Mar 1999) 70 for: {CSHEP2D}: Cubic {Shepard Method for Bivariate Interpolation of Scattered Data by: R. J. Renka size: 103 kB 1999) 74 for: {TSHEP2D}: Cosine Series {Shepard} Method for Bivariate Interpolation of Scattered Data by: R. J. Renka and R. Brown size: 106 kB (Mar 1999) 78 for: Accuracy Tests of {ACM} Algorithms for Interpolation of Scattered Data in the Plane by: R. J. Renka and R. Brown size: 96 kB = > . . . . piecewise polynomial functions mimic > whta is known from splines in 1D. > you need then a triangulation (tedrahedral decomposition in 3d)Delaunay, I suppose? > here are some sources > hope this helps > http://www.netlib.org/toms-- Anton Sherwood, http://www.ogre.nu/ => How do you interpolate smoothly (i.e. with continuous rst derivative) > when the sample points are scattered randomly in two or more dimensions?Generally speaking, the constrained minimization problemIntegral F[Some functional of the curvature tensor of a surface] du dvsubject toG(uk,vk)is unsolved and this eld is wide open.Except for no G constraints and F=K_Gauss in which case the solutionis a minimal surface.For practical use I have found [professor Richard] Frankes LTPS[Local Thin Plate Splines] method quite useful.The method is based on Duchon(?)s Thin Plate Splines which interpolatethe data points but have a logarithmic term and which are C1 continuousand C2 except at the data points.Frankes LOCAL Thin Plate Spline blends local thin plate splinesinto a continuous interpolant by using cartesian products of Hermiteblending functions which are also C1 and C2 except at the patch lines.Yes, the original code is in Fortran. I have a C version where I haveremoved the special cases for the boundary regions and improved thehandling of empty regions and which I have tied into a contouringmethod.The results are quite satisfying for small number of data points (20-100) but less so for large datasets (5000) although the methodworks.Jentje Goslinga => What you want is called scattered data interpolation;> many links are given in > http://www.mat.univie.ac.at/~neum/surface.html(You know the very rst link on that page - POVRAY 2.2 with equipotential surface - is dead? PoVRay is up to 3.5 now anyway.)-- Anton Sherwood, http://www.ogre.nu/ =I need to generate a random orthogonal basis of dimension N.I know, I can use a regular random number generator, to generate NxN numbersand then orthogonalize the N vectors. To do that however, I will have tostore NxN numbers. N gets really large for my problem.I wonder if there is a random number generator out there, that will generatenumbers in sequences, each of size N, such that all the sequences areautomatically orthogonal to each others.Alien+ values ofz of the shape 1/2+it or 1/4+it, and t can besmall. I need gmp multiprecision (which I usethrough mpfr). So Lanczos formula seemsto be what I need, but if I could nd in advance, Amities, Olivier =I read recently about these two largest named numbers, and eversince have been wondering how much bigger Grahams Number is thanSkewes.Is there any way of expressing how many times the one can go into theother? =Are there any good hash functions for permutations?By good I mean: 1. Each permutation is uniquely mapped to a small number. 2. The hash of a permutation can be quickly computed.What I need to do is select 1000 distinct random permutationsfrom the set of n! possible permutations where n >= 7. My algorithm for doing this:1. Select a random permutation.2. Has this permutation already been selected? If so, repeat. Since I need to compare the currently selected permutationwith all of the previously selected permutations, I am looking for a fast way to compare. Of course, the worst case in comparing 2 permutations is (n - 1) comparisons. Obviously, I can afford to live with a few hash collisions.-- K. Banerjee =Fike, A permutation generation method The Computer Journal, 18-1, Feb 75, 21-22.> Are there any good hash functions for permutations?> By good I mean: > 1. Each permutation is uniquely mapped to a small number. > 2. The hash of a permutation can be quickly computed.> What I need to do is select 1000 distinct random permutations> from the set of n! possible permutations where n >= 7. My > algorithm for doing this:> 1. Select a random permutation.> 2. Has this permutation already been > selected? If so, repeat. > Since I need to compare the currently selected permutation> with all of the previously selected permutations, I am looking > for a fast way to compare. Of course, the worst case in > comparing 2 permutations is (n - 1) comparisons. > Obviously, I can afford to live with a few hash collisions.> -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. ---Randomness comes in bunches. =Knuth describes a cute algorithm for generating a randompermutation. I dont remember if it is his own or not. I have nottested the pseudo-code below.for (i=1; i<=n; i++) a[i] = i;for (i=1; i<=n; i++) { j = random((i+1)..n); swap a[i] and a[j];}> Are there any good hash functions for permutations?> By good I mean: > 1. Each permutation is uniquely mapped to a small number. > 2. The hash of a permutation can be quickly computed.> What I need to do is select 1000 distinct random permutations> from the set of n! possible permutations where n >= 7. My > algorithm for doing this:> 1. Select a random permutation.> 2. Has this permutation already been > selected? If so, repeat. > Since I need to compare the currently selected permutation> with all of the previously selected permutations, I am looking > for a fast way to compare. Of course, the worst case in > comparing 2 permutations is (n - 1) comparisons. Obviously, I can \ afford to live with a few hash collisions.> -- > K. Banerjee-- Use of tools distinguishes Man from Beast. And UNIX users from WINDOZE lusers. => Are there any good hash functions for permutations?> By good I mean:> 1. Each permutation is uniquely mapped to a small number.> 2. The hash of a permutation can be quickly computed.> What I need to do is select 1000 distinct random permutations> from the set of n! possible permutations where n >= 7. My> algorithm for doing this:> 1. Select a random permutation.> 2. Has this permutation already been> selected? If so, repeat.> Since I need to compare the currently selected permutation> with all of the previously selected permutations, I am looking> for a fast way to compare. Of course, the worst case in> comparing 2 permutations is (n - 1) comparisons.> Obviously, I can afford to live with a few hash collisions.> --> K. BanerjeeNot a hash function, but something that might be useful:Ranking function for permutations,i.e. assigning a unique integer in the range [0,n!-1] to each of theRanking and Unranking Permutations in Linear Time, available athttp://www.csr.uvic.ca/~fruskey/Publications/ RankPerm.htmlBTW, Knuths draft edition ofTAOCP Vol.4, Pre-fascicle 2b (generating all permutations) is stillavailable for download: http://www-cs-faculty.stanford.edu/~knuth/fasc2b.ps.gzHTHHugo =>> Are there any good hash functions for permutations?>> By good I mean:>> 1. Each permutation is uniquely mapped to a small number.> 2. The hash of a permutation can be quickly computed.>> What I need to do is select 1000 distinct random permutations> from the set of n! possible permutations where n >= 7. My> algorithm for doing this:>> 1. Select a random permutation.> 2. Has this permutation already been> selected? If so, repeat.> Since I need to \ compare the currently selected permutation with all of the previously \ selected permutations, I am looking> for a fast way to compare. Of course, the worst case in> comparing 2 permutations is (n - 1) comparisons.> Obviously, I can \ afford to live with a few hash collisions.>>> -->> K. Banerjee> Not a hash function, but something that might be useful:> Ranking function for permutations,> i.e. assigning a unique integer in the range [0,n!-1] to each of the> Ranking and Unranking Permutations in Linear Time, available at> http://www.csr.uvic.ca/~fruskey/Publications/RankPerm.html>[. ..] Frank Ruskey sent me the following note:<>Hugo Pfoertner =I seem to have some sort of a mental block and cant quite see how toimplement numerical integration in a dynamic simulation I am working on.The problem is one where I can establish initial value forces, and know themass, hence expect to derive an acceleration at each time step from F=Ma, andthen was planning to use the 4th order Runge-Kutta method (RKG4) to derivevelocity, and then again to derive displacment. In looking at a fewreferences re RKG4 (and other numerical methods) the problem is usuallydescribed in terms of a function of y and t ie f(t,y) and the RKG4 solution,using time steps of h, is described (generally) as; m1 = f(t,y) m2 = f(t + 0.5*h, y + 0.5*h*m1) m3 = f(t + 0.5*h, y + 0.5*h*m2) m4 = f(t + h, y + h*m3) t = t + h y = y + h*(m1 + 2*(m2 + m3) + m4)/6 Now I know Im getting old and maybe a tad silly, but I dont quite see howto apply this logic to my problem, which doesnt have a nice function, simplyan initial acceleration values at time t0 (say). The forces are derived fromquite complex relationships which I would nd impossible to describe interms of time.Some very brief background, in case it helps understanding. The dynamicsimulation is that of a motor vehicle where the engine torque-rpm curve is(sort of) arbitrary, the driver logic is also a bit that way, the tractionsurface is quite variable, as is the path and hence the resistance forces,but given the start conditions, a particular driver logic, and a 3d path, itis possible to derive the forces that are being applied to the vehicle. So the question is, to anyone nice enough to try to help, how to I apply RKG4in a situation like this? Is there a better way of doing it?-- Terry Duell => I seem to have some sort of a mental block and cant quite see how to> implement numerical integration in a dynamic simulation I am working on.> The problem is one where I can establish initial value forces, and know the> mass, hence expect to derive an acceleration at each time step from F=Ma, and> then was planning to use the 4th order Runge-Kutta method (RKG4) to derive> velocity, and then again to derive displacment. In looking at a few> references re RKG4 (and other numerical methods) the problem is usually> described in terms of a function of y and t ie f(t,y) and the RKG4 solution,> using time steps of h, is described (generally) as;> m1 = f(t,y)> m2 = f(t + 0.5*h, y + 0.5*h*m1)> m3 = f(t + 0.5*h, y + 0.5*h*m2)> m4 = f(t + h, y + h*m3)> t = t + h> y = y + h*(m1 + 2*(m2 + m3) + m4)/6> Now I know Im getting old and maybe a tad silly, but I dont quite see how> to apply this logic to my problem, which doesnt have a nice function, simply> an initial acceleration values at time t0 (say). The forces are derived from> quite complex relationships which I would nd impossible to describe in> terms of time.> Some very brief background, in case it helps understanding. The dynamic> simulation is that of a motor vehicle where the engine torque-rpm curve is> (sort of) arbitrary, the driver logic is also a bit that way, the traction> surface is quite variable, as is the path and hence the resistance forces,> but given the start conditions, a particular driver logic, and a 3d path, it> is possible to derive the forces that are being applied to the vehicle.> So the question is, to anyone nice enough to try to help, how to I apply RKG4> in a situation like this? Is there a better way of doing it?IMO your problem is the second order of your differential equation.Since fromF(t)= m * a(t) = m * dv/dt = m * d^2s/dt^2 (mass m is constant)you haved^2s/dt^2 = F(t)/m .In addition, every ordinary differential equation of second order can betreated as a system of two ordinary differential equations of rst order, henceds/dt = v(t)dv/dt = F(t)/mMeans, you might apply Runge-Kutta to this set of equations, which suitquite well to the formulation of the RK algorithm. In practice, youcompute s(t) and v(t) synchronously by[formula above]s = s + h*(m1_s + 2*(m2_s + m3_s) + m4_s)/6v = v + h*(m1_v + 2*(m2_v + m3_v) + m4_v)/6Axel-- Dr. rer.nat. Axel HuttWeierstra-Institute for Applied Analysis and StochasticsMohrenstr. 39, 10117 Berlin, Germanyhttp://www.wias-berlin.de/people/hutt => Now I know Im getting old and maybe a tad silly, but I dont quite see how> to apply this logic to my problem, which doesnt have a nice function, simply> an initial acceleration values at time t0 (say). The forces are derived from> quite complex relationships which I would nd impossible to describe in> terms of time.> Some very brief background, in case it helps understanding. The dynamic> simulation is that of a motor vehicle where the engine torque-rpm curve is> (sort of) arbitrary, the driver logic is also a bit that way, the traction> surface is quite variable, as is the path and hence the resistance forces,> but given the start conditions, a particular driver logic, and a 3d path, it> is possible to derive the forces that are being applied to the vehicle.> So the question is, to anyone nice enough to try to help, how to I apply RKG4> in a situation like this? Is there a better way of doing it?Ive formulated several models regarding e.g. combustion engines. The secret is to model all the small parts separately with input and output signals and couple them in a system of differential equations between various variables. Formulate the system as a state-space and then just integrate. This is very easy in e.g. MATLAB Simulink (expensive), or Scilabs Scicos, which is free.Otherwise youll have to derive it all by hand and then apply an ODE-solver. There are several good solvers out there already, so there shouldnt be any need for writing your own.-- /S => I seem to have some sort of a mental block and cant quite see how to> implement numerical integration in a dynamic simulation I am working on.> The problem is one where I can establish initial value forces, and know the> mass, hence expect to derive an acceleration at each time step from F=Ma, and> then was planning to use the 4th order Runge-Kutta method (RKG4) to derive> velocity, and then again to derive displacment. In looking at a few> references re RKG4 (and other numerical methods) the problem is usually> described in terms of a function of y and t ie f(t,y) and the RKG4 solution,> using time steps of h, is described (generally) as;> m1 = f(t,y)> m2 = f(t + 0.5*h, y + 0.5*h*m1)> m3 = f(t + 0.5*h, y + 0.5*h*m2)> m4 = f(t + h, y + h*m3)> t = t + h> y = y + h*(m1 + 2*(m2 + m3) + m4)/6> Now I know Im getting old and maybe a tad silly, but I dont quite see how> to apply this logic to my problem, which doesnt have a nice function, simply> an initial acceleration values at time t0 (say). The forces are derived from> quite complex relationships which I would nd impossible to describe in> terms of time.> Some very brief background, in case it helps understanding. The dynamic> simulation is that of a motor vehicle where the engine torque-rpm curve is> (sort of) arbitrary, the driver logic is also a bit that way, the traction> surface is quite variable, as is the path and hence the resistance forces,> but given the start conditions, a particular driver logic, and a 3d path, it> is possible to derive the forces that are being applied to the vehicle.> So the question is, to anyone nice enough to try to help, how to I apply RKG4> in a situation like this? Is there a better way of doing it?Runge-Kutta, euler or heun solve diferential equations likedy/dx(xi)=f(xi,yi), so with a initial conditions you can solve forevery point. And i think that u have not a difencial eq.Pablo =Terry, by denition any dynamics problem can be expressed in terms ofrst-order differential equationsdy/dt = f(t, y)Your function need not be a function only of t, but also (and typically)of y. Because dynamics is basicallya second-order description (based on acceleration), we typically have towritey = [x, v]where x is all the position-like parameters, and v all thevelocity-like. The derivative of x is v (orexpressible in terms of v), and the derivative of v is youracceleration. You need to write f(t,y) so thatgiven test values of t and y, it returns with dy/dt.A word of warning: Dont try to separate the integration of the positionand velocity components. It soundslike thats what youre contemplating, but that doesnt work. f(t,y)must compute all the derivatives simultaneously. The integrator will then properly integrate themsimultaneously.Jack> I seem to have some sort of a mental block and cant quite see how to> implement numerical integration in a dynamic simulation I am working on.> The problem is one where I can establish initial value forces, and know the> mass, hence expect to derive an acceleration at each time step from F=Ma, and> then was planning to use the 4th order Runge-Kutta method (RKG4) to derive> velocity, and then again to derive displacment. In looking at a few> references re RKG4 (and other numerical methods) the problem is usually> described in terms of a function of y and t ie f(t,y) and the RKG4 solution,> using time steps of h, is described (generally) as;> m1 = f(t,y)> m2 = f(t + 0.5*h, y + 0.5*h*m1)> m3 = f(t + 0.5*h, y + 0.5*h*m2)> m4 = f(t + h, y + h*m3)> t = t + h> y = y + h*(m1 + 2*(m2 + m3) + m4)/6> Now I know Im getting old and maybe a tad silly, but I dont quite see how> to apply this logic to my problem, which doesnt have a nice function, simply> an initial acceleration values at time t0 (say). The forces are derived from> quite complex relationships which I would nd impossible to describe in> terms of time.> Some very brief background, in case it helps understanding. The dynamic> simulation is that of a motor vehicle where the engine torque-rpm curve is> (sort of) arbitrary, the driver logic is also a bit that way, the traction> surface is quite variable, as is the path and hence the resistance forces,> but given the start conditions, a particular driver logic, and a 3d path, it> is possible to derive the forces that are being applied to the vehicle.> So the question is, to anyone nice enough to try to help, how to I apply RKG4> in a situation like this? Is there a better way of doing it?> --> Terry Duell =still ghting with it.Now my question is as follow.I have a 2x2 matrix A with determinant equal to 1 a_11 = 1 a_12 = -1/k a_21 = m f*f a_22 = 1-m*f*f/kEvething is real.I would like to compute the determinant as a function of f.Now the problem arise when I compute the determinant of A^6which again is always one, but if I compute it explicitly i get nonsensewhen f is too large.Any idea ?(I need to compute such kind of determinants because in some cases it isno more one).Luciano--Luciano RibichiniMax Planck Insitute for Gravitational PhysicsCallinstr. 38D 30167 Hannover Germany => still ghting with it.> Now my question is as follow.> I have a 2x2 matrix A with determinant equal to 1> a_11 = 1> a_12 = -1/k> a_21 = m f*f> a_22 = 1-m*f*f/k> Evething is real.> I would like to compute the determinant as a function of f.> Now the problem arise when I compute the determinant of> A^6 which again is \ always one, but if I compute it explicitly i get nonsense> when f is too large.> Any idea ?> (I need to compute such kind of determinants because in some cases it is> no more one).> Luciano> --> Luciano Ribichini> Max Planck Insitute for Gravitational Physics> Callinstr. 38> D 30167 Hannover Germany det( A^6) == [ det(A) ] ^6If you do it by explicitly raising the matrix to the power 6, you will en-counter roundoff problems. Try it for A^2 and see what you have to compute.Lots of small differences between large numbers, if mf^2/k is large. Indet( A^2 ) I see this term to the 3rd power being added and subtractedto/from unity. But (say m=k=1) if f = 1 + (max_fp_#)^{1/6}, then f^6overows. It gets worse as the powers increase. For HW gure out thelargest power of f in det(A^6) and then ask how large f can be beforethe calculation overows on your machine. It might surprise you. What every computer scientist should know about oating point arithmetic by D. Goldberg.You can download it in *.pdf format from http://www.phys.virginia.edu/classes/551.jvn.fall01/-- Julian ^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. => still ghting with it.> Now my question is as follow.>> I have a 2x2 matrix A with determinant equal to 1> a_11 = 1> a_12 = -1/k> a_21 = m f*f> a_22 = 1-m*f*f/k>> Evething is real.> I would like to compute the determinant as a function of f.> Now the problem arise when I compute the determinant of> A^6>> which again is always one, but if I compute it explicitly i getnonsense> when f is too large.>> Any idea ?> (I need to compute such kind of determinants because in some casesit is> no more one).>IIRC, for square A(nxn), det(A) = Product (eigenvalue_i). In the caseof the posted A, det(A)=1 as one can explicitly demonstrate. I gatherthat youre calculating det(A) using an algorithm that involves powersof A. As Professor Julian correctly points out, this can bebook, Accuracy and chapter on matrix power computations.HTH,Gerry T. => still ghting with it.> Now my question is as follow.> I have a 2x2 matrix A with determinant equal to 1> a_11 = 1> a_12 = -1/k> a_21 = m f*f> a_22 = 1-m*f*f/k> Evething is real.> I would like to compute the determinant as a function of f.> Now the problem arise when I compute the determinant of> A^6> which again is always one, but if I compute it explicitly i get nonsense> when f is too large.> Any idea ?> (I need to compute such kind of determinants because in some cases it is> no more one).> Luciano> --> Luciano Ribichini> Max Planck Insitute for Gravitational Physics> Callinstr. 38> D 30167 Hannover GermanyWhat is the reason to rst compute A^6 and then det(A^6)?All software implementations of det(A) have some means tocircumvent the overow problem, see e.g.Class D3a1 ... Determinants of general real nonsymmetric matriceshttp://gams.nist.gov/serve.cgi/Class/D3a1/http:// gams.nist.gov/serve.cgi/Module/NAPACK/DET/11837/ andhttp://gams.nist.gov/serve.cgi/ModuleComponent/11837/ Fullsource/NETLIB/det.fOf course you must factor the matrix before computing the determinant:http://gams.nist.gov/serve.cgi/ModuleComponent/ 11854/Fullsource/NETLIB/fact.fHugo Pfoertner =Consider an unknown matrix A (NxN), where yn=A*xn, xn input vector, and ynthe output vector.A is unknown, but the user can supply xn to a black box and obtain yn byforward calculations.My target is to estimate the principal left singular vectors of the matrixA.principal, i.e. corresponding to those with relatively large singularvalues.The direct approach to do that, is to obtain a set of N pairs (xn,yn) andrequire that xn are independent, and hence one can back-calculate A.However, I wonder if there is another way that would reduce the number offorward calculations << N, since I dont need the whole spectrum of leftsingular vectors, and in addition A is severelly ill-conditioned, so inprincipal r<Consider an unknown matrix A (NxN), where yn=A*xn, xn input vector, and yn >the output vector. >A is unknown, but the user can supply xn to a black box and obtain yn by >forward calculations. >My target is to estimate the principal left singular vectors of the matrix >A. >principal, i.e. corresponding to those with relatively large singular >values. >The direct approach to do that, is to obtain a set of N pairs (xn,yn) and >require that xn are independent, and hence one can back-calculate A. >However, I wonder if there is another way that would reduce the number of >forward calculations << N, since I dont need the whole spectrum of left >singular vectors, and in addition A is severelly ill-conditioned, so in >principal r< valuesfor this a matrix-vector-multiply routine is enough and you have one.look athttp://www.netlib.org/svdpackit does exactly what you wantpeter =I have to inverse a big ill conditioned matrix. So I would like to know if it exists numerical software working in quadruple precision.Karim. =An immediate question that must come here is: do you really need theinverse? What will you do with it? Solving a system of linear equations? Howbig is your matrix, what is its condition number?Zdenek HurakKarim Bernardet p.92se v diskusn.92m pr.92spevku>> I have to inverse a big ill conditioned matrix. So I would like to know it> exists numerical software working in quadruple precision.>> Karim.>---Odchoz.92 zpr.87va neobsahuje viry.Zkontrolov.87no antivirovm syst.8emem AVG (http://www.grisoft.cz).Verze: 6.0.489 / Virov.87 b.87ze: 288 =>An immediate question that must come here is: do you really need the>inverse? What will you do with it? Solving a system of linear equations? How>big is your matrix, what is its condition number?>>Zdenek Hurak>>Karim Bernardet p.92se v diskusn.92m pr.92spevku>I have to inverse a big ill conditioned matrix. So I would like to know if>>itexists numerical software \ working in quadruple precision.>>Karim.>>--->Odchoz.92 zpr.87va neobsahuje viry.>Zkontrolov.87no antivirovm syst.8emem AVG (http://www.grisoft.cz).>Verze: 6.0.489 / Virov.87 b.87ze: inverse it using scalapack (Ax=b with different b corresponding to the identity matrix) but it says that double precision is not enough. The size of matrix is about 8000 and the condition number is very big.I now try to regularize it using SVD decomposition + ltering.Karim. = >I have to inverse a big ill conditioned matrix. So I would like to know if it >exists numerical software working in quadruple precision. >rst you should ask yourself why this is so:is your original problem well dened with an unique solution?is your problem setup o.k.?maybe you use a hamsted formulation:e.g. you want to t a parabola to data, use the forma0+a1*x+a2*x^2, but your data x(i) are in the range 1000000then using the normal equations you end up with a singular matrix (conditionnumber 10^24) . had you used a0+a1*(x-1000000)+a2*(x-1000000)^2 .......all these questions because it is quite seldom that you have exact input data matrix A of such an illconditioning and need now the inverse ...(if your matrix has a single roundoff in it, what should quadruple precisionhelp with computing an inverse which will be totally different from what youexpect?)if nothing of this applies:http://www.netlib.org/toms:le toms/719ref TOMS 19,3for multiprecision translation and execution of Fortran programsby Bailyou can feed a normal fortran code for matrix inversion in, it translates itto a multiple precision version and executes it automatically.hope this helpspeter =http://crd.lbl.gov/~dhbailey/mpdist/mpdist.html =Im trying to gure out a way to plot (in my Delphi 6 programmedfunction-plotting program) the graph of any implicit function, i.e. thegraph of x - 2xy = 0.I dont want to substitute the coordinates of all the visible points in thexoy-grid to see if there should be drawn a point or not ;)Real programs such as Mathematica can draw the graphs of implicitfunctions. Id like to know how this is done. Is there a known algorithm?(prefereably in Delphi/Pascal).TIAHans Roos. =>> Im trying to gure out a way to plot (in my Delphi 6 programmed> function-plotting program) the graph of any implicit function, i.e. the> graph of x - 2xy = 0.> I dont want to substitute the coordinates of all the visible points in the> xoy-grid to see if there should be drawn a point or not ;)> Real programs such as Mathematica can draw the graphs of implicit> functions. Id like to know how this is done. Is there a known algorithm?> (prefereably in Delphi/Pascal).>> TIA>> Hans Roos.I have because oferrors. Hopefully, the clutter will not appear in the newsgroup. This is acorrected response.Anyway, one method of plotting implicit curves is to form a parametricrepresentation in terms of an auxiliary variable, usually s (distance) or t(time). On my website you will nd a paper describing a method for scanconverting parametric curves: http://www.crbond.com/graphic.htmJust scan down the list for the algorithm for scan converting parametriccurves.HTH--There are two things you must never attempt to prove: the unprovable -- and theobvious.--Democracy: The triumph of popularity over principle.--http://www.crbond.com => Im trying to gure out a way to plot (in my Delphi 6 programmed> function-plotting program) the graph of any implicit function, i.e. the> graph of x - 2xy = 0.> Anyway, one method of plotting implicit curves is to form a parametric> representation in terms of an auxiliary variable, usually s (distance)or t> (time). On my website you will nd a paper describing a method for scan> converting parametric curves: http://www.crbond.com/graphic.htmCan every Implicit representation be expressed in a parametricrepresentations ?I am wondering if in some cases you will no nd some ambiguous cases, forwhich the parametric representations is difcult to extract ? expressed in a parametric> representations ?> I am wondering if in some cases you will no nd some ambiguous cases, for> which the parametric representations is difcult to extract ?Well, I only describe one method for plotting implicit functions -- a ratherlarge class of problems can be treated this way. If the method doesnt apply,another method may be needed.--There are two things you must never attempt to prove: the unprovable -- and theobvious.--Democracy: The triumph of popularity over principle.--http://www.crbond.com => ==> Find the greatest positive integer p and smallest integer q , > such that> (1+ 1/pi)^{pi + p/3} < pi < (1+ 1/pi)^{pi + q/2} .> It is true that (p,q) in { (3,4) , (3,5), (4,5) } ?> ==According to Maple, (1+ 1/pi)^{pi + p/3} = pi for p=3.002, so p=3 isyour answer for p. And pi = (1+ 1/pi)^{pi + q/2} for q=2.001, so q=3also.(1+1/Pi)^(Pi+1) = 3.140968877, p=3 works but(1+1/Pi)^(Pi+4/3) = 3.444050170, p=4 is too big;(1+1/Pi)^(Pi+3/2) = 3.606387487, q=3 works but(1+1/Pi)^(Pi+1) = 3.140968877, q=2 is too small. smallest integer q , such that (1+ 1/pi)^{pi + p/3} < pi < (1+ 1/pi)^{pi + q/2} .It is true that (p,q) in { (3,4) , (3,5), (4,5) } ? ==Ive been inverting a particular 18x18 matrix by SVD. As acheck, I gured yesterday that Id give LU decomposition a crackat it and see that the two agreed.They didnt. A closer look showed that the SVD routine wasgiving a good inverse, while LU decomposition wasnt. An evencloser look, and it appears that LU decomposition itself isfailing after a few (say, ve) rows.Im using the implementation of LU decomposition with partialpivoting from Numerical Recipes. The matrix is not particularlysingular (singular values range over about two orders ofmagnitude), and anyway its my understanding that inversion by LUdecomposition is supposed to work *particularly* well on nearly-singular matrices.I understand that, without partial pivoting, the algorithm isextremely unstable, but it surprises me that it would be sounstable with a relatively moderate matrix. Is this a knownproblem with LU decomposition? ------------------------------------------------------------- ---- . . . Except when they dont, Because sometimes they wont. - Dr. Seuss-------------------------------------------------------- ---------Jason Cooper jcooper@acs.ucalgary.ca =from your vague description I guess that there is something goingwrong either in the LU decomposition code or in your calling program.I would suggest to use LAPACK (it has also to be called correctly, ofcourse) or use something like Matlab, Mathematica, Maple, ... andcompare the results.Alois => r*e^(-r*|x|)*p(x,t) (1)> for large r by spatial and temporal integration and I tried> to do the integration via an adaptive Gauss-Kronrod 21-point> integration rule with implied extrapolation (GSL routine> gsl_integration_qags).> It turns out, that everything works ne in case of a smooth> function p(x,t). However, if> p(x,t)=V(x,t) ,> I get roundoff errors or others most of the time. I am> confused about this difference to a smooth function for> p(x,t). My suspicion is, that one reason is the discreteness> of V(x,t), which is unfortunate for the extrapolation steps in> qags. Is this the reason? Are there other routines suitable> for solving (1) ?> Axel> Something seems quite wrong in the way you have posed the problem.What is the variable of integration? Where is the dependence on x?What are the boundary conditions?I suspect you mean something like dV(x,t)/dt = int_a^b dx { r*e^(-r*|x-x|)*p(x,t) } .When p(x,t) is an arbitrary function there is no problem. However,when p = V (the unknown function) you have an integro-differentialequation. For large r, the kernel is almost a delta-function, withsupport only for x approx x, so that the equation becomes dV(x,t)/dt approx const. * V(x,t)which is just a simple exponential with respect to time.-- Julian V. NobleProfessor Emeritus of ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. => Something seems quite wrong in the way you have posed the problem.> What is the variable of integration? Where is the dependence on x?> What are the boundary conditions?> I suspect you mean something like> dV(x,t)/dt = int_a^b dx { r*e^(-r*|x-x|)*p(x,t) } .that is my problem, sorry for my loose formulation.> When p(x,t) is an arbitrary function there is no problem. However,> when p = V (the unknown function) you have an integro-differential> equation. For large r, the kernel is almost a delta-function, with> support only for x approx x, so that the equation becomes> dV(x,t)/dt approx const. * V(x,t)> which is just a simple exponential with respect to time.well, that is trivial. My problem is the integration of an integro-differential equation with nite r, that does not approximates the delta-distribution quite well.Axel =due to the kind reply of J.Noble, I realized toformulate my problem more precisely:I have problems with an integro-differential equationdV(x,t)/dt = int_a^b r*e^(-r*|x-x|)*V(x,t) dx ,which I tried to solve by a method probably assuminga smooth function. The boundary conditions might be periodicor natural.Does anybody know any package/piece of code to solvethis equation in an optimized way? I have already tried ahand-made algorithm applying a uniform discretization, butthat is very slow indeed.Axel solvedV(x,t)/dt=int_a^b r*e^(-r*|x|)*p(x,t) (1)for large r by spatial and temporal integration and I triedto do the integration via an adaptive Gauss-Kronrod 21-pointintegration rule with implied extrapolation (GSL routinegsl_integration_qags).It turns out, that everything works ne in case of a smoothfunction p(x,t). However, ifp(x,t)=V(x,t) ,I get roundoff errors or others most of the time. I amconfused about this difference to a smooth function forp(x,t). My suspicion is, that one reason is the discretenessof V(x,t), which is unfortunate for the extrapolation steps inqags. Is this the reason? Are there other routines suitablefor solving (1) ?Axel-- Dr. rer.nat. Axel HuttWeierstra-Institute for Applied Analysis and StochasticsMohrenstr. 39, 10117 Berlin, Germanyhttp://www.wias-berlin.de/people/hutt =I have a bit of a problem in trying to analyze some data:I have several hundred sets of around 4 data points each.Every point in every set should conform to some function y=f(x)+C.However, in each set, C will be different.Is there some way of nding what the f(x) part is for all of the data points and somehow ignoring the C difference?The most important part of this analysis is to nd the x values of the local and global maximums and minimums, the y values are less important.Ed => I have several hundred sets of around 4 data points each.> Every point in every set should conform to some function y=f(x)+C.> However, in each set, C will be different.> Is there some way of nding what the f(x) part is for all of the data> points and somehow ignoring the C difference?Choose a set of basis functions phi_l(x) and solve the least square problem for y_ik = sum_l phi_l(x_ik) a_l + c_k for all i,k (where data with the same C have the same index k)for the unknowns a_l and c_k. You need to make sure that no linear combination of the Phi_l(x) is constant, to have all coefcients identiable. And the total numberof data points should be much larger than the number of coefcientsa_l, c_k you want to estimate.Arnold Neumaier = >I have a bit of a problem in trying to analyze some data: I have several \ hundred sets of around 4 data points each. >Every point in every set should conform to some function y=f(x)+C. >However, in each set, C will be different. >Is there some way of nding what the f(x) part is for all of the data >points and somehow ignoring the C difference? >The most important part of this analysis is to nd the x values of the >local and global maximums and minimums, the y values are less important. >Ed >do you have an idea concerning f(x)? around 4 data points this calls not for a comlicated thing, maybe a polynomial of second order?then you can do the following:for every data set (x(i,j),y(i,j)) j=1 to number of points in set i, i=number of setyou dene the residuals:res(i,j)=y(i,j)- ( a0+a1*(x(i,j)-xm(i))+a2*(x(i,j)-xm(i))^2+c(i)) with xm(i)=mean value of the x(i,j) in set i c(i) constant for set i but the same a0,a1,a2 for all sets i .now determine a0,a1,a2 and the c(i) such that the sum of squares of all residuals is simultaneously minimized with respect to a0,a1,a2 and the c(i).this gives a large linear least squares problem where the matrix is verysparse but has three dense columns (the columns corresponding to thevariables a0,a1,a2) by an adapted version of the QR solution method you can take advantage of this (make the dense columns the last ones.apply for every subset corresponding to one small dataset the householder transformation transforming the vector of all ones(related to the constant c(i))(as long as there are points in this set) to the rst unit vector. (of course, the same tranformation has to be applied also to the corresponding part of the dense columns and to the y-vector. then renumber the total of rows such that you obtain in the left upperpart the unit matrix, with three dense columns appended. then nish the qr-decomposition by transforming the remaining three dense columns and backsubstitute.hope this helpspeter =STPT is an invaluable tool if you need to import or extract data from anytext (ASCII) le. STPT makes dealing with ASCII data easy and efcientthrough a patented model and a graphical user interface where parameters canbe dened and the extraction operation tested. The data which may benumeric (tables, lists, numbers) or non-numeric (words, sentences) can beimported as an array into languages such as MATLAB, Visual Basic, C or C++where it can be used in calculations. STPT can also be used to automaticallypopulate tables in MS Word, worksheets in MS Excel or databases in MSAccess.A free 10-day trial can be found at http://www.stiwww.com/Parser.htm =Folks:I am sure you have seen this before:I have a tet. in 3d given by 4 nodes P1-P4. I am also gioven a point Pin 3d, and I want to nd ifP is inside this tet. Right now, I just want something that will alwayswork, rather than the most efcient algorithm.Suresh--Its all there when you look,Into Santas book. =>I have a tet. in 3d given by 4 nodes P1-P4. I am also gioven a point P>in 3d, and I want to nd if>P is inside this tet. Right now, I just want something that will always>work, rather than the most efcient algorithm.The most straightforward way I know is to compute the barycentriccoordinates of P with respect to P1,...,P4. Let A be thematrix[P1 P2 P3 P4; 1 1 1 1]then the barycentric coordinates of P with respect to P1,...,P4 are given by inv(A)*[P; 1] (using matlab notation). If the barycentric coordinates are all positive, then P is inside the tetrahedron. (If youre further interested... If the ith component of the barycentricand if it is negative, then it is on the opposite side. In fact, the ithcomponent is a constant times the signed perpendicular distance from P tocheers,Rick =Right. To generalize slightly to the case where the points P1...Pn are not necessarily independent. Set up a linear program and solve for feasibility: the tableau is the same as below.>>I have a tet. in 3d given by 4 nodes P1-P4. I am also gioven a point P>>in 3d, and I want to nd if>>P is inside this tet. Right now, I just want something that will always>>work, rather than the most efcient algorithm.> The most straightforward way I know is to compute the barycentric> coordinates of P with respect to P1,...,P4. Let A be the> matrix> [P1 P2 P3 P4;> 1 1 1 1]> then the barycentric coordinates of P with respect to P1,...,P4 are given > by inv(A)*[P; 1] (using matlab notation). If the barycentric coordinates > are all positive, then P is inside the tetrahedron. > (If youre further interested... If the ith component of the barycentric> and if it is negative, then it is on the opposite side. In fact, the ith> component is a constant times the signed perpendicular distance from P to> cheers,> Rick> -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. ---Randomness comes in bunches. =>I have a tet. in 3d given by 4 nodes P1-P4. I am also gioven a point P>in 3d, and I want to nd if>P is inside this tet. Right now, I just want something that will always>work, rather than the most efcient algorithm.> The most straightforward way I know is to compute the barycentric> coordinates of P with respect to P1,...,P4. Let A be the> matrix> [P1 P2 P3 P4;> 1 1 1 1]> then the barycentric coordinates of P with respect to P1,...,P4 are given > by inv(A)*[P; 1] (using matlab notation). If the barycentric coordinates > are all positive, then P is inside the tetrahedron. > (If youre further interested... If the ith component of the barycentric> and if it is negative, then it is on the opposite side. In fact, the ith> component is a constant times the signed perpendicular distance from P to> cheers,> RickCorrect, but he also needs to check that det(A)>0, ie., hispoint numbering must insure a positive tet volume. If the volume isnegative the coordinate condition reverses. =>>>I have a tet. in 3d given by 4 nodes P1-P4. I am also gioven a point P>>in 3d, and I want to nd if>>P is inside this tet. Right now, I just want something that will always>work, rather than the \ most efcient algorithm.>>> The most straightforward way I know is to compute the barycentric>> coordinates of P with respect to P1,...,P4. Let A be the>> matrix>> [P1 P2 P3 P4;>> 1 1 1 1]>> then the barycentric coordinates of P with respect to P1,...,P4 are given >> by inv(A)*[P; 1] (using matlab notation). If the barycentric coordinates >> are all positive, then P is inside the tetrahedron. >>> (If youre further interested... If the ith component of the barycentric>> and if it is negative, then it is on the opposite side. In fact, the ith>> component is a constant times the signed perpendicular distance from P to>>> cheers,>> Rick>>Correct, but he also needs to check that det(A)>0, ie., his>point numbering must insure a positive tet volume. If the volume is>negative the coordinate condition reverses.No, it is unnecessary to check the orientation of the tetrahedron (more generally, simplex), and the coordinate condition remains the same. For example, Let A1=[P1 P2 P3 P4; 1 1 1 1] and let [b1 b2 b3 b4] = inv(A)*[x; 1], then for A2 = [P2 P1 P3 P4; 1 1 1 1], we have inv(A2)*([x; 1] = [b2 b1 b3 b4].By denition, the barycentric coordinates [b1 ... b4] (for d=3) of apoint x with respect to P1,...P4 isx = sum_{i=1,...,4} bi*Pi and sum_{i=1,...,4} bi = 1. It would be disturbing if computing the barycentric coordinates dependedon a specic orientation of P1,...,P4.(The above assumes, as another poster pointed out, that P1,...,P4 areafnely independent, i.e. are not all in the same plane)cheers, Rick =If one considers that life is based on Carbon 6, then one can also assumethat the hidden messages that determine human destiny can also be exploredby using some kind of decoding method using 6.If one assume that all letters A to Z of English language are based onCarbon 6 code, so then for example A=6, B=12, C=18 etc and Z=156.If any words that have historical meaning, and their numbers add up to anyof the three fundamental 3 digit number 222, 444, 666, this may help todecipher some hidden messages.For example: 666= Computer = 18+ 90+ 78+ 96+ 126+120+30+108 ( the modern civilizationbased on microchip? the ultimate limit of human knowledge?)similarly 444= Cross, English, Jesus, Lucifer, Messiah,and 222= India, ( some hidden energy source? Enlightenment? )and also note the sum of all numbers for letters A to Z is = 2025 (Spooky?)So, can anybody help us decipher any more words, which may guide Humancivilization towards the Holy Grail?R2 D2 & 4 Horsemen =Does it work for Latin too? Or Yiddish?Numerology is not mathematics. Its not even a (exact) science. You might want to move this discussion to sci.numerology to wherefollow-ups have been set.| If one considers that life is based on Carbon 6, then one can also assume| that the hidden messages that determine human destiny can also beexplored| by using some kind of decoding method using 6.|| If one assume that all letters A to Z of English language are based on| Carbon 6 code, so then for example A=6, B=12, C=18 etc and Z=156.| If any words that have historical meaning, and their numbers add up to any| of the three fundamental 3 digit number 222, 444, 666, this may help to| decipher some hidden messages.|| For example:|| 666= Computer = 18+ 90+ 78+ 96+ 126+120+30+108 ( the modern civilization| based on microchip? the ultimate limit of human knowledge?)|| similarly 444= Cross, English, Jesus, Lucifer, Messiah,|| and 222= India, ( some hidden energy source? Enlightenment? )|| and also note the sum of all numbers for letters A to Z is = 2025(Spooky?)||| So, can anybody help us decipher any more words, which may guide Human| civilization towards the Holy Grail?|| R2 D2 & 4 Horsemen|| =I believe this is an LP problem, but am not yet familiar enough with LPconcepts to be certain. If anyone cares to assist, here goes:I have an avocado farm and want to maximize my proceeds each harvest period.I have three supermarket chains that bid on my fruit, and they will give mea price for each avocado.In order to keep my business going, I have agreed to the following:I will sell at least 80% of my harvest (by weight) to Chain A.I will sell at least 10% of my harvest to Chain B.Additionally, I have agreed that at least 20% of the avocados sold to Chain A will be 3oz. or less.I have agreed that no more than 10% of the avocados sold to Chain A will bemore than 6 oz.I have 1,000 avocados with a weight and a price from each of the threechains and I want to know where to send each avocado to maximize my salesand satisfy my agreements.So, is this an LP problem? And, if so, would anyone care to recommend agood resource to learn how to program such a problem in C? = >I believe this is an LP problem, but am not yet familiar enough with LP >concepts to be certain. If anyone cares to assist, here goes: >I have an avocado farm and want to maximize my proceeds each harvest period. >I have three supermarket chains that bid on my fruit, and they will give me >a price for each avocado. >In order to keep my business going, I have agreed to the following: >I will sell at least 80% of my harvest (by weight) to Chain A. >I will sell at least 10% of my harvest to Chain B. >Additionally, >I have agreed that at least 20% of the avocados sold to Chain A will be 3 >oz. or less. >I have agreed that no more than 10% of the avocados sold to Chain A will be >more than 6 oz. >I have 1,000 avocados with a weight and a price from each of the three >chains and I want to know where to send each avocado to maximize my sales >and satisfy my agreements. >So, is this an LP problem? And, if so, would anyone care to recommend a >good resource to learn how to program such a problem in C? > isnt it an integer LP? you will hardly sell 0.27 avocados.avocado i selled to chain j is in 0 1 (true false) these are thevariables. and from this you can, knowing weight of each avocado and the prices the chains pay you can indeedsolve it this way. but in truth: isnt it a stochastic lp? will you indeed know (you make your contracts beforehand i guess) the weight of each avocado?hope this helpspeter =It is a MIP (Mixed Integer Problem).You can use lp_solve to solve this. Seehttp://groups.yahoo.com/group/lp_solve/It is a free library and can be called from C.Peter> I believe this is an LP problem, but am not yet familiar enough with LP> concepts to be certain. If anyone cares to assist, here goes:> I have an avocado farm and want to maximize my proceeds each harvestperiod.> I have three supermarket chains that bid on my fruit, and they will giveme> a price for each avocado.>> In order to keep my business going, I have agreed to the following:>> I will sell at least 80% of my harvest (by weight) to Chain A.> I will sell at least 10% of my harvest to Chain B.>> Additionally,>> I have agreed that at least 20% of the avocados sold to Chain A will be 3> oz. or less.> I have agreed that no more than 10% of the avocados sold to Chain A willbe> more than 6 oz.>> I have 1,000 avocados with a weight and a price from each of the three> chains and I want to know where to send each avocado to maximize my sales> and satisfy my agreements.> So, is this an LP problem? And, if so, would anyone care to recommend a> good resource to learn how to program such a problem in C?>> => I believe this is an LP problem, but am not yet familiar enough with LP> concepts to be certain. If anyone cares to assist, here goes:> I have an avocado farm and want to maximize my proceeds each harvest period.> I have three supermarket chains that bid on my fruit, and they will give me> a price for each avocado.> In order to keep my business going, I have agreed to the following:> I will sell at least 80% of my harvest (by weight) to Chain A.> I will sell at least 10% of my harvest to Chain B.> Additionally,> I have agreed that at least 20% of the avocados sold to Chain A will be 3> oz. or less.> I have agreed that no more than 10% of the avocados sold to Chain A will be> more than 6 oz.> I have 1,000 avocados with a weight and a price from each of the three> chains and I want to know where to send each avocado to maximize my sales> and satisfy my agreements.> So, is this an LP problem? And, if so, would anyone care to recommend a> good resource to learn how to program such a problem in C?> I think that the problem needs to be more clearly stated. Are there onlysize restrictions for Chain A? By promising to sell at least 80% of the harvest to Chain A, the problemwith size restrictions could be infeasible. Or is it that at least 80%(by weight) of the avocados actually sold should go to Chain A?What are the prices paid by the Chains? Are they sold by weight? In sizebands? There seem to be three major size bands; do the prices depend onthe size bands?The problem could be formulated as a standard LP, but you might getsolutions involving fractions, but this probably wouldnt matter in apractical situation.Johanneshttp://sizetter.com =I know that Pochhammer symbol is (x)_n and it is equal to x*(x+1)*(x+2)..(x+n-1).When all xs are different, I have x_0*(x_1+1)*(x_2+2)...(x_n+n) for (n+1) factors or x_0*(x_1+1)*(x_2+2)...(x_(n-1)+(n-1)) for n factors.If you know the expansion of x_0*(x_1+1)*(x_2+2)...(x_n+n), please let me know.References of this expansion are appreciated.John =Ornarose Inc. has an opening for a Research SoftwareDeveloper in Northern New Jersey or Chicago, IL. This is ashort-term (5 to 6 month) position with an SBIR-supportedstartup company. Developer will be responsible forimplementation and testing of algorithms for supervised andunsupervised learning, prediction, classication, etc.Work includes data set cleanup, experimentation, andmeasuring resource demands of algorithms. Requirements: * B.S. or advanced degree in computer science,statistics, or related eld. 5+ years professionalsoftware development experience. 2+ years professionalexperience with one or more of the following: machinelearning, data mining, statistics, pattern recognition,numerical optimization, numerical analysis. Experience withtext processing is also desirable, as is experience withdesigning and running computational experiments. * C and Perl prociency. C++ prociency is desirable. * Excellent verbal and written communication skills inEnglish. Demonstrated ability to meet deadlines andcommunicate effectively when working from home.Interested parties should send a resume (leads also welcome)in ASCII or opportunity employer. Because the position beginsimmediately, the candidate must be eligible to work legallyin the United =as if it is commonly known. Unfortunatly I have not got the faintest clue what it means, and have trouble nding texts on the subject. I only need a supercial explanation and maybe some references to available texts.Please advice.Christopher Grinde => as if it is commonly known.Its pretty simple. A Krylov sequence is the sequence r_{i+1}=Ar_i.Most iterative methods for linear systems are based on this ideasomehow, and you can also use it for eigenvalue problems.The subspace simply refers to the fact that you only form a nitepart of this innite series. You then try solving a linear system bytaking r_1=Ax_1-b, and setting x_{i+1} = x_1 + [something in the span ofthe rst i Krylov vectors].V.-- => as if it is commonly known. Unfortunatly I have not got the faintest > clue what it means, and have trouble nding texts on the subject. I > only need a supercial explanation and maybe some references to > available texts.> Please advice.http://www.maths.lth.se/na/courses/EDA410/EDA410-01/ manus/manus9.pdfBob Kolker =Christopher Grinde> as if it is commonly known. Unfortunatly I have not got the faintest >> clue what it means, and have trouble nding texts on the subject. I >> only need a supercial explanation and maybe some references to >> available texts.>> Please advice.>> http://www.maths.lth.se/na/courses/EDA410/EDA410-01/manus/ manus9.pdf> Bob Kolker>>-- =Your explanation was very if it is commonly known.>> Its pretty simple. A Krylov sequence is the sequence r_{i+1}=Ar_i.>> Most iterative methods for linear systems are based on this idea> somehow, and you can also use it for eigenvalue problems.>> The subspace simply refers to the fact that you only form a nite> part of this innite series. You then try solving a linear system by> taking r_1=Ax_1-b, and setting x_{i+1} = x_1 + [something in the span of> the rst i Krylov vectors].>> V.>-- =Im trying to use call SGESV from a C program (well, actually Squeak calling a C library, but anyway). It seems to work, except that the resulting contents of the ipiv vector dont make sense to me.Starting with an A and B matrix:[-4.0 -1.0 -10.00.0 9.0 1.02.0 5.0 7.0]respectively, everything works out ok.However, starting with an A matrix of:[2.0 5.0 7.00.0 9.0 1.0-4.0 -1.0 -10.0]and the correspondingly permuted B matrix, the resultingvalues in ipiv are {3, 2, 3}. I was under the impressionthat each index must appear once and only once. I searchedonline in vain to nd documentation of the (presumably)bijective mapping between LAPACK permutation vectors andpermutation matrices.Am I interpreting the ipiv vector incorrectly? I shouldalso mention that I am calling the function as follows: sgesv_(&n, &nrhs, a, &n, ipiv, b, &n, &info);and that the value stored into info is 0.Joshua => However, starting with an A matrix of:> [2.0 5.0 7.0> 0.0 9.0 1.0> -4.0 -1.0 -10.0]> and the correspondingly permuted B matrix, the resulting> values in ipiv are {3, 2, 3}. I was under the impression> that each index must appear once and only once. I searched> online in vain to nd documentation of the (presumably)> bijective mapping between LAPACK permutation vectors and> permutation matrices.> Am I interpreting the ipiv vector incorrectly? In order, the three values of ipiv mean that the rst equation wasinterchanged with the third equation (the rst 3), then the secondequation was interchanged with the second equation (the 2) (no actualinterchange), and nally, the third equation was interchanged withthe third equation (the second 3) (again, no actual interchange).Dave =>>However, starting with an A matrix of:>>[2.0 5.0 7.0>>0.0 9.0 1.0>>-4.0 -1.0 -10.0]>>and the correspondingly permuted B matrix, the resulting>>values in ipiv are {3, 2, 3}. I was under the impression>>that each index must appear once and only once. I searched>>online in vain to nd documentation of the (presumably)>>bijective mapping between LAPACK permutation vectors and>>permutation matrices.>>Am I interpreting the ipiv vector incorrectly? > In order, the three values of ipiv mean that the rst equation was> interchanged with the third equation (the rst 3), then the second> equation was interchanged with the second equation (the 2) (no actual> interchange), and nally, the third equation was interchanged with> the third equation (the second 3) (again, no actual interchange).> DaveAh, I thought it might be that, but couldnt satisfy myself that it wasDo you know where I might nd this documented?Joshua => In order, the three values of ipiv mean that the rst equation was> interchanged with the third equation (the rst 3), then the second> equation was interchanged with the second equation (the 2) (no actual> interchange), and nally, the third equation was interchanged with> the third equation (the second 3) (again, no actual interchange).> Dave> Ah, I thought it might be that, but couldnt satisfy myself that it was> Do you know where I might nd this documented?The comments in the code describe it, but rather succinctly.* IPIV (output) INTEGER array, dimension (N)* The pivot indices that dene the permutation matrix P;* row i of the matrix was interchanged with row IPIV(i).DAve =Could anybody set me on the right track to nd a routine for numericallyperforming a Laplace transform? I use Maple, Matlab or C (++).Dennis =Could someone please explain to me the geometric meaning of a column space? I just dont see it. In class we derived it from the vector space of Ax = b where A is an (m by n) matrix and x and b are vectors of size m and n respectively. And to be honest Im not sure whether I understand the geometric meaning behind this vector space either.Anybody ?ThanxDoug => Could someone please explain to me the geometric meaning of a column > space? I just dont see it. In class we derived it from the vector > space of Ax = b where A is an (m by n) matrix and x and b are vectors of > size m and n respectively. And to be honest Im not sure whether I > understand the geometric meaning behind this vector space either.> Anybody ?> Thanx> Doug> Perhaps to help me clear this up, someone could answer me this:In a system Ax=b where A is a (1 by n) matrix, and x a vector of length n and b merely one value; for any imaginable linear system A of this sort, we can come up with values for all x such that any value for b can be achieved. Does this always hold?ThanxDoug =>> Could someone please explain to me the geometric meaning of a column >> space? I just dont see it. In class we derived it from the vector >> space of Ax = b where A is an (m by n) matrix and x and b are vectors of >> size m and n respectively. And to be honest Im not sure whether I typo: n and m>> understand the geometric meaning behind this vector space either.>>> Anybody ?>> Thanx>> Doug>>Perhaps to help me clear this up, someone could answer me this:>In a system Ax=b where A is a (1 by n) matrix, and x a vector of length >n and b merely one value; for any imaginable linear system A of this >sort, we can come up with values for all x such that any value for b can >be achieved. Does this always hold?Doug,Ill try to guess what could have been bothering you ...The column space of the matrix A could be better understood as therange space or image of the linear transformation represented bythe matrix A under the canonical basis, i.e. [1,0,...], [0,1,0,...],[0,0,1,0,...], ...Why do we talk about column space? This can be said to be due tothe convention of writing the matrix-vector-multiplication as(1) A x = b (a matrix A right-multiplied by a column vector x).In this case, we say that the vector x in the n-dim. space issent by the matrix A to the m-dimensional vector b and b is a linear combination of the columns of A.(We are using the equality sign here seriously now.)Sometimes, people would also write an equation as(2) y B = c (a matrix B left-multiplied by a row-vector y)yielding an n-dimensional row-vector c. (Please check the dimensionality).In this case, we are considering a linear transformation sending m-dimensional vectors to n-dimensional vectors and the row-vector c is a linear combination of the row-vectors in the matrix B with the coeffcients the coordinates of the vector y. (Why coordinates? Please recall that we are using the canonical basis whenever we write a linear transformation as a matrix with numbers. Also, please convince yourself of these ideas of a linear combination on a piece of paper --- partition the matrix and the vector, so that we could see that directly as a linear combination. Please also note that there are at least two ways of partitioning the matrix and the vector for the multiplication. We choose the proper one we need. This can also be found in some linear algebra texts.)The default range space or image of the linear transformationin question (or implied) is now the row space.So, we have column space and row space. How are they connected?In (1) and (2), weve used different matrices and vectors to avoidconfusion. We consider now the matrix A only and write it as(3) A x = bto which I myself am used. *By convention*, the range space or imageof A is the column space. But the row space of A exists and is alsoa very important fundamental subspace (Gilbert Strang) associated tothe matrix A. We can rewrite (3) (A x) = x A = band talk about the row-vectors x, b and A (the transpose of A) ANDthe row space as the range space since we are now sending therow-vector x via A to b. *To avoid confusion*, wed better talkabout the column space or the raw space rather than the rangesapce.These two subspaces of A are well-connected, but Im still consideringif I should try that immediately ..... A full geometrical interpretationof these concepts is well-exposed by the so-calledsingular value decomposition (SVD) both theoretically and practically.There are two further important subspaces which are called thenullspace (kernel) and the left nullspace (left kernel). Byusing the word nullspace, we usually mean the convention ofwriting A x = 0i.e. the column-convention. On the other hand, the solution set of y B = 0is called the left nullspace (matrix B *LEFT*-multiplied by a row-vectory), because the solution set can be shown to be a subspace. So, the usualname nullspace may also be thought of as theRIGHT-nullspace with RIGHT implied. We have mentioned so farfour fundamental subspace. These 4 subspaces are always associated with every single matrix a priori and are very tightly connected. The complexof the relations between these subspaces are IMHO best exposed by the SVDby drawing a picture (see e.g. that one by Gilbert Strang. I have made anothermyself, but I think there would be too many symbols in it for you now). Would somebody provide a link to the diagram by Strang for Doug?So, I have not even started with possible geometrical interpretations ....and I think considering a (1 by n) matrix A would be a good approachto see something geometrically, but I hesitate in this moment towrite more. If you respond, Ill surely try to reply (please Cc a copy to remind me.)Rudi---- =I am trying to solve an linear integer minimization problem with linearconstraints involving about 150 variables. Is there free software I can usefor a problem of that size?John =Yes there are. Take a look at lp_solve. You can nt it athttp://groups.yahoo.com/group/lp_solve/It is completely free.Peter> I am trying to solve an linear integer minimization problem with linear> constraints involving about 150 variables. Is there free software I canuse> for a problem of that size?>> John>> = >I am trying to solve an linear integer minimization problem with linear >constraints involving about 150 variables. Is there free software I can use >for a problem of that size? >John >yes, lpsolve from will do thatftp://ftp.ics.ele.tue.nl/pub/projects/lp_solvehope this helpspeter =I have experimental data (x_i, y_i) where the y_i have Gaussian noise with a mean of 0 and known standard deviation, s. Since the x_i are irregularly spaced, I am interpolating the y_i onto a grid of equally-spaced x. This is being done by linear interpolation; i e, y_interpolated = y_0 + (y_1 - y_0) * (x_interpolated-x_0)/(x_1 - x_0), using the points (x_0,y_0) and (x_1,y_1) which bracket each x_interpolated on the equally-spaced grid.My question is, what is the standard deviation of the resulting y_interpolated?A simple-minded analysis suggests it is unchanged at s, since y_interpolated = y_0 * f + y_1 *(1-f), and the variances of y_0 and y_1 are s^2, so the variance of y_interpolated is s^2*f + s^2*(1-f) = s^2.Numerical experiments do not bear this out, however; calculating the standard deviation of interpolated data gives a value about 0.8s. What is wrong with the above analysis? => I have experimental data (x_i, y_i) where the y_i have Gaussian noise> with a mean of 0 and known standard deviation, s. Since the x_i are> irregularly spaced, I am interpolating the y_i onto a grid of> equally-spaced x. This is being done by linear interpolation; i e,> y_interpolated = y_0 + (y_1 - y_0) * (x_interpolated-x_0)/(x_1 - x_0),> using the points (x_0,y_0) and (x_1,y_1) which bracket each> x_interpolated on the equally-spaced grid.> My question is, what is the standard deviation of the resulting> y_interpolated?> A simple-minded analysis suggests it is unchanged at s, since> y_interpolated = y_0 * f + y_1 *(1-f), and the variances of y_0 and y_1> are s^2, so the variance of y_interpolated is s^2*f + s^2*(1-f) = s^2. Var[y_interp] = s_0^2 * f^2 + s_1^2 * (1-f)^2> Numerical experiments do not bear this out, however; calculating the> standard deviation of interpolated data gives a value about 0.8s. Not surprising.> What is wrong with the above analysis? 0.707 = 1/sqrt(2) <= sqrt [ f^2 + (1-f)^2 ] <= 1-- Julian V. NobleProfessor Emeritus of ^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt Brecht, Galileo. =I have Macsyma 2.4. It is a very powerful symbolic tool, and it comes withmany interesting demos and examples. However, there are almost 0 freelyavailable resources, archives, or discussions on the Internet.Any ideas ?Andrew Isaverdian => I have Macsyma 2.4. It is a very powerful symbolic tool, and it comes with> many interesting demos and examples. However, there are almost 0 freely> available resources, archives, or discussions on the Internet.> Any ideas ?>> Andrew IsaverdianMaxima is a close relative of Macsyma which is freely available. To getstarted on third party Maxima/Macsyma code, try:http://maxima.sourceforge.net/3rdpartycode.shtmlMike Thomas. SLATEC functions (available from NETLIB)to my calculations on IBM PC, using C/C++. However, I am not surehow to choose the machine constants suitable for C and forPentium MMX processor that I use. The SLATEC function D1MACH(I) providesthe constants for FORTRAN double precision variables, but at someadditional assumptions that I suppose may not hold for the aboveconditions (I am including relevant fragments of this function below).I need machine constants for oat, double, and long double C variables,preferably for the whole Pentium family. Does anyone knowL.B. DOUBLE PRECISION FUNCTION D1MACH (I)CC D1MACH can be used to obtain machine-dependent parameters for theC local machine environment. It is a function subprogram with oneC (input) argument, and can be referenced as follows:CC D = D1MACH(I)CC where I=1,...,5. The (output) value of D above is determined byC the (input) value of I. The results for various values of I areC discussed below.CC D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude.C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude.C D1MACH( 3) = B**(-T), the smallest relative spacing.C D1MACH( 4) = B**(1-T), the largest relative spacing.C D1MACH( 5) = LOG10(B)CC Assume double precision numbers are represented in the T-digit,C base-B formCC sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )CC where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), andC EMIN .LE. E .LE. EMAX.CC The values of B, T, EMIN and EMAX are provided in I1MACH asC follows:C I1MACH(10) = B, the base.C I1MACH(14) = T, the number of base-B digits.C I1MACH(15) = EMIN, the smallest exponent E.C I1MACH(16) = EMAX, the largest exponent E.CCC MACHINE CONSTANTS FOR THE IBM PCC ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISIONC ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087.CC DATA SMALL(1) / 2.23D-308 /C DATA LARGE(1) / 1.79D+308 /C DATA RIGHT(1) / 1.11D-16 /C DATA DIVER(1) / 2.22D-16 /C DATA LOG10(1) / 0.301029995663981195D0 /*----------------------------------------------------------- --------*| Dr. Leslaw Bieniasz, || Institute of Physical Chemistry of the Polish Academy of Sciences,|| Molten Salts Department, ul. Zagrody 13, 30-318 Cracow, Poland. || tel./fax: +48 (12) 266-03-41 |*----------------------------------------------------------- --------*| Interested in Computational Electrochemistry? || Visit my web site: http://www.cyf-kr.edu.pl/~nbbienia |*----------------------------------------------------------- --------* some SLATEC functions (available from NETLIB) >to my calculations on IBM PC, using C/C++. However, I am not sure >how to choose the machine constants suitable for C and for >Pentium MMX processor that I use. The SLATEC function D1MACH(I) provides >the constants for FORTRAN double precision variables, but at some >additional assumptions that I suppose may not hold for the above >conditions (I am including relevant fragments of this function below). >I need machine constants for oat, double, and long double C variables, >preferably for the whole Pentium family. Does anyone knowfor double precision the constants in d1mach for the ibm pc are correct alsofor the pentium . for single precision, you nd the same in r1mach (also in slatec and in blas)netlib/port contains a routine s1mach for testing the consistency of theseconstants. for true long double you can compute them yourself by doubling /halng the increment for one (until 1+term=1 ) resp. until you obtain underlow or overow. you are sure that long double on your system is notidentical with double? (because of storage mode)?hope this helpspeter =You can nd a brief Acrobat le (ent15.pdf) with the IEEE machine constantsfor single and double precision at: I am trying to adapt some SLATEC functions (available from NETLIB)> to my calculations on IBM PC, using C/C++. However, I am not sure> how to choose the machine constants suitable for C and for> Pentium MMX processor that I use. The SLATEC function D1MACH(I) provides> the constants for FORTRAN double precision variables, but at some> additional assumptions that I suppose may not hold for the above> conditions (I am including relevant fragments of this function below).> I need machine constants for oat, double, and long double C variables,> preferably for the whole Pentium family. Does anyone know>> L.B.--There are two things you must never attempt to prove: the unprovable -- andthe obvious.--Democracy: The triumph of popularity over principle.--http://www.crbond.com NETLIB)> to my calculations on IBM PC, using C/C++. However, I am not sure> how to choose the machine constants suitable for C and for> Pentium MMX processor that I use. The SLATEC function D1MACH(I) provides> the constants for FORTRAN double precision variables, but at some> additional assumptions that I suppose may not hold for the above> conditions (I am including relevant fragments of this function below).> I need machine constants for oat, double, and long double C variables,> preferably for the whole Pentium family. Does anyone know> L.B.> You can also try going the C route directly. Something like:#include #include #include f2c.hdoublereald1mach_(integer *I){ switch (*I) { case 1: return DBL_MIN; case 2: return DBL_MAX; case 3: return DBL_EPSILON; case 4: return 2.*DBL_EPSILON; case 5: return M_LN2/M_LN10; default: return 2.*DBL_EPSILON; /* Should dump to console also. */ }} =hii have a number of questions, hopefully someone knows the answers:1. i nd programming mathematica to be more concise and rapidly developed than c++, but with a signicant speed hit (100x). is the similar to what others experience?2. im considering using matlab. what is the speed hit of matlab versus c++?3. what is the decrease in development time of matlab versus c++ (ierapid prototyping)?4. on c++, i use clapack, but i have the g77 compiler available. is clapack slower/faster than f77-based lapack?5. what is the highest level, easiest to use, numerical library available for c++? i am currently using tnt (which is dated), and ublas (which is not very useful) as => 1. i nd programming mathematica to be more concise and rapidly> developed than c++, but with a signicant speed hit (100x). is the> similar to what others experience?No I think mathematica is better for rapid small analysis an playing around.As long as there are no real OO routines inside C++ programming is moreconsise. And sadly mathematica is very slow, but it is not really anumerical tool.> 2. im considering using matlab. what is the speed hit of matlab versusc++?> Just one early comment. As you surely know, ATLAS architecture.> In a sense, it is similar to hardware-xed BLAS, but it adapts> to the hardware by itself. Matlab in its latest releases uses ATLAS> for performing basic arithmetic operations. That is why Matlab> (or several of its latest releases) is so fast!clue: using ATLAS (C) or MATLAB should achieve similar performance> 3. what is the decrease in development time of matlab versus c++ (ie> rapid prototyping)?Dont know that> 4. on c++, i use clapack, but i have the g77 compiler available. is> clapack slower/faster than f77-based lapack?Dont know either> 5. what is the highest level, easiest to use, numerical library> available for c++? i am currently using tnt (which is dated), and ublas> (which is not very useful) as simple wrappers around clapack.I like uBlas. Why do you think it is not very useful ?Patrick => 2. im considering using matlab. what is the speed hit of matlab versusc++?Well it depends what do you write in C++. E.g. convolution, ML is using ain a standard discrete way x[i]=sum(y(t-tau)*y(t)) ML is even faster. As youcan see its more a problem of the used algorithm. The computation is onething. The other the used parser for the ML scripts. Anyway, Mathworks doeshave an own compiler where you can use your own functions (s- and mexfunctions).> 3. what is the decrease in development time of matlab versus c++ (ie> rapid prototyping)?Developing time is faster by using ML. I have to write code in C++ sinceits used in other programs. Using ML I get results very quickly - thanporting to C++.> 5. what is the highest level, easiest to use, numerical library> available for c++? i am currently using tnt (which is dated), and ublas> (which is not very useful) as simple wrappers around clapack.IMO mtl but, this is the slowest lib too. Maybe these link may help:http://www.opencascade.org/upload/87/. If you need highest performance thereisnt any way arround fortran and their libs...Olaf =>> 3. what is the decrease in development time of matlab versus c++ (ie>> rapid prototyping)?To some extent it depends on your programming skills. Perhaps you are exceptionally fast C++ programmer and learning Matlab will be a waste of time for you. But most computational people are not like that. Matlab is a lot more comfortable, perfectly documented, with systematically organized libraries.It is perhaps illustrative to remind that MATLAB stands for MATrix LABoratory. It is a tool for developing algorithms, testing them, ne-tuning them and solving small to medium size problems. It is not a tool for number crunching applications (especially owing to simple memory management). On the other hand, I was surprised (almost shocked) the other day when I learnt that there are not adequate comprehensive, reliable, well-documented C++ libraries for LA, on which it would be possible to base a project easily. There are some nice projects but they still need some time to develop (after some search I have settled with uBLAS).>> 4. on c++, i use clapack, but i have the g77 compiler available. is>> clapack slower/faster than f77-based lapack?As far as I know, CLAPACK is just a C port of the Fortran version, the original LAPACK. >> 5. what is the highest level, easiest to use, numerical library>> available for c++? i am currently using tnt (which is dated), and ublas>> (which is not very useful) as simple wrappers around clapack.> I like uBlas. Why do you think it is not very useful ?In my opinion, uBLAS is a very promissing project. However, it provides the functionalities of BLAS only (basic arithmetic operations with vectors and matrices). Use of bindings to some other libraries (ATLAS, LAPACK) is thus inevitable, if you want to solve equations, perform SVD, LU, and so on. The idea of full OOP is thus not fullled. Not yet:-))) Perhaps one day we will have uLAPACK :-))))))))))Zdenek Hurak => Well it depends what do you write in C++. E.g. convolution, ML is using a> convolution in a standard discrete way x[i]=sum(y(t-tau)*y(t)) ML is even> faster. As you can see its more a problem of the used algorithm. The> computation is one thing. The other the used parser for the ML scripts.> Anyway, Mathworks does have an own compiler where you can use your own> functions (s- and mex functions).This is not correct. You claim that Matlab is fast for convolution because it uses FFTW. This product is free and written in C, so I can use it from my C/C++ code easily. I dont need matlab to achieve this great performance for convolution.If you write the convolution the standard way, Matlab used free ATLAS for manipulations with vectors and matrices like inner product, outer product and so on. I can easily write my own C/C++ code using ATLAS and I will achieve the same performance. >> 3. what is the decrease in development time of matlab versus c++ (ie>> rapid prototyping)?> Developing time is faster by using ML. I have to write code in C++ since> its used in other programs. Using ML I get results very quickly - than> porting to C++.Exactly. I believe that this is a very nice thing about Matlab. Developing and implementing algorithms is so easy and convenient.Zdenek => 2. im considering using matlab. what is the speed hit of matlab versus c++?On aggregate operations (i.e., vector and matrix operations, notexpressed as explicit for-loops), Matlab, Octave, R, and doubtlessother packages will pass the arguments through to optimized Fortranor assembler functions, so using the package is usually faster than a C++ program, often much faster.By the way, Octave provides the same basic functionality as Matlab,and it is free (as in free beer & free speech). The main differenceis that Matlab has a lot of add-on packages; if you dont need them,you might as well save yourself the expense of Matlab.> 3. what is the decrease in development time of matlab versus c++ (ie> rapid prototyping)?Using a package (Matlab, Octave, R, etc) that allows you to workdirectly with vectors and matrices is much, much faster for developmentthan C++ -- I would say at least 10 times faster.> 4. on c++, i use clapack, but i have the g77 compiler available. is > clapack slower/faster than f77-based lapack?I think your best strategy is to nd the most mature and, secondarily,most optimized version of LAPACK, and then gure out how to link it in.I dont know for sure, but I doubt if CLAPACK is most mature ormost optimized on any platform.> 5. what is the highest level, easiest to use, numerical library > available for c++? i am currently using tnt (which is dated), and ublas > (which is not very useful) as simple wrappers around clapack.Dont know for sure, but heres a place to look: http://math.haifa.ac.il/msoftware.htmlFor what its worth,Robert Dodier--``If I have not seen as far as others, it is because giants were standing on my shoulders. -- Hal Abelson =Just my opinion:I have never used MATLAB, but I know many swear by it (not at it :-) ). It has good plotting routines and is more rapid development than C++. You do take a speed hit in MATLAB, but you can link to C++ compiledprograms, but I understand that is tricky. I would suggest, however,that you look into MATLAB since, if you have the money for it there aremany packages of code ready to use. Its a pretty powerful workingenvironment.I use the combination Python and C++. Python is about the speed ofMATLAB. It is more rapid in development than even MATLAB is from whatI hear. It certainly is light years more rapid than C++. Its aclean, beautiful language that has many numerical packages and someplotting packages. I guarantee that once you use Python youll notwant to code in anything else. The whole coding cycle is marvelouslyclean and fast (no edit, compile, link, debug cycle -- only edit andrun -- you debug as you go with lots of feed back). Python can alsolink to C++ code. There are several packages out there that do thissemi-automatically and I am now looking into them. So you can get upand running very quickly in Python, then nd the bottlenecks and re-dothem later in C++.Downside: Python is not nearly as integrated as MATLAB from what Iveseen. Plotting is spotty, but there are some good packages. Nevertheless, youll not have many of the plotting niceties of MATLAB. You will have to put together your Python packages as they come frommany sources (although the main Python language package is in onelump).Upside (beside tremendous ease of use): Python is FREE as are many ofthe packages including the numeric packages. The user community is bigand supportive. Its an established language that is also continuingto grow. Many in science and math are nding Python good for RAD sonumeric support is good.You can see that Im prejudiced, but that also comes with 25+ years ofcoding in Fortran, C, C++, Mathematica, BASIC, APL (many years ago),machine language, and now Python. Its worth a look.-- Lou Pecora - My views are my own. => Well it depends what do you write in C++. E.g. convolution, ML is usinga> convolution in a standard discrete way x[i]=sum(y(t-tau)*y(t)) ML iseven> faster. As you can see its more a problem of the used algorithm. The> computation is one thing. The other the used parser for the ML scripts.> Anyway, Mathworks does have an own compiler where you can use your own> functions (s- and mex functions).>> This is not correct. You claim that Matlab is fast for convolution because> it uses FFTW. This product is free and written in C, so I can use it from> my C/C++ code easily. I dont need matlab to achieve this greatperformance> for convolution.>> If you write the convolution the standard way, Matlab used free ATLAS for> manipulations with vectors and matrices like inner product, outer product> and so on. I can easily write my own C/C++ code using ATLAS and I will> achieve the same performance.Well, thats right. The problem behind this is, that ML imo mostly does havethe best algorithm behind, because this is their/Mathworks job. If you arenot friend with Mathematics and FFT etc. you will have to pay the penalty.The same with kalman etc. or the use of banded/special matrizes - ML doesall this internally for you and, its quite fast. With other words, youspend more time on programming etc. as on using ML for circumstances.Further more, the decision depends on the goal, what do you want to do. As Imentioned before I have to use both.Olaf =I am trying to nd test suite for statistical testing ofrandom sequences (like DIEHARD Statistical Tests) in Matlab. Any hints?Michal => I am trying to nd test suite for statistical testing of> random sequences (like DIEHARD Statistical Tests) in Matlab. Any hints?> MichalStephen Wolfram has suggested a random number generatorbased on a particular Cellular automaton some years agoand he has accompanied the algorithm with quite a richcollection of tests that might be interesting for you.Sorry I have no precise references right now.Ahoj-- Pavel PokornyMath Dept, Prague Institute of Chemical Technologyhttp://staffold.vscht.cz/mat/Pavel.Pokorny =I have a nice problem, but I am not even sure if there exist a solution ;-)Given two matrices X and Y with X=A.t*A and Y=A*A.t, is there a possibilityto determine A?Patrick = I have a nice problem, but I am not even sure if there exist a solution ;-) >Given two matrices X and Y with X=A.t*A and Y=A*A.t, is there a possibility >to determine A? >yes. take the spectral decomposition of X and Y: X=V*diag(lambda_i)*V V unitary Y=U*diag(lambda_i,0,...0)*U U unitary (I assume Y is the larger one, hence if you know this relation its additioanleigenvalues will all be zero ) then A=U*S*V where S is matrix of the same shape as A, with its entries on the diagonal sqrt(lambda_i), and zero otherwise.look up svdhope this helps and wasnt homeworkpeter =>>I have a nice problem, but I am not even sure if there exist a solution ;-)>>Given two matrices X and Y with X=A.t*A and Y=A*A.t, is there a possibility>to determine A?>> yes. take the spectral decomposition of X and Y:> X=V*diag(lambda_i)*V V unitary> Y=U*diag(lambda_i,0,...0)*U U unitary> (I assume Y is the larger one, hence if you know this relation its additioanl> eigenvalues will all be zero ) then> A=U*S*V> where> S is matrix of the same shape as A, with its entries on the diagonal> sqrt(lambda_i), and zero otherwise.But note that A is not unique - since each square root can be taken with arbitrary sign, and if there are multiple eigenvalues then U and V involvefree essential choices, too. So one can determine _some_ A, but not _the_ A.Arnold Neumaier => yes. take the spectral decomposition of X and Y:> X=V*diag(lambda_i)*V V unitary> Y=U*diag(lambda_i,0,...0)*U U unitary> (I assume Y is the larger one, hence if you know this relation itsadditioanl> eigenvalues will all be zero ) then> A=U*S*V> where> S is matrix of the same shape as A, with its entries on the diagonal> sqrt(lambda_i), and zero otherwise.> look up> svdAs far as I could follow now, it is even better than I excpected while thesvd is the only real linear algebra function I am using in my algorithmsso far. I will dig in a bit deeper towmorrow....> hope this helps and wasnt homeworkYes it really helps, and for sure it was no homework (this time is over;- ) ) but if I regard my gaps I should do some additional homework.Patrick => But note that A is not unique - since each square root can be taken with> arbitrary sign, and if there are multiple eigenvalues then U and V involve> free essential choices, too.>> So one can determine _some_ A, but not _the_ A.This is what I expected, but in this context there is (I hope so!) noproblem as long as the Nullspace is still the same. This I will try out thenext days.Patrick = My 2 dimensional diffusion code is allowing the distributionfunction to go negative in regions. I am using rkc.f from netlib todo the time stepping and a library of higher order nite differenceapproximations which I have taken from W. E. Schiessers webpage athttp://www.lehigh.edu/~wes1/ctp/ctp.html for the spacial differencing. I should probably also note that my problem has variable diffusioncoefcients both on and off diagonal which are all about the samesize. In an earlier version of the same code, a colleague used ahopscotch method and an implicit method which I dont remember. Inall cases we get regions where the distribution is negative and whichgrow quickly. This, of course, is not physical. In the hopscotch andimplicit code the negative distributions could be avoided by settingthem to zero at each timestep, but we did not know at what cost, so wehave tried the method of lines and higher order differencingcoefcients. Again we get the same problems. Are there methods thatwill avoid this? I am wondering about symplectic methods, if Iunderstand correctly, they conserve phase space. I dont know thatthat precludes regions where the distribution goes negative, but doesanyone have experience with these? Is it worth trying a new method, orshould I just continue trying different approaches to this method. (For instance, I have been trying different boundary conditions). Isthere any more information I need to give to make my question clearer?Shawn = > My 2 dimensional diffusion code is allowing the distribution >function to go negative in regions. I am using rkc.f from netlib to >do the time stepping and a library of higher order nite difference >approximations which I have taken from W. E. Schiessers webpage at >http://www.lehigh.edu/~wes1/ctp/ctp.html for the spacial differencing. > I should probably also note that my problem has variable diffusion >coefcients both on and off diagonal which are all about the same >size. In an earlier version of the same code, a colleague used a >hopscotch method and an implicit method which I dont remember. In >all cases we get regions where the distribution is negative and which >grow quickly. This, of course, is not physical. In the hopscotch and >implicit code the negative distributions could be avoided by setting >them to zero at each timestep, but we did not know at what cost, so we >have tried the method of lines and higher order differencing >coefcients. Again we get the same problems. Are there methods that >will avoid this? I am wondering about symplectic methods, if I >understand correctly, they conserve phase space. I dont know that >that precludes regions where the distribution goes negative, but does >anyone have experience with these? Is it worth trying a new method, or >should I just continue trying different approaches to this method. >(For instance, I have been trying different boundary conditions). Is > there any more information I need to give to make my question clearer? >Shawnspurious oscillations can be introduced by:1) the space grid too coarse: if you write the differencing done as a star, it must satisfy the M-matrix property (diagonally weakly dominent L-structure)2) poor discretization of the boundary conditions. 3) effects from changing space-grids.4) effects of the time integrator. using higher order methods makes the occurence of spurious oscillations even moreprobable, hence this will be not the way to go. there are special methods which avoid this (osher-engquist type, ENO (essentially nonoscillating) ) but in the 2D case these are not easy to apply. since rkc is quite stable and order two, I would guess that most probably the effects come from the boundary conditions, which are the most troublesome part of the job.hope this helpspeter => My 2 dimensional diffusion code is allowing the distribution>function to go negative in regions. I am using rkc.f from netlib to>do the time stepping and a library of higher order nite difference>approximations which I have taken from W. E. Schiessers webpage at > spurious oscillations can be introduced by:> 1) the space grid too coarse: if you write the differencing done as > a star, it must satisfy the M-matrix property (diagonally weakly > dominent L-structure)> 2) poor discretization of the boundary conditions. > 3) effects from changing space-grids.> 4) effects of the time integrator. > using higher order methods makes the occurence of spurious oscillations even more> probable, hence this will be not the way to go.> there are special methods which avoid this (osher-engquist type,> ENO (essentially nonoscillating) ) but in the 2D case> these are not easy to apply. > since rkc is quite stable and order two, I would guess that most> probably the effects come from the boundary conditions, which are the> most troublesome part of the job.> hope this helps> peter I have suspected the boundaries as the problem, and have had somesuccess by simplifying them, making them straight lines, and evenusing simpler initial distributions. I would have thought that thesimple condition f=0 at the boundary would be benign, but apparentlynot. I still have not been able to get the spurious oscillations out. Unfortunately, I do not yet know what you mean in (1) in your list. What would be a good source to read more on this? = >> spurious oscillations can be introduced by: >> 1) the space grid too coarse: if you write the differencing done as >> a star, it must satisfy the M-matrix property (diagonally weakly >> dominent L-structure) >> 2) poor discretization of the boundary conditions. >> 3) effects from changing space-grids. >> 4) effects of the time integrator. > I have suspected the boundaries as the problem, and have had some >success by simplifying them, making them straight lines, and even >using simpler initial distributions. I would have thought that the >simple condition f=0 at the boundary would be benign, but apparently >not. I still have not been able to get the spurious oscillations out. >Unfortunately, I do not yet know what you mean in (1) in your list. >What would be a good source to read more on this? you have a twodimensional parabolic problem , variable coefcients,a mixed derivative sayu_t = a* u_xx + b* u_xy + c* u_yy + d *u_x + e*u_y +f(u,x,y).rst of all, you should approximate u_xy by the star (I assume delta_x=delta_y ; otherwise it becomes more complicated) (1/(2*h^2)) [ -1 1 0 ; 1 -2 1 ; 0 1 -1 ]; (in matlab matrix notation)if b<0 resp. (1/(2*h^2)) [ 0 -1 1 ; -1 2 -1 ; 1 -1 0 ]; if b>0. I assume that a>|b| and c>|b| next, the stepsize h should be such that a>|b| +(h/2)*|d| and c>|b|+(h/2)*|e|.if d or e are large, this introduces sharp restricitions on the stepsize. if this is the case, then you must use noncentral discretization, those ofengquist-osher type. otherwise spurious oscillations will always occur.hope this helpspeter =>> spurious oscillations can be introduced by:>> 1) the space grid too coarse: if you write the differencing done as >> a star, it must satisfy the M-matrix property (diagonally weakly >> dominent L-structure)>> 2) poor discretization of the boundary conditions. >> 3) effects from changing space-grids.>> 4) effects of the time integrator.> I have suspected the boundaries as the problem, and have had some>success by simplifying them, making them straight lines, and even>using simpler initial distributions. I would have thought that the>simple condition f=0 at the boundary would be benign, but apparently>not. I still have not been able to get the spurious oscillations out. >Unfortunately, I do not yet know what you mean in (1) in your list. >What would be a good source to read more on this?> you have a twodimensional parabolic problem , variable coefcients,> a mixed derivative > say> u_t = a* u_xx + b* u_xy + c* u_yy + d *u_x + e*u_y +f(u,x,y).> rst of all, you should approximate u_xy by the star > (I assume delta_x=delta_y ; otherwise it becomes more complicated) > (1/(2*h^2)) [ -1 1 0 ; 1 -2 1 ; 0 1 -1 ]; (in matlab matrix notation)> if b<0 resp.> (1/(2*h^2)) [ 0 -1 1 ; -1 2 -1 ; 1 -1 0 ]; > if b>0. > I assume that a>|b| and c>|b| > next, the stepsize h should be such that > a>|b| +(h/2)*|d| and c>|b|+(h/2)*|e|.> if d or e are large, this introduces sharp restricitions on the stepsize. if > this is the case, then you must use noncentral discretization, those of> engquist-osher type. otherwise spurious oscillations will always occur.> hope this helps> peterThis does help. I have a couple further questions, probably the mostimportant is for a recommendation of a text that would show me how toderive the above star for unequal step sizes, one that would alsoexplain how this damps the spurious oscillations. I am also a bitcurious about these engquist-osher type discretizations, I haventlooked at the d and e yet, but I am a little concerned that theymight be too large. I looked online and found references to thepaper ENGQUIST, B. & OSHER, S. (1981) One-sided differenceapproximations for nonlinear conservation laws. Math. Comp., 36,321[CapitalEth]35. It appears from this that you are suggesting I use anapproximation for a nonlinear system. Are equations with d or etoo large considered non-linear, or is the method just good fordifcult problems? (If this was also covered in the recommended textthat might be helpful). I should ask if my assumption that yourstatement d or e being too large means that the inequalities youcite are not able to be met. I am also assuming that the otherderivatives (including nite differencing of coefcients a, band c to get d and e) can be done with simple centereddifference approximations?Currently I only have access to some very old texts on PDEs, at leastas long as I stick with nite differencing. None of them mention theM-matrix property or what an L-structure is. I looked at your webpage hoping to nd the information there. It looks as if you havemade a lot of very useful tools and information public! UnfortunatelyI was not able to read the papers in the Teaching section as I donot know German.Shawn = >> say >>> u_t = a* u_xx + b* u_xy + c* u_yy + d *u_x + e*u_y +f(u,x,y). >> rst of all, you should approximate u_xy by the star >> (I assume delta_x=delta_y ; otherwise it becomes more complicated) >> (1/(2*h^2)) [ -1 1 0 ; 1 -2 1 ; 0 1 -1 ]; (in matlab matrix notation) >> if b<0 resp. >> (1/(2*h^2)) [ 0 -1 1 ; -1 2 -1 ; 1 -1 0 ]; >> if b>0. >>> I assume that a>|b| and c>|b| >>> next, the stepsize h should be such that >> a>|b| +(h/2)*|d| and c>|b|+(h/2)*|e|. >> if d or e are large, this introduces sharp restricitions on the stepsize. if think about a=c=2, b=1, d=-4000 e=8000 this is a convection dominated problem, and better handled like a hyperbolicequation. (there several ways to do this, for example there are splittingsteps possible: rst take only the hyperbolic part, then diffuse (forone time step) >This does help. I have a couple further questions, probably the most >important is for a recommendation of a text that would show me how to >derive the above star for unequal step sizes, one that would also >explain how this damps the spurious oscillations. I am also a bit >curious about these engquist-osher type discretizations, I havent >looked at the d and e yet, but I am a little concerned that they >might be too large. I looked online and found references to the >paper ENGQUIST, B. & OSHER, S. (1981) One-sided difference >approximations for nonlinear conservation laws. Math. Comp., 36, >321[CapitalEth]35. It appears from this that you are suggesting I use an >approximation for a nonlinear system. Are equations with d or e >too large considered non-linear, or is the method just good for >difcult problems? (If this was also covered in the recommended textthe nonlinearities cause the same effect as a strong convection >that might be helpful). I should ask if my assumption that your >statement d or e being too large means that the inequalities you >cite are not able to be met. I am also assuming that the other >derivatives (including nite differencing of coefcients a, b >and c to get d and e) can be done with simple centered >difference approximations? not necessarily, might be better by upwind schemes, decreasing the order >Currently I only have access to some very old texts on PDEs, at least >as long as I stick with nite differencing. None of them mention the >M-matrix property or what an L-structure is. I looked at your web >page hoping to nd the information there. It looks as if you have >made a lot of very useful tools and information public! Unfortunately >I was not able to read the papers in the Teaching section as I do >not know German. >L -structure is : diagonal strictly positive, nondiagonal elements nonpositive (<=0)M-matrix is a L-matrix for which the inverse exists and is componentwise positive.you can read M here as monotonic:If Ax = u, Ay = v, (u, v are the loads) and v>=u componentwiseand A is a M-matrix, then y >= x componentwise as it should . x,y deformation.the M-matrix property also gives the necessary stability property for thediscretization which is required for convergence as the gridsize goes to zero.a good text for your problem might beMikhail Shashkov: conservative nite difference methods on general gridscrc press , 1996, ISBN 0-8493-7375-1 and Thomas, J.W.Numerical partial differential equations: Finite difference methods. (English)Texts in Applied Mathematics. 22. New York, NY: Springer-Verlag. xv, 445 p(1995)hope this helpspeter =I am writing a computer application that performs a Michaelis-Mentenregression. Could anybody help me in obtaining the source code or the algorithmthat performs that regression? Gio =repulsive forces to each other, and can freely move on the torus.Does anyone know any work related to this kind of problems?I am especially interested in questions like: How to nd the minimumenergy conguration?; How to characterize the systems stabilities? etc.James X. Li => repulsive forces to each other, and can freely move on the torus.> Does anyone know any work related to this kind of problems?> I am especially interested in questions like: How to nd the minimum> energy conguration?; How to characterize the systems stabilities? etc.I have done some work on:http://www.members.shaw.ca/goslinga/sphere/sphere.htmIf the minimization problem on the torus is similar to the one on thesphere be prepared for many meta-stable solutions for larger numberof points and for multiple optima for given number of points.I thought about but never tackled the torus since the sphere problemis already very hard. Good Luck!> James X. LiJentje Goslinga information. I have another question:do you know works which address the problem from a statisticalaspect? The system I am trying to study can easily have manythousand points. Approaching this problem analytically (like mostof the reference do) seems to be impossible. I thought about applyingGibbs eld theory, but couldnt work out any useful concept.James X. Li>> repulsive forces to each other, and can freely move on the torus.>> Does anyone know any work related to this kind of problems?>> I am especially interested in questions like: How to nd the minimum> energy conguration?; How to characterize the systems stabilities?etc.>> I have done some work on:>> http://www.members.shaw.ca/goslinga/sphere/sphere.htm>> If the minimization problem on the torus is similar to the one on the> sphere be prepared for many meta-stable solutions for larger number> of points and for multiple optima for given number of points.> I thought about but never tackled the torus since the sphere problem> is already very hard. Good Luck!>>> James X. Li>> Jentje Goslinga> the very helpful information. I have another question:> do you know works which address the problem from a statistical> aspect?Sorry, no. I believe some of the earlier analytical methods userandom search vectors and there is a statistical aspect to these.Imagine a conguration which is almost optimal with only a singlenegative direction (energy problem). There must be an epsilonenvironment of negative search directions but nding a negativesearch vector using random search may become statistically difcultand practically impossible as the order increases.Andrew Sherwood is using a topological approach to obtain certainsymmetry congurations. He also has a sphere FAQ at:http://www.ogre.nu/sphere.htmHe seems to have a mechanism to generate congurations which isnot based on analytical calculations. You should check with himwhether this is applicable to your order of problem.[I am sending this reply as a bcc to him.]> The system I am trying to study can easily have many> thousand points. Approaching this problem analytically (like most> of the reference do) seems to be impossible.I am working on an ultra fast Divide and Conquer eigensolver inC and expect to be able to crack 300-400 points on a 2GHz+ machinebut thousands of points is clearly out off my reach.These eigensolvers are cubic order and worse in practice due tocache effects and so on. I use methods for constrained optimization.On my web site I have collections of accurate optimum solutionsfor the geometric problem which I believe to be fairly completefor orders to 96. It would be possible to run statistics on those.> I thought about applying> Gibbs eld theory, but couldnt work out any useful concept.> James X. Li>>repulsive forces to each other, and can freely move on the torus.>>Does anyone know any work related to this kind of problems?>>I am especially interested in questions like: How to nd the minimum>energy conguration?; How to characterize the systems stabilities?> etc.>>I have done some work on:>>http://www.members.shaw.ca/goslinga/sphere/sphere.htmIf the \ minimization problem on the torus is similar to the one on the>>sphere be prepared for many meta-stable solutions for larger number>>of points and for multiple optima for given number of points.>>I thought about but never tackled the torus since the sphere problem>>is already very hard. Good Luck!>>James X. Li>>Jentje GoslingaI wish you all the best of luck,Jentje = > Andrew SherwoodAlmost. :P > is using a topological approach to obtain certain > symmetry congurations. He also has a sphere FAQ at: > http://www.ogre.nu/sphere.htmNot a FAQ but a classied collection of links(which include a faq or two).Perhaps this new publicity will attract some new links, eh? > He seems to have a mechanism to generate congurations which is > not based on analytical calculations. . . .What Ive written is a tree-search for topologically-distinct triangulations of the sphere, or equivalently `cubic polyhedra. (It is still quite crude, and lately my attention has been elsewhere.)Rather than moving the vertices around, my code builds an abstract model of the edges.My oldest webpage (http://www.ogre.nu/soccer.html) reects an earlier version, which handled only pentagons and hexagons up to a total of 31 faces (equivalent to triangulations with up to 31 vertices of degree 5 and 6). I believe that code is lost.Im not the rst to do any of this. See for example:Fullerene Structure Galleryhttp://www.cochem2.tutkie.tut.ac.jp/Fuller/ Fuller.htmlCounting Polyhedrahttp://home.att.net/~numericana/data/polyhedra.htm-- Anton Sherwood, http://www.ogre.nu/ => and repulsive forces to each other, and can freely move on the> torus. Does anyone know any work related to this kind of problems?> I am especially interested in questions like: How to nd the> minimum energy conguration?; How to characterize the systems> stabilities? etc.In a molecular system the nearest minimum can be found by usingconjugate gradient to minimize the poential energy The followingfurther searchingThese guys found some minima of an LJ liquid[10] F.H. Stillinger and T.A. Weber, Phys. Rev. A 25, 978 (1982)energy minima using conjugate gradientSchrder T.B., and Dyre J.C., Hopping in a Supercooled Binary Lennard-Jones Liquid,J. Non-Cryst. Solids 235-237(1-3), 331-334 (1998). stationary pointsAuthors: L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala, F. SciortinoComments: 4 pages, 3 gures, to be published on PRLSubj-class: Disordered Systems and Neural NetworksJournal-ref: Phys. Rev. Lett. 85, 5356 (2000)See alsoDoye and Wales j chem phys 116, 3777 I hope that this helps Nieks L Ellegaard -- Niels L Ellegaard http://dirac.ruc.dk/~gnalle/ information. I have another question: do> you know works which address the problem from a statistical aspect?> The system I am trying to study can easily have many thousand> points. Approaching this problem analytically (like most of the> reference do) seems to be impossible. I thought about applying Gibbs> eld theory, but couldnt work out any useful concept.For nding all the mnima of a small system. Here are random startingreference)Angelani L, Parisi G, Ruocco G, et al.PHYS REV LETT 81 (21): 4648-4651 NOV 23 Doye and Wales J Chem Phys 116, p 3777At low temperature you may be able to use Bolzmann factors to rescaleremember correctly they do soemthing similar or they may refer topeople that doHeuer A, Buchner SWhy is the density of inherent structures of a Lennard-Jones-typesystem Gaussian?J PHYS-CONDENS MAT 12 (29): 6535-6541 JUL 24 2000-- Niels L Ellegaard http://dirac.ruc.dk/~gnalle/ helpful information. I have another question:> do you know works which address the problem from a statistical> aspect? The system I am trying to study can easily have many> thousand points. Approaching this problem analytically (like most> of the reference do) seems to be impossible. I thought about applying> Gibbs eld theory, but couldnt work out any useful concept.I have recently worked on clustering on a n-torus. Since attractivexed points are associated to clusters in data space, you can getyour positions on the torus by clustering.A basic reference for the relation of xed points to clusters isPhys.Rev.E. 61(5),4691 (2000) ,while you nd a reference for clustering on a n-torus athttp://www.wias-berlin.de/publications/preprints/843/Hope that helpsAxel information. I have another question:> do you know works which address the problem from a statistical> aspect? The system I am trying to study can easily have many> thousand points. Approaching this problem analytically (like most> of the reference do) seems to be impossible. I thought about applying> Gibbs eld theory, but couldnt work out any useful concept.I have recently worked on clustering on a n-torus. Since attractivexed points are associated to clusters in data space, you can getyour positions on the torus by clustering.A basic reference for the relation of xed points to clusters isPhys.Rev.E. 61(5),4691 (2000) ,while you nd a reference for clustering on a n-torus athttp://www.wias-berlin.de/publications/preprints/843/Hope that helpsAxel =>>and repulsive forces to each other, and can freely move on the>>torus. Does anyone know any work related to this kind of problems?>>I am especially interested in questions like: How to nd the>>minimum energy conguration?; How to characterize the systems>>stabilities? etc.> In a molecular system the nearest minimum can be found by using> conjugate gradient to minimize the poential energy The following ^^^^^^^^^^^^^^^^^^> further searching> These guys found some minima of an LJ liquidYou have no idea what you are talking about:Firstly, even for unconstrained problems with widely disparateeigenvalues, simplistic rst order methods like CG will notwork for problems involving many variables.The original poster mentions thousands of points which translatesinto three times as many variables.Secondly, this is a constrained optimization problem not anunconstrained one:If it is formulated in terms of cartesian coordinates there are3*n variables and n+1 constraints, the requirement of the pointsbeing on the surface of the torus and 1 rotation constraint.If it is formulated in terms of the coordinates of the surface ofthe torus, there are 2*n variables and still 1 rotation constraint.[In this case absolute second derivatives should be used in secondorder methods.]For the sphere problem which has 3 rotation constraints I havefound the solutions for order up to order 96 using the secondorder method which I describe in my paper on the sphere problem.Of course I have tried gradient methods but they just do notwork for this order of problem and they seem to converge to onlyseven decimals after ten thousands of iterations.My results are about 10-11 decimals accurate.I also describe a systematic approach towards nding all theoptima of a problem.Note also that even using methods for constrained optimizationthere are still optimum solutions for which the constrainedproblem has zero eigenvalues because the solutions admitsadditional degrees of freedom.> I hope that this helpsYour misinformation is of no help.> Nieks L Ellegaard > Jentje Goslinga =usual scenario: The coefcients of p(x) in form (1) are> rounded real numbers; in this case, multiple roots are a myth,> and clusters of simple roots are a reality with probability 1.> (Find out about Weierstasss Preparation Theorem which tells you> how the perturbations of multiple roots travel in a> reworks-like fashion from their centre.)>If the polynomial does have multiple roots, but the coefcients arerounded, the roots can still be calculated accurately. For details, see myrecently paper Computing multiple roots of inexact polynomials(x-1)^20 is rather easy..>> => If the polynomial does have multiple roots, but the coefcients are> rounded, the roots can still be calculated accurately. But how do you know that a root is multiple and not several close ones?Arnold Neumaier =Thats up to the user to decide. My method tells you how close yourpolynomial is to a polynomial having exact multiple roots, what is thesensitivity of the results, what is the estimated forward error, etc.For example, if you expand p(x) = (x-1)^20 * (x-2)^15 * (x-3)^10 * (x-4)^5and round the coefcients to 16 digits using IEEE standard doubleprecision. Then send the polynomial to my software MultRoot: >> z = multroot(p)THE PEJORATIVE CONDITION NUMBER: 76.7706THE BACKWARD ERROR: 8.80e-016THE ESTIMATED FORWARD ROOT ERROR: 1.35e-013 computed roots multiplicities 4.000000000000001 5 3.000000000000003 10 1.999999999999998 15 1.000000000000000 20In other words, the given (inexact) polynomial is 8.8e-16 away from an exactpolynomial with above exact roots, andthe condition number is 76 (extremely stable!). The simple rootscalculated by Matlab have a conditition number somewhere near 1e+15.So, which roots do you trust? It is totally up to you.Zhonggang Zeng>> If the polynomial does have multiple roots, but the coefcients are> rounded, the roots can still be calculated accurately.>> But how do you know that a root is multiple and not several close ones?>> Arnold Neumaier => Thats up to the user to decide. My method tells you how close your> polynomial is to a polynomial having exact multiple roots, what is the> sensitivity of the results, what is the estimated forward error, etc.> For example, if you expand> p(x) = (x-1)^20 * (x-2)^15 * (x-3)^10 * (x-4)^5> and round the coefcients to 16 digits using IEEE standard double> precision. Then send the polynomial to my software MultRoot:>> z = multroot(p)> THE PEJORATIVE CONDITION NUMBER: 76.7706> THE BACKWARD ERROR: 8.80e-016> THE ESTIMATED FORWARD ROOT ERROR: 1.35e-013> computed roots multiplicities> 4.000000000000001 5> 3.000000000000003 10> 1.999999999999998 15> 1.000000000000000 20> In other words, the given (inexact) polynomial is 8.8e-16 away from an exact> polynomial with above exact roots, and> the condition number is 76 (extremely stable!). The simple roots> calculated by Matlab have a conditition number somewhere near 1e+15.> So, which roots do you trust? It is totally up to you.The le README.txt of your matlab implementation available fromhttp://www.neiu.edu/~zzeng/zroot.htmsays:At the current version, the code works for polynomials withmoderate root multiplicities (such as 30).I tried to be even more moderate, searching only for double roots(except in one case, after I suspected that the method needs high multiplicity to work...)p=poly([2:8,2]);z = multroot(p)did not locate the double root at 2 but gave a success message, allmultiplicities equal to 1.p=poly([1:6,1:6]);z = multroot(p)did not nd any of the six double root but returned an error message(The computation is unsuccessful ...).The same happened for the sixfold rootsn=5;p=poly([1:n,1:n,1:n,1:n,1:n,1:n]);z = multroot(p)with lengthy warnings returned in addition.p=poly([1:15,7.9999]);z = multroot(p)found a double root at 14.696665562196047-0.000000000381809 ialthough there is no root near 14.69. For comparison, p=poly([1:15,7.9999]);format short;z = roots(p)with Matlabs built in root nder got all roots correct to several decimal places: 15.0000 14.0000 12.9999 12.0002 10.9997 10.0003 8.9997 8.0092 7.9909 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000Thus there is no obvious reason to trust your roots.Arnold Neumaier =Your examples are very good. They are in spirit of the Wilkinsonspolynomial that are ill-conditioned in the sense that those polynomials arenear several polynomials with different multiplicity structure. In myfollow-up work, Ill try to address some of these issues. Even in currentversion, there are control parameters you can adjust to get differentresults. You are absolutely right, there is no obvious reason to trust anyset of roots more than the others. My method at least provides some optionsother methods do not.>> Thats up to the user to decide. My method tells you how close your> polynomial is to a polynomial having exact multiple roots, what is the> sensitivity of the results, what is the estimated forward error, etc.>> For example, if you expand> p(x) = (x-1)^20 * (x-2)^15 * (x-3)^10 * (x-4)^5> and round the coefcients to 16 digits using IEEE standard double> precision. Then send the polynomial to my software MultRoot:>> z = \ multroot(p)>> THE PEJORATIVE CONDITION NUMBER:76.7706> THE BACKWARD ERROR: 8.80e-016> THE ESTIMATED FORWARD ROOT ERROR: 1.35e-013>> computed roots multiplicities>> 4.000000000000001 5> 3.000000000000003 10> 1.999999999999998 15> 1.000000000000000 20>> In other words, the given (inexact) polynomial is 8.8e-16 away from anexact> polynomial with above exact roots, and> the condition number is 76 (extremely stable!). The simple roots calculated by \ Matlab have a conditition number somewhere near 1e+15.>> So, which roots do you trust? It is totally up to you.> The le README.txt of your matlab implementation available from> http://www.neiu.edu/~zzeng/zroot.htm> says:> At the current version, the code works for polynomials with> moderate root multiplicities (such as 30).>> I tried to be even more moderate, searching only for double roots> (except in one case, after I suspected that the method needs> high multiplicity to work...)> p=poly([2:8,2]);z = multroot(p)> did not locate the double root at 2 but gave a success message, all> multiplicities equal to 1.>> p=poly([1:6,1:6]);z = multroot(p)> did not nd any of the six double root but returned an error message> (The computation is unsuccessful ...).>> The same happened for the sixfold roots> n=5;p=poly([1:n,1:n,1:n,1:n,1:n,1:n]);z = multroot(p)> with lengthy warnings returned in addition.>> p=poly([1:15,7.9999]);z = multroot(p)> found a double root at 14.696665562196047-0.000000000381809 i> although there is no root near 14.69. For comparison,> p=poly([1:15,7.9999]);format short;z = roots(p)> with Matlabs built in root nder got all roots correct to> several decimal places:> 15.0000> 14.0000> 12.9999> 12.0002> 10.9997> 10.0003> 8.9997> 8.0092> 7.9909> 7.0000> 6.0000> 5.0000> 4.0000> 3.0000> 2.0000> 1.0000>> Thus there is no obvious reason to trust your roots.>> Arnold Neumaier =Although for computing mean, we always normalize summation by N, inmany statistical literatures I see two different variations forcomputing Covariance Matrix, one is normalized by N and the other byN-1. I was wondering why it is so and what is the difference and whichone is a better estimation, say for Covariance Matrix estimation formodeling data by a normal distribution.Thnx--Hossein Mobahi =True and unbiased variance:http://www.visualstatistics.net/Workbook/True_and_ Unbiased_Variance.htmZdenek Hurak>> Although for computing mean, we always normalize summation by N, in> many statistical literatures I see two different variations for> computing Covariance Matrix, one is normalized by N and the other by> N-1. I was wondering why it is so and what is the difference and which> one is a better estimation, say for Covariance Matrix estimation for> modeling data by a normal distribution.>> Thnx>> --Hossein Mobahi---Odchoz.92 zpr.87va neobsahuje viry.Zkontrolov.87no antivirovm syst.8emem AVG (http://www.grisoft.cz).Verze: 6.0.491 / Virov.87 b.87ze: 290 =in my opinion, the one normalized by N is a maximum likelihoodestimate on the real (unknown) covariance matrix, the second estimate(normalized by N-1) has minimal variance.in my opinion, the one normalized by N is a maximum likelihood>estimate on the real (unknown) covariance matrix, the second estimate>(normalized by N-1) has minimal variance.Minimal expected squared error would require dividing by N+1. Dividing by N-1 gives the unbiased estimator, and also extends to other cases where more complicated proceduresare used, such as the result of a regression on k predictors,where one divides by N-k. Subtracting the sample mean isa regression on one predictor.-- This address is for information only. I do not claim that these viewsare those of the Statistics Department or of Purdue University.Herman Rubin, Deptartment of Statistics, Purdue University