Subject: Re: GCD of a list of polynomial coefficient p[x] = 5x+10; k=CoefficientList[p[x], x]; GCD@@k 5 This is the same as Apply[GCD, k] 5 === > Subject: GCD of a list of polynomial coefficient > I wand to get a GCD of a list of polynomial coefficients but > had no luck, i.e. p[x]=5x+10, k=CoefficientList[p[x],x]={5,10}, > GCD[k]={5,10}. No result. I cant remove the curly braces from > GCD[{...}] in order GCD to work. please help === Subject: Re: Re: Beware of adding 0.0 The range of numbers 0.0 represents (if that were the situation at all) has nothing to do with it, since even ZERO digits of precision is enough to prevent the rrormentioned. p = SetPrecision[314159265358979323, 25] s = SetPrecision[p + 0.0`0, 25] 3.14159265358979323`25.*^17 3.14159265358979323`25.*^17 The roblemoccurs for a very different reason: because adding a machine-precision number to another number, in Mathematica, causes arithmetic to be done in machine floating point. If I write down arithmetic that involves machine-precision numbers, I expect this to occur and I probably want this to occur, because it's the fastest way to do arithmetic. p loses digits from being converted to machine precision, hence the roblem Anyway, the problem is entirely artificial. If I know x to 25 digits and add y, where y is known to achine precision what is the correct precision on an arbitrary machine? Unknown, correct? (Machine precision differs from machine to machine, hence the name.) If x and y have been given specific precisions, Mathematica does the arithmetic just fine, so far as I know. That's one reason for saying the problem is artificial. A second reason is that (IMHO) there are hardly any numbers you'd ever know to 25 digits, unless you also know them PERFECTLY, either because they result from exact expressions or they're defined constants, etc. If you know a number exactly, enter it as such! That may or may not be a valid analysis for the constant p, but you clearly meant 0.0 to be exactly 0 (or why all the complaints?). So enter it that way: p = SetPrecision[314159265358979323, 25] s = SetPrecision[p + 0, 25] 3.14159265358979323`25.*^17 3.14159265358979323`25.*^17 That's at least as easy as the other way of entering it, and it gets you the answer you seem to want. Languages are allowed to have syntax, and Mathematica syntax says 0.0, 0, and 0.0`25 are three different things. Why not just get used to it? Bobby > Paul, since you are copying this from my posting to sci.math.symbolic > I took the example from my own paper, A Review of Mathematica, > which appeared in J. Symbolic Computation,(1992) and a version of which > is online. You might also read Mark Soufroniou's INRIA paper (1999?) > explaining why this behavior is apparently intentional, even if it is > hard to explain (shown by the responses you got). > The idea that 0.0 should represent all numbers in [-1,1]* 2^{-precision} > is generally regarded with skepticism by computational scientists. Some > of the other consequences of this are discussed by Soufroniou and others. > RJF >> In Mathematica 5.0, I get a very strange result when I add 0.0 to a number, >> why is this? See below: >> p = SetPrecision[314159265358979323, 25] >> q = SetPrecision[314159265358979323., 25] >> r = SetPrecision[314159265358979323.00000000000000000000, 25] >> s = SetPrecision[p + 0.0, 25] >> {p == q, q == r, r == s, r == p} >> {Tan[s], N[Tan[p]], Tan[q], Tan[r]} >> !(3.14159265358979323`25.*^17) >> Out[47]= >> !(3.14159265358979323`25.*^17) >> Out[48]= >> !(3.14159265358979323`25.*^17) >> Out[49]= >> !(3.14159265358979328`25.*^17) *****NOTICE the 28 instead of 23 at the >> end, why does this happen? ******* >> Out[50]= >> {True, True, False, True} >> Out[51]= >> {1.599808, -1.12979, -1.1297927, -1.1297927} -- DrBob@bigfoot.com www.eclecticdreams.net === Subject: Re: Beware of adding 0.0 helping us understand something, why is it that when we give the following input: In[69]:= p = 314159265358979323 q = 314159265358979323. r = 314159265358979323.00000000000000000000; s = p + 0.00000000000000000000; {p, q, r, s} {N[Tan[p]], N[Tan[q]], N[Tan[r]], N[Tan[s]]} And get the following output: Out[69]= 314159265358979323 Out[70]= !(3.14159265358979323`17.497149872694134*^17) Out[73]= !({314159265358979323, 3.14159265358979323`17.497149872694134*^17, 3.14159265358979323`37.49714987269414*^17, 3.141592653589793`*^17}) Out[74]= {1.60125, ComplexInfinity, -1.12979, 1.60125} Then, the output for Tan[p] and Tan[q] give two different answers? Even though they are almost exactly the same numbers (only different in precision), If it can calculate the Tan[314159265358979323] then why is it not able to give an approximate value for Tan[3.14159265358979323`17.497149872694134*^17], Also, if it calculates Tan[314159265358979323] to be 1.60125, then why does it give Tan[3.14159265358979323`37.49714987269414*^17] with such a severe variation (-1.12979)? > Paul, since you are copying this from my posting to sci.math.symbolic > I took the example from my own paper, A Review of Mathematica, > which appeared in J. Symbolic Computation,(1992) and a version of which > is online. You might also read Mark Soufroniou's INRIA paper (1999?) > explaining why this behavior is apparently intentional, even if it is > hard to explain (shown by the responses you got). > The idea that 0.0 should represent all numbers in [-1,1]* 2^{-precision} > is generally regarded with skepticism by computational scientists. Some > of the other consequences of this are discussed by Soufroniou and others. > RJF > In Mathematica 5.0, I get a very strange result when I add 0.0 to a number, > why is this? See below: > p = SetPrecision[314159265358979323, 25] > q = SetPrecision[314159265358979323., 25] > r = SetPrecision[314159265358979323.00000000000000000000, 25] > s = SetPrecision[p + 0.0, 25] > {p == q, q == r, r == s, r == p} > {Tan[s], N[Tan[p]], Tan[q], Tan[r]} > !(3.14159265358979323`25.*^17) > Out[47]= > !(3.14159265358979323`25.*^17) > Out[48]= > !(3.14159265358979323`25.*^17) > Out[49]= > !(3.14159265358979328`25.*^17) *****NOTICE the 28 instead of 23 at the > end, why does this happen? ******* > Out[50]= > {True, True, False, True} > Out[51]= > {1.599808, -1.12979, -1.1297927, -1.1297927} === Subject: Re: Beware of NSolve Old version 2.2 does not report roots 1 and 4 by NSolve. May be later on improvements (complex roots) have associated problems, may be elimination of extraneous roots by plotting is a necessary check. > Run v. 4.2 on Mac: > f=5/432 - 11/(27*Sqrt[70]*Sqrt[19 - 1890*x]) + x/(2*Sqrt[38/35 - 108*x]); > Solve[f==0,x] returns 2 real roots: > {{x -> (-171 - 25*Sqrt[105])/30240}, {x -> (-171 + 25*Sqrt[105])/30240}} > NSolve[f==0,x] returns 4 real roots: > {{x -> -0.10481082961146104}, {x -> -0.014126116704662378}, > {x -> 0.002816592895138556}, {x -> 0.0003796126802330315}} > Roots 1 and 4 are incorrect. (Just plot f) > Had a similar problem with a quartic 3 months ago. This is a > simpler example. === Subject: Re: minimal power S=A t^(-5/2)+B t^(-3/2)+C t^(-1/2)+D t^(1/2)+F t^(3/2); Exponent[S,t,Min] -(5/2) Min[Cases[S, a_.*t^x_:>x]] -(5/2) === > Subject: minimal power > Hi ! > I have for example following expression: > S = A t^(-5/2) + B t^(-3/2) + C t^(-1/2) + D t^(1/2) + F t^(3/2); > Eith command xponent[S,t]I can receive maximal power of S. > Why I may calculate minimal power of S ? > Nodar Shubitidze > Joint Institute for Nuclear Research > Dubna, RUSSIA === Subject: Re: Simply derivative question, Math 5. The input form of the exponential constant is E not e D[(1 + c*e^t)/(1 - c*e^t),t]//Simplify (2*c*e^t*Log[e])/(c*e^t - 1)^2 D[(1 + c*E^t)/(1 - c*E^t),t]//Simplify (2*c*E^t)/(c*E^t - 1)^2 === > Subject: Simply derivative question, Math 5. > Hello- I am trying to take the derivative of the following function: > (1 + c*e^t)/(1 - c*e^t) > with respect to t. > !(Simplify[D[(1 + c e^t)/(1 - c e^t), t]]) > The answer I get includes an xtraLog[e] in the numerator. Am I not > using the program correctly or am I not understanding the answer > correctly? Other? > -- === Subject: Re: A functional measure of roughness The measure of roughness based on the 3 point Bezier which is like a Lyapunov exponent average, but more sensative: It turns out to be like a second derivative: Measure[n]=Sum[Log[1+Abs[f''(x(i))]/4],{i,1,n}]/n Limit[Measure[n],n-> Infinity]=rho The roughness measure is relate to the Lyapunov exponent average by k=Sum[Log[f'(x(i))],{i,1,n}]/n/Sum[Log[1+Abs[f''(x(i))]/4],{i,1,n}]/n k is close to 2 for the primes. So I have actually found two measures. > In thinking of a way to get a better than Lyapunov , Hausdorff or Kolmogorov > measure of dimension , I thought of this: > F(curve)=0 if smooth and continuous > F(curve)<>0 if rough or discontinuous > The best measure of dimensional roughness (Mandelbrot's way of > expressing it) is the > Lyapunov exponent (or maybe the Hurst exponent?). > Box counting or capacity/ entropy dimension of the Kolmogorov type > is too big most of the time > while Hausdorff being very cut-off measure like > is usually too small. > The trouble with Lyapunov is that it depends on a derivative > and unless you are talking about a fractional derivative, > many fractal functions are of the Weierstrass fractal type > where the classical derivative doesn't exist. > I did some work on Bezier functions in IFS in the past > and fractional partial derivatives of an angular sort as well. > I came to realize that the three point Bezier function of an iterative > sequence in n: > Bezier[p,n]=p^2*f(n+2+2*p*(1-p)*f(n+1)+(1-p)^2*f(n) > is such that if smooth and continuous: > f(n+1)=Bezier[1/2,n]=f(n+2)/4+f(n+1)/2+f(n)/4 > So that the function : > delta[n]=f(n+2)/4+f(n+1)/2+f(n)/4-f(n+1) > is a measure of the roughness. > Putting this measure in an Lyapunov average type function: > Measure[n]=Sum[Log[1+delta[i]],{i,1,n}]/n > I tried this out by comparing it to a known rough set, the primes > and it's Lyapunov integer difference average. > In this experiment the new Bezier roughness measure performs better than the > Lyapunov equivalent over the same range in detecting roughness. > Respectfully, Roger L. Bagula > tftn@earthlink.net, 11759Waterhill Road, Lakeside,Ca 92040-2905,tel: > 619-5610814 : > URL : http://home.earthlink.net/~tftn > URL : http://victorian.fortunecity.com/carmelita/435/ Respectfully, Roger L. Bagula tftn@earthlink.net, 11759Waterhill Road, Lakeside,Ca 92040-2905,tel: 619-5610814 : URL : http://home.earthlink.net/~tftn URL : http://victorian.fortunecity.com/carmelita/435/ === Subject: Re: Simply derivative question, Math 5. That's because your e is NOT the base of natural logarithms (which is E). Steve Luttrell > Hello- I am trying to take the derivative of the following function: > (1 + c*e^t)/(1 - c*e^t) > with respect to t. > !(Simplify[D[(1 + c e^t)/(1 - c e^t), t]]) > The answer I get includes an xtraLog[e] in the numerator. Am I not > using the program correctly or am I not understanding the answer > correctly? Other? > -- === Subject: Re: Simply derivative question, Math 5. Ted, Use capital E instead of small e. Simplify[D[(1 + c*E^t)/(1 - c*E^t), t]] (2*c*E^t)/(-1 + c*E^t)^2 David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Hello- I am trying to take the derivative of the following function: (1 + c*e^t)/(1 - c*e^t) with respect to t. !(Simplify[D[(1 + c e^t)/(1 - c e^t), t]]) The answer I get includes an xtraLog[e] in the numerator. Am I not using the program correctly or am I not understanding the answer correctly? Other? -- === Subject: Re: Simply derivative question, Math 5. As the Mathematica Book (printed or electronic) explains, ALL built-in objects in Mathematica have names that begin with an upper-case letter. So the Napierian is represented by E, not e. (Or, if you want it to look more like conventional mathematical notation, type Esc e Esc to obtain the Mathematica stylized e; in fact, you'll see that stylized e in the result.) So Mathematica was treating e as an unknown constant and gave the correct answer. When e = E, then Log[e] = 1, so no xtrafactor: D[(1 + c E^t)/(1 - c E^t), t] // Simplify By the way, you don't need an explicit multiplication operator after the c here -- just a space (so Mathematica can see you have two single-letter names rather than a two-letter name cE). It's REALLY worth a couple hours of your time to read the early material in the Mathematica Book! It will save you a lot of time avoiding such basic difficulties. > Hello- I am trying to take the derivative of the following function: > (1 + c*e^t)/(1 - c*e^t) > with respect to t. > !(Simplify[D[(1 + c e^t)/(1 - c e^t), t]]) > The answer I get includes an xtraLog[e] in the numerator. Am I not > using the program correctly or am I not understanding the answer > correctly? Other? -- Murray Eisenberg murray@math.umass.edu Mathematics & Statistics Dept. Lederle Graduate Research Tower phone 413 549-1020 (H) University of Massachusetts 413 545-2859 (W) 710 North Pleasant Street fax 413 545-1801 Amherst, MA 01003-9305 === Subject: Re: Simply derivative question, Math 5. Both these answers are correct: D[(1 + c*E^t)/(1 - c*E^t), t] (c*E^t)/(1 - c*E^t) + (c*E^t*(1 + c*E^t))/(1 - c*E^t)^2 D[(1 + c*e^t)/(1 - c*e^t), t] (c*e^t*Log[e])/(1 - c*e^t) + (c*e^t*(1 + c*e^t)*Log[e])/(1 - c*e^t)^2 Built-in constants start with capitals; hence e is undefined (unless you define it). Bobby > Hello- I am trying to take the derivative of the following function: > (1 + c*e^t)/(1 - c*e^t) > with respect to t. > !(Simplify[D[(1 + c e^t)/(1 - c e^t), t]]) > The answer I get includes an xtraLog[e] in the numerator. Am I not > using the program correctly or am I not understanding the answer > correctly? Other? -- DrBob@bigfoot.com www.eclecticdreams.net === Subject: Re: minimal power odar Shubitidze schrieb im Newsbeitrag > Hi ! > I have for example following expression: > S = A t^(-5/2) + B t^(-3/2) + C t^(-1/2) + D t^(1/2) + F t^(3/2); > Eith command xponent[S,t]I can receive maximal power of S. > Why I may calculate minimal power of S ? > Nodar Shubitidze > Joint Institute for Nuclear Research > Dubna, RUSSIA What about Limit[Log[t, S], t -> 0] -5/2 ?? -- Peter Pein, Berlin to write to me, start the subject with [ === Subject: Re: minimal power Firstly, it is not a good idea to use D as a variable name because it is reserved for doing differentiation. The safe bet is to use variable names that start with a lowercase character. Define your series (you can use e rather than f with no problems) s = a/t^(5/2) + b/t^(3/2) + c/t^2^(-1) + d*t^(1/2) + e*t^(3/2); If this were a polynomial then CoefficientList[s, t] would extract the various powers, but here you have to be cleverer. There are lots of ways to do what you want, but one that will be generally useful here and elsewhere is to evaluate FullForm[s], which gives Plus[Times[a, Power[t, Rational[-5, 2]]], Times[b, Power[t, Rational[-3, 2]]], Times[c, Power[t, Rational[-1, 2]]], Times[d, Power[t, Rational[1, 2]]], Times[e, Power[t, Rational[3, 2]]]] This tells you that the various powers of t are represented as Power[t, Rational[x_, y_]], where x_ and y_ are patterns that match to any arguments of Rational. Now use this to extract the various powers of t in s (extracting only the powers and omitting the base t): p = Cases[s, Rational[x_, y_], 3]; This gives p as the list {-(5/2), -(3/2), -(1/2), 1/2, 3/2}. {Min[p],Max[p]} then gives you the minimum and maximum powers of t. This approach can be generalised. You start by looking at FullForm[expression] and then design an appropriate Cases expression that pulls out the pieces that you want. Steve Luttrell > Hi ! > I have for example following expression: > S = A t^(-5/2) + B t^(-3/2) + C t^(-1/2) + D t^(1/2) + F t^(3/2); > Eith command xponent[S,t]I can receive maximal power of S. > Why I may calculate minimal power of S ? > Nodar Shubitidze > Joint Institute for Nuclear Research > Dubna, RUSSIA === Subject: Re: minimal power Nodar, S = A t^(-5/2) + B t^(-3/2) + C t^(-1/2) + D t^(1/2) + F t^(3/2); Exponent[S, t, Min] -(5/2) Check out Exponent in Help. David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Hi ! I have for example following expression: S = A t^(-5/2) + B t^(-3/2) + C t^(-1/2) + D t^(1/2) + F t^(3/2); Eith command xponent[S,t]I can receive maximal power of S. Why I may calculate minimal power of S ? Nodar Shubitidze Joint Institute for Nuclear Research Dubna, RUSSIA === Subject: Re: 2005 Geonetry Calendars The website doesn't show nearly enough of the calendar to make me think, or even suspect, that I'd want to buy it. Bobby > I have made two high quality math calendars for 2005, > one with minimal surfaces and the other with polyhedra. > They are available at > http://www.lulu.com/content/66100/ > http://www.lulu.com/content/66755/ > === Subject: Re: Technical Publishing Made Easy with New Wolfram Publicon Software This appears to be an elaborate waste of binary bits. Rather than make Mathematica do pagination right (and a few other simple things), they made a new stand-alone LaTex derivative with no computational capability. MUCH of the content I'd likely put into Publicon, if I used it, would originate in Mathematica. But conversion is a one-way street. Note that Publicon doesn't support footnotes; something every word processor does do, and something every technical document needs. On the PLUS side, it's cheap--except in terms of the learning curve. The online tour makes using it look very involved. Bobby > Technical Publishing Made Easy with New Wolfram Publicon > Software > Wolfram Publicon, a powerful new publishing tool based on the > underlying document technology of Mathematica, is now available > to purchase as a download for Windows and Mac OS X. > Created for the growing number of academic researchers, > students, and industry professionals who need to create > precisely formatted technical documents in XML and other > structured data formats, Publicon incorporates many exciting > features including inline math and chemistry typesetting, > publisher-specific style sheets, and a scrolling WYSIWYG > interface ideal for online presentation. > With Publicon, users can compose more engaging technical > documents that intuitively incorporate complex scientific > research. Mathematica users will especially appreciate > Publicon's unique ability to understand and identify math. All > Mathematica work, including dynamic 2D and 3D plots, can be > pasted directly into Publicon documents. Publicon will preserve > the mathematical content so the work may be evaluated at any > time in Mathematica. > Heralded as a ajor advanceby Open Access publisher BioMed > Central, Publicon was built to take the guesswork and hassle out > of formatting technical documents for publication. Combining > ease of use with cutting-edge technology, Publicon is the first > choice for composing structured technical documents for > electronic or print publication. > For more information, please visit: > http://www.wolfram.com/publicon === Subject: Re: Technical Publishing Made Easy with New Wolfram Publicon Software Does anyone know if Publicon comes with serious documentation or a tutorial for technical wordprocessing? Something along the line of the many TeX tutorials. > Technical Publishing Made Easy with New Wolfram Publicon > Software > Wolfram Publicon, a powerful new publishing tool based on the > underlying document technology of Mathematica, is now available > to purchase as a download for Windows and Mac OS X. > Created for the growing number of academic researchers, > students, and industry professionals who need to create > precisely formatted technical documents in XML and other > structured data formats, Publicon incorporates many exciting > features including inline math and chemistry typesetting, > publisher-specific style sheets, and a scrolling WYSIWYG > interface ideal for online presentation. > With Publicon, users can compose more engaging technical > documents that intuitively incorporate complex scientific > research. Mathematica users will especially appreciate > Publicon's unique ability to understand and identify math. All > Mathematica work, including dynamic 2D and 3D plots, can be > pasted directly into Publicon documents. Publicon will preserve > the mathematical content so the work may be evaluated at any > time in Mathematica. > Heralded as a ajor advanceby Open Access publisher BioMed > Central, Publicon was built to take the guesswork and hassle out > of formatting technical documents for publication. Combining > ease of use with cutting-edge technology, Publicon is the first > choice for composing structured technical documents for > electronic or print publication. > For more information, please visit: > http://www.wolfram.com/publicon Garry Helzer gah@math.umd.edu === Subject: Re: 3D Screensavers for Mathematica users I tried one of the ultimonitorscreen savers, and it wasn't. Multimonitor, that is. The left screen froze, and the right screen didn't know it. This is very typical of the screen-savers I've seen all over the Internet. None of them work. Bobby > Let me introduce new Multimonitor 3D Math Screensavers at > http://www.graphtree.com/ > The Screensavers are designed for Windows 98/ME/XP, with automatic > install and uninstall. > The Screensavers are made using mathematical objects. These objects > are constructed in your video memory card (no objects are loaded from > your hard drive). Thus the screensavers are very small in size and > efficient. > The Screensavers are multimonitor compatible - up to 5 monitors. > There are the following new Math Screensavers: > 1. 3DThings - Infinit number of 3D Math objects > 2. MathLogo - One of the best screensaver for Mathematica 4 and > Mathematica 5 users! > Besides that for each screensaver appropriate mathematical desktop > backgrounds (wallpapers) are available. === Subject: Conditonal sum I am new to mathematica and could really need a tip on how to implement a conditional sum. Pseudocode would be something like this EndSum=a; if j in Sm: EndSum=EndSum+f[j] end I have tried doing this by using an If-clause inside the sum-function, but this seems little elegant and is nor very suitable for symbolic calculations. Any advice would be highly appreciated Christopher Grinde === Subject: Re: Function Solve from its Property > Regarding versatililty in handling functions. Starting with a > characteristic property definition of any function, how does > Mathematica e.g., Solve[f[x_] +or- f[y_]==f[x *or/ y]] to provide its > solution f[u_]=C Log[u] ? Jos8e Manuel Guti8errez, and Andr8es Iglesias in The Mathematica Journal 5(1), 1995: 82-86. The associated package, which I have not tested, is available at http://cdec.unican.es/ftp/herramientas/fsolve/Packages/ Paul -- Paul Abbott Phone: +61 8 9380 2734 School of Physics, M013 Fax: +61 8 9380 1014 The University of Western Australia (CRICOS Provider No 00126G) 35 Stirling Highway Crawley WA 6009 mailto:paul@physics.uwa.edu.au AUSTRALIA http://physics.uwa.edu.au/~paul === Subject: Re: Do-loop conversion > I have a problem that involves converting a Do-loop structure into > faster-executing version using functional programming. This is an optics > problem -- specifically, calculating the angle dependence of the X-ray > reflectivity of an n-layered medium on an infinitely thick substrate. The > calculation relies on a well-known recursion relation: > X(j) = [ r(j) + X(j+1) p^2(j+1) ]/[1+ r(j) X(j+1) p^2(j+1)] > where X(j) is the ratio of the reflected and transmitted X-ray amplitudes > at the bottom of layer j, p(j+1) is a phase factor (given below), and r(j) > is the Fresnel coefficient for reflection at the interface between layers j > and j+1: > r(j) = [k(z,j) - k(z,j+1)]/[k(z,j) + k(z,j+1)] > where the wavevector k(z,j) = (2*Pi/wavelength)Sqrt[n^2(j) - Cos(phi)^2], You have k(z,j) on the left-hand side but wavelength and phi on the right-hand side. To turn this into a Mathematica function the relationship between (implicit) variables needs to be made clear. > n(j) is the complex index of refraction for X-rays for the layer n(j) = 1 - > delta(j) + I*beta(j), and phi is incident angle of the X-rays. The phase > factor mentioned above is given by > p(j) = r(j) exp[-2 (k(j) k(j+1)) s(j+1)] with s being the roughness at > interface j+1. Shouldn't these be k(z,j) and k(z,j+1)? > The recursion layer works because with an infinitely thick substrate, the > reflected amplitude coming up from the substrate = 0, so at a given angle > of incidence, you work from the bottom up starting with X(bottom) = 0 and > add the amplitudes at each interface until you get to the top surface. Entering the recurrence formula for X(j), x[j_] := (x[j + 1] p[j + 1]^2 + r[j])/(r[j] x[j + 1] p[j + 1]^2 + 1) and the bottom condition, X[5] = 0; you can compute the results for intermediate layers, say Factor[X[3]] or for the top layer (depending on your numbering system), Factor[X[1]] in terms of r[j] and p[j]. These results are valid for _any_ angle phi and are not all that complicated. You could use dynamic programming (x[j_] := x[j] = ...) to make this more efficient, if required. > The various functions above -- wavetabc, leftvectmapc, thtestcomp, > roughrevcomp, coeffcomp, phasecomp, fcoeffcomp, intensitycomp, and > newdatacomp -- are complied functions that take structural parameters for > each layer (index of refraction parameters delta and beta, layer thickness z, I assume? > interfacial roughness) by reading them out of a list {params} that is input > at the start of the program. Both r[j] and p[j] depend (implicitly) on z, lambda, and phi. Computing these functions should be very fast for specific parameter values. > As I said, the above Do-loop works just fine by calculating at each angle > phi the final intensity at that angle, normalizing it with a background and > maximum peak factor (also in the {params} list and applied by the > newdatacomp line) and appending it to the intensity datafile > newc. Executing this loop for a five-layer system and a couple of hundred > angular points takes around 0.2 seconds -- not bad, except that I need to > use this in a genetic algorithm minimization routine which needs many > hundreds of these curves to be calculated with GA- dictated changes to the > initial parameter matrix. So: if anyone can give me some guidance on how to > eliminate the Do-loop over the angular range so I can speed things up, I > would be very grateful. If I understand what you're trying to do, I think that you can leave phi as a symbolic parameter. It appears that all other parameters will be numerical so that the final result, although complicated, should be manageable. Paul -- Paul Abbott Phone: +61 8 9380 2734 School of Physics, M013 Fax: +61 8 9380 1014 The University of Western Australia (CRICOS Provider No 00126G) 35 Stirling Highway Crawley WA 6009 mailto:paul@physics.uwa.edu.au AUSTRALIA http://physics.uwa.edu.au/~paul === Subject: Re: Do-loop conversion Do loops are not always the wrong way to go and, in this case, the real problem could be this statement: newc = AppendTo[newc, newdata] In the first place, it can be replaced with simply AppendTo[newc, newdata] The Set statement isn't necessary, since AppendTo already changes the value of its first argument. Secondly, AppendTo is slow. Consider these three timings: n=10000 SeedRandom[1] new1={}; Timing@Do[AppendTo[new1,Random[]],{i,1,n}] 10000 {0.938 Second,Null} SeedRandom[1] new2={}; Timing[Do[new2={new2,Random[]},{i,1,n}];new2=Flatten@new2;] {0.031 Second,Null} SeedRandom[1] Timing[new3=Reap[Do[Sow@Random[],{i,1,n}]][[-1,1]];] {0.016 Second,Null} new1 == new2 == new2 True The second method only works for you if ewdatais a scalar (not a List). The speed difference gets much MORE dramatic as n gets larger. But, if n=200, there may or may not be a noticeable difference. In that case, other Mathgroupers may have ideas for you. Offhand, I doubt that getting rid of the Do loop will speed things up for you in this case. Bobby > Mathematica colleagues, > I have a problem that involves converting a Do-loop structure into > faster-executing version using functional programming. This is an optics > problem -- specifically, calculating the angle dependence of the X-ray > reflectivity of an n-layered medium on an infinitely thick substrate. The > calculation relies on a well-known recursion relation: > X(j) = [ r(j) + X(j+1) p^2(j+1) ]/[1+ r(j) X(j+1) p^2(j+1)] > where X(j) is the ratio of the reflected and transmitted X-ray amplitudes > at the bottom of layer j, p(j+1) is a phase factor (given below), and r(j) > is the Fresnel coefficient for reflection at the interface between layers j > and j+1: > r(j) = [k(z,j) - k(z,j+1)]/[k(z,j) + k(z,j+1)] > where the wavevector k(z,j) = (2*Pi/wavelength)Sqrt[n^2(j) - Cos(phi)^2], > n(j) is the complex index of refraction for X-rays for the layer n(j) = 1 - > delta(j) + I*beta(j), and phi is incident angle of the X-rays. The phase > factor mentioned above is given by > p(j) = r(j) exp[-2 (k(j) k(j+1)) s(j+1)] with s being the roughness at > interface j+1. > The recursion layer works because with an infinitely thick substrate, the > reflected amplitude coming up from the substrate = 0, so at a given angle > of incidence, you work from the bottom up starting with X(bottom) = 0 and > add the amplitudes at each interface until you get to the top surface. The > Mathematica code that works just fine as a Do-loop is: > Do[ > wavevectmap := wavetabc[revdelta, revbeta, phi]; > leftvectmap := leftvectmapc[wavevectmap]; > thtest := thtestcomp[revrough, leftvectmap, wavevectmap]; > maproughrev := roughrevcomp[thtest]; > coeffmap := coeffcomp[wavevectmap, leftvectmap, maproughrev]; > phasemap := phasecomp[wavevectmap, revthick]; > coeffs := fcoeffcomp[coeffmap, phasemap]; > recurr[prior_, {fc_, pf_}] := (fc + (prior*pf^2))/(1 + > (fc*prior*pf^2)); > intensity := intensitycomp[Fold[recurr, 0, coeffs ]]; > newdata = newdatacomp[phi, intensity, params]; > newc = AppendTo[newc, newdata], > {phi, phistart, phiend, phiinc}] > The various functions above -- wavetabc, leftvectmapc, thtestcomp, > roughrevcomp, coeffcomp, phasecomp, fcoeffcomp, intensitycomp, and > newdatacomp -- are complied functions that take structural parameters for > each layer (index of refraction parameters delta and beta, layer thickness, > interfacial roughness) by reading them out of a list {params} that is input > at the start of the program. I've left out these functions for brevity; if > needed, I'll post them as well. > As I said, the above Do-loop works just fine by calculating at each angle > phi the final intensity at that angle, normalizing it with a background and > maximum peak factor (also in the {params} list and applied by the > newdatacomp line) and appending it to the intensity datafile > newc. Executing this loop for a five-layer system and a couple of hundred > angular points takes around 0.2 seconds -- not bad, except that I need to > use this in a genetic algorithm minimization routine which needs many > hundreds of these curves to be calculated with GA- dictated changes to the > initial parameter matrix. So: if anyone can give me some guidance on how to > eliminate the Do-loop over the angular range so I can speed things up, I > would be very grateful. === Subject: Re: Re: Re: FindMinimum and the minimum-radius circle The recent messages about the problem of finding the minimum radius circle are actually a discussion on the numerical aspects of the problem. The following remark is a little bit outside that scope. It is possible to construct that circle directly. I have to admit that I did not follow this thread very carefully, so if someone else already gave a similar construction, I strongly apologize. The idea is that the smallest circle is a circle that passes through three points or passes through two diametrical points (that has been remarked). A little bit more can be said: when the smallest circle passes through three points, then the center of the circle must be inside the triangle. So the construction runs as follows. First construct a circle that contains all points and passes through one point. Then move the center of the circle in the direction of that point until it passes through a second point. When these points are diametric, we are ready. Otherwise, move the center of the circle to the connecting line until it passes through a third point. Then test. If the center of the circle is outside the triangle, move it in the direction of the largest side of the triangle taking care that all points remain inside the circle, and so on. Here is an implementation: smallestcircle[points_] := Block[{m, p1, p2, p3, v, t, t0, pts = points}, m = Mean[points]; p1 = p2 = p3 = points[[Ordering[points, -1, (#1 - m) . (#1 - m) < (#2 - m) . (#2 - m) & ][[1]]]]; While[Length[Union[{p1, p2, p3}]] < 3 || (p1 + p2)/2 == m || (p1 - p2) . (p1 - p2) > (p1 - p3) . (p1 - p3) + (p2 - p3) . (p2 - p3), v = (p1 + p2)/2 - m; t = 1; p = p3; pts = Reap[Scan[If[v . (#1 - p1) < 0, Sow[#1]; t0 = ((m - #1) . (m - #1) - (m - p1) . (m - p1))/(2* v . (#1 - p1)); If[t0 < t, t = t0; p = #1]] & , pts]][[2]]; If[pts == {}, Break[], pts = pts[[1]]]; m = m + t*v; p3 = p; {p1, p2, p3} = If[(p3 - p1) . (p3 - p1) >= (p2 - p1) . (p2 - p1), {p1, p3, p2}, {p3, p2, p1}]]; Circle[m, Norm[m - p1]]] For exact coordinates we find the exact circle. This function is faster than the numerical ones: In[561]:= SeedRandom[1]; points = Partition[Table[Random[], {1000}], 2]; Timing[NMinimize[ Sqrt[Max[((x - #1[[1]])^2 + (y - #1[[2]])^ 2 & ) /@ points]], {x, y}]] Timing[smallestcircle[points]] Out[563]= {4.625*Second, {0.6628153793262672, {x -> 0.48683054240909296, y -> 0.46145830394871523}}} Out[564]= {0.14000000000000057*Second, Circle[{0.48684175891608683, 0.46146943938871915}, 0.6628150360165743]}