This may happen on the first request due to CSS mimetype issues. Try clearing your browser cache and refreshing.

- Home
- Browse
- Projects
- Maxima -- GPL CAS based on DOE-MACSYMA
- Mailing Lists

Update of /cvsroot/maxima/htdocs/docs/macref/share1 In directory sc8-pr-cvs1:/tmp/cvs-serv18373 Added Files: airy.usg asymp-2.html asymp-3.html asymp-4.html asymp.html asymp.usg asympa.html dblint.dem dblint.html dblint.usg dimen.dem dimen.usg facexp.html facexp.usg hypgeo.html hypgeo.usg inteqn.usg intsce.usg invert.usg kach.dem nchrpl.dem qq.dem qq.usg rncomb.dem rncomb.usg sets.html sets.usg Log Message: macref share1 --- NEW FILE: airy.usg --- May 2, 1981 4:52 pm by Leo P. Harten (LPH) The Airy equation diff(y(x),x,2)-x*y(x)=0 has two linearly independent solutions, taken to be Ai(x) and Bi(x). This equation is very popular as an approximation to more complicated problems in many mathematical physics settings. Do LOAD("AIRY"); to get the functions AI(x), BI(x), DAI(x), DBI(x) . The file SHARE1;AIRY FASL (by LPH@...) contains routines to compute the Ai(x), Bi(x), d(Ai(x))/dx, and d(Bi(x))/dx functions. The result will be a floating point number if the argument is a number, and will return a simplified form otherwise. An error will occur if the argument is large enough to cause an overflow in the exponentials, or a loss of accuracy in sin or cos. This makes the range of validity about -2800 to 1.e38 for AI and DAI, and -2800 to 25 for BI and DBI. The GRADEF rules are now known to MACSYMA: diff(AI(x),x)=DAI(x), diff(DAI(x),x)=x*AI(x), diff(BI(x),x)=DBI(x), diff(DBI(x),x)=x*BI(x). The method is to use the convergent Taylor series for abs(x)<3., and to use the asymptotic expansions for x<-3. or x>3. as needed. This results in only very minor numerical discrepancies at x=3. or x=-3. More accuracy can be had if you request from LPH. For details, please see Abramowitz and Stegun's Handbook of Mathematical Functions, section 10.4 (hardcover ed.) and Table 10.11 . To get the floating point Taylor expansions of the functions here, do ev(TAYLOR(AI(x),x,0,9),infeval); for example. Please also try SHARE;BESSEL FASL (by CFFK) for the AIRY function there. Leo P. Harten (LPH) --- NEW FILE: asymp-2.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/asymp.dm2 was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: asymp-3.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/asymp.dm3 was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: asymp-4.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/asymp.dmo was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: asymp.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/asymp.dm1 was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: asymp.usg --- Preliminary Version The ASYMP package June 1, 1982 ASYMP A package for the evaluation of bounds on Feynman Diagrams. by William E. Caswell (WEC@...) Anthony D. Kennedy (ADK@...) Preliminary Version The ASYMP package June 1, 1982 I. Introduction. ASYMP is a package for determining the asymptotic behavior of Feynman integrals. Given a topological description of a Feynman diagram as a set of lines and vertices, together with information about the mass of the virtual particle corresponding to each line and the momentum entering at each external leg, it will tell one the leading asymptotic behavior of that graph as some sets of masses get much larger than others. As this package is very unlikely to be of use to people who are not familiar with Feynman diagrams and other basic aspects of perturbative quantum field theory, we will refrain from describing the basics here and refer the interested reader to any of the standard textbooks on the subject instead. Perhaps this is also the appropriate place to mention the limitations of the package. These are of two kinds, those which are fundamental limitations of the formalism and methods used in the package itself, and those which are just features which could be added easily if they are ever needed. In the first category we stress that the bounds are obtained for individual Feynman graphs, and not for sums of them; in other words the asymptotic behavior of a green's function might be quite different from that of the graphs which contribute to it, because there may be "miraculous" cancellations. Such cancellations occur in many interesting theories, in particular gauge theories, but they are best dealt with by means of Ward identities rather than explicit calculation. Another mathematical limitation is that the actual behavior of a graph is only bounded by the result given -- in reality the graph might have a smaller asymptotic growth: the bounds obtained are usually fairly good, however. In the second class of limitations we should mention that (1) the package currently deals only with boson fields, (2) allows only a trivial dependence of the vertices upon momenta and masses, (3) tries to compute 1/0 for IR divergent graphs [which is honest, in a way], 1 Preliminary Version The ASYMP package June 1, 1982 and (4) returns INF for a UV divergent graph [which is correct]. All of these are simple to generalize in the program, and if one need to get around these limitations, please contact the authors. A slightly harder problem to circumvent is related to point (4) above, namely (5) one cannot currently specify that a UV divergent graph is to be subtracted in a certain way: part of the problem is that there are many different subtraction schemes (minimal, zero momentum Taylor series, etc.) and how to specify which method one wants is not clear, but it would also require a fair amount of thought to make the program renormalize automatically. Any suggestions or comments on this, or any other aspect of ASYMP, would be appreciated. II. Simple Example. The easiest way to see how ASYMP works is to look at the simplest example, the one-loop three-point function in (phi)^3 theory. First of all we must load the ASYMP package into a MACSYMA: (C1) loadfile(asymp,fasl,dsk,share1)$ ASYMP: version of 11:54pm Saturday, 4 July 1981 (C2) graph1:diagram(line(a,b,1,m),line(b,c,2,m),line(c,a,3,mm), extline(a,4,p),extline(b,5,q),extline(c,6,-p-q))$ 1 Loop Diagram (C3) bound(graph1,[[m,p,q],mm,inf]); MM LOG(--) M (D3) ------- 2 MM 2 Preliminary Version The ASYMP package June 1, 1982 First of all, in line (C1) we have loaded up the FASL (compiled) version of the ASYMP package. It identifies itself by telling us the date on which it was born. We then define the desired Feynman diagram as GRAPH1 using the DIAGRAM function. DIAGRAM takes an arbitrary number of arguments, each of which is a pseudo-function describing a part of the graph. Currently, there are two such pseudo-functions, LINE and EXTLINE. Logically enough LINE describes an internal line; if we type LINE(LONDON,PARIS,RHUBARB,5*M[PLANCK]) we are defining a line from a vertex called LONDON to a vertex called PARIS corresponding to a particle of mass 5*M[PLANCK]. A couple of points are to be noted, (1) the vertices can be names, numbers, or anything one want as long as it is a valid argument to a hashed array, (2) the factor of 5 in the mass is pointless, as numerical factors are ignored in asymptotic bounds. The third argument, RHUBARB, is a name for the line, which is solely there for debugging purposes: internally ASYMP will invent its own name for the line. This argument is, like rhubarb, best left by the side of the plate and ignored. EXTLINE describes an external leg to our Feynman diagram. EXTLINE(ROME, CELERY,- P+2*Q) says that there is an external leg attached to our graph at vertex ROME carrying momentum 2*Q-P into the graph. It is one's own responsibility to ensure over all momentum conservation. The second argument, CELERY, has great similarities to RHUBARB and is also best forgotten (well, it has to be there, but it seems to serve no other useful role in life). OK, so we have now defined our graph. DIAGRAM sets up tables of lines containing their masses etc., assigns internal loop-momenta, routes all momenta through the graph, and tells one the number of loops in the diagram. If we had been nosey and typed a ; rather than a $ at DIAGRAM, it would have returned a list of the form [G000002653,G000005532,G000007771]. The Goo's are internal line-names of no interest to one, other than that they are used by later programs to index the tables set up by DIAGRAM and its cohorts. The only point of interest is that GRAPH1 is now a list of variable names, in other words it behaves just like any other MACSYMA variable, which is not too surprising because it IS just like any other MACSYMA variable. In line (C2) we get down to the real business of the day. We use the function BOUND to find the asymptotic behavior of GRAPH1 when the Euclidean momenta p and q and the mass m are much smaller than the mass mm, and both are much smaller than INF (of course: the need to put in INF by hand is just a foible 3 Preliminary Version The ASYMP package June 1, 1982 of the program, so don't forget it!). To put it another way, we set up three mass scales, which we shall call m, mm, and INF, such that any mass of order m is asymptotically bounded by (or, in everyday terms, much less than) any mass of order mm, and in turn mm << INF. The second argument to BOUND, therefore, is a list of mass-scales, each of which is either a mass/momentum or a list of masses and/or momenta of the same scale. The result of BOUND is that GRAPH1 is bounded by (an implicit constant) times log(mm/m)/mm^2, at least for mm/m large enough. For further examples look at the files SHARE1;ASYMP DEMOUT, SHARE1;ASYMP DEMOU1, etc., which are the output from the demo files SHARE1;ASYMP DEMO, SHARE1;ASYMP DEMO1, etc. III. Method. For a long write up, see the paper "The Asymptotic Behavior of Feynman Integrals," Maryland Physics publication #PP-81-188. IV. Syntax. As the syntax has been described in section II, we just summarize it below: DIAGRAM(<pseudo-function>,<pseudo-function>,...); LINE(<from-vertex>,<to-vertex>,<name>,<mass>); EXTLINE(<to-vertex>,<name>,<momentum>); 4 Preliminary Version The ASYMP package June 1, 1982 BOUND(<diagram>,[<mass-scale>,<mass-scale>,...,INF]); <mass-scale> :: mass | momentum | [<mass-scale>,<mass- scale>,...] V. Notes. For further information please send mail to ADK@... or WEC@... 5 --- NEW FILE: asympa.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/asympa.usg was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: dblint.dem --- /* Example of use of double integrator (DBLINT) */ /* F must be a function of two variables, R and S must be functions of one variable, and A and B must be floating point numbers (all of the functions must be translated or compiled prior to using DBLINT) */ /* Get the fasl file containing DBLINT(F,R,S,A,B) */ LOAD('DBLINT); /* To get the double integral of exp(x-y^2) over the region bounded by y=1 and y=2+x and x=0 to x=1 */ /* Define the integrand as a function of two variables */ F(X,Y):=(MODE_DECLARE([X,Y],FLOAT),EXP(X-Y^2)); /* Define the lower and upper limits on the inner (y in this case) integral as a function of the outer variable (x in this case) */ R(X):=(MODE_DECLARE(X,FLOAT),1.0); S(X):=(MODE_DECLARE(X,FLOAT),2.0+X); /* Now translate these functions for the sake of efficiency */ TRANSLATE(F,R,S); /* Call the DBLINT function with quoted arguments for function names, and floating point values for the endpoints of the outer (x) integration */ DBLINT_ANS:DBLINT('F,'R,'S,0.0,1.0); /* Compare with the exact answer */ INTY:RISCH(EXP(X-Y^2),Y); XINT:EV(INTY,Y:2+X)-EV(INTY,Y:1); DINT:RISCH(XINT,X); ANS:EV(DINT,X:1.0)-EV(DINT,X:0.0),NUMER; /* Relative error check here */ (ANS-DBLINT_ANS)/ANS; /* Quite reasonable */ /* Of course, the DBLINT routine will still work even when we choose an area over which the closed-form integral fails to be expressible in terms of standard transcendental functions */ S1(X):=(MODE_DECLARE(X,FLOAT),2.0+X^(3/2)); TRANSLATE(S1); DBLINT('F,'R,'S1,0.0,1.0); /* To see what happens when this integral is attempted in closed-form, you may BATCH(DBLINT,DEMO1,DSK,SHARE1), and that will also show sample use of FLOATDEFUNK and QUANC8 */ "end of dblint.dem"$ --- NEW FILE: dblint.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/dblint.dm1 was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: dblint.usg --- The file SHARE1;DBLINT FASL contains a double-integral routine which was written in top-level macsyma and then translated and compiled to machine code. It uses the Simpson's Rule method (see note below) in both the x and y directions to calculate /B /S(X) | | | | F(X,Y) DY DX . | | /A /R(X) The function F(X,Y) must be a translated or compiled function of two variables, and R(X) and S(X) must each be a translated or compiled function of one variable, while A and B must be floating point numbers. The routine has two global variables which determine the number of divisions of the x and y intervals: DBLINT_X and DBLINT_Y, both of which are initially 10, and can be changed independently to other integer values [there are 2*DBLINT_X+1 points computed in the x direction, and 2*DBLINT_Y+1 in the y direction]. The routine subdivides the x axis and then for each value of x it first computes r(x) and s(x); then the y axis between r(x) and s(x) is subdivided and the integral along the y axis is performed using Simpson's Rule; then the integral along the x axis is done using Simpson's Rule with the function values being the y-integrals. This procedure may be numerically unstable for a great variety of reasons, but is reasonably fast: avoid using it on highly oscillatory functions and functions with singularities (poles or branch points in the region). The y integrals depend on how far apart r(x) and s(x) are, so if the distance s(x)-r(x) varies rapidly with x, there may be substantial errors arising from truncation with different step-sizes in the various y integrals. One can increase dblint_x and dblint_y in an effort to improve the coverage of the region, at the expense of computation time. The function values are not saved, so if the function is very time-consuming, you will have to wait for re-computation if you change anything (sorry). It is required that the functions F, R, and S be either translated or compiled prior to calling DBLINT. This will result in orders of magnitude speed improvement over interpreted code in many cases! Ask LPH (or GJC) about using these numerical aids. The file SHARE1;DBLINT DEMO can be run in batch or demo mode to illustrate the usage on a sample problem; the file SHARE1;DBLNT DEMO1 is an extension of the DEMO which also makes use of other numerical aids, FLOATDEFUNK and QUANC8. Please send all bug notes and questions to LPH@... Note: Simpson's Rule specifies that /X[2*N] | | F(X) DX = H/3* (F(X[0]) + | /X[0] 4*(F(X[1])+F(X[3])+...+F(X[2*N-1])) + 2*(F(X[2])+F(X[4])+...+F(X[2*N-2])) + F(X[2*N])) in one dimension, where H is the distance between the equally spaced X[N]'s, and DBLINT_X=N. The error in this formulation is of order H^5*N*DIFF(F(X),X,4) for some X in (X[0],X[2*N]). --- NEW FILE: dimen.dem --- /* First load the necessary file: */ LOAD("dimen.mc")$ /* It is conjectured that for thermistors there is a physical relationship between the voltage drop, current, ambient temperature, room-temperature resistance, convective heat transfer coefficient, and a constant, BETA, having the dimension of temperature. First, to see if the dimension of BETA is already known: */ GET(BETA, 'DIMENSION); /* It is not. To establish it: */ DIMENSION(BETA=TEMPERATURE); /* To automatically determine a set of dimensionless variables sufficient to characterize the physical relation: */ NONDIMENSIONALIZE([VOLTAGE, CURRENT, TEMPERATURE, RESISTANCE, HEATTRANSFERCOEFFICIENT, BETA]); /* We learn that the relation may be expressed as a function of only the above 3 variables rather than a function of the six physical quantities. Evidently dimensions were preestablished for all but the last of these particular input quantities, but an appropriate error message would have informed us if this were not so. An extensive set of dimensions have been prestablished, as may be seen from the listing of DIMEN >. As another example, there is thought to be a relation between the viscosity, average velocity, molecular mass, and repulsion coefficient of a gas. The repulsive force between two molecules is believed to be of the form K/DISTANCE^N, with unknown N, so K must have the following dimensions: */ DIMENSION(K=MASS*LENGTH^(N+1)/TIME^2) $ /* To get the computation time in milliseconds to be printed automatically: */ CPUTIME: TRUE $ /* To do a dimensional analysis of the gas viscosity problem: */ NONDIMENSIONALIZE([VISCOSITY, K, MASS, VELOCITY]); /* The physical relation must be expressible as a function of this one dimensionless variable, or equivalently, this variable must equal a constant. Consequently, physical measurements may be used to determine N. It turns out to be in the range 7 to 12 for common gases. As a final example, suppose that we conjecture a relation between the deflection angle of a light ray, the mass of a point mass, the speed of light, and the distance from the mass to the point of closest approach: */ NONDIMENSIONALIZE([ANGLE, MASS, LENGTH, SPEEDOFLIGHT]); /* We learn that there cannot be a dimensionless relation connecting all of these quantities and no others. Let us also try including the constant that enters the inverse-square law of gravitation: */ NONDIMENSIONALIZE([ANGLE, MASS, LENGTH, SPEEDOFLIGHT, GRAVITYCONSTANT]); /* Altermatively, for astrophysics problems such as this,we may prefer to use a dimensional basis in which the gravity constant is taken as a pure number, eliminating one member from our dimensional basis: */ %PURE: CONS(GRAVITYCONSTANT, %PURE); /* Note that the latter two of the above constants are pure numbers by default, respectively eliminating TEMPERATURE and CHARGE from the basis, but the user may include all five of TEMPERATURE, CHARGE, MASS, LENGTH, and TIME in the basis by resetting %PURE to []. Altermatively, the user may wish to include SPEEDOFLIGHT in %PURE for relativistic problems or PLANCKSCONSTANT for quantum problems. For dimensional analysis it doesn't really matter which basis member is eliminated by each pure constant, but in fact the latter two respectively eliminate LENGTH and TIME, whereas GRAVITYCONSTANT eliminates MASS. To proceed with our analysis: */ NONDIMENSIONALIZE([ANGLE, MASS, LENGTH, SPEEDOFLIGHT]); --- NEW FILE: dimen.usg --- File DIMEN > contains functions for automatic dimensional analysis, and file DIMEN DEMO contains a demonstration. Usage is of the form NONDIMENSIONALIZE(list of physical quantities); The returned value is a sufficient list of nondimensional products of powers of the physical quantities. A physical relation between only the given physical quantities must be expressible as a relation between the nondimensional quantities. There are usually fewer nondimensional than physical quantities, which reduces the number of experiments or numerical computations necessary to establish the physical relation to a specified resolution, in comparison with the number if all but one dependent physical variable were independently varied. Also, the absence of any given physical quantity in the output reveals that either the quantity is irrelevant or others are necessary to describe the relation. The program already knows an extensive number of relations between physical quantities, such as VELOCITY=LENGTH/TIME. (CPUTIME plays the role of the customary MACSYMA global variable TIME.) The user may over-ride or supplement the prespecified relations by typing DIMENSION(equation or list of equations); where each equation is of the form indeterminate=expression, where expression is a product or quotient of powers of none or more of the indeterminates CHARGE, TEMPERATURE, LENGTH, TIME, or MASS. To see if a relation is already established type GET(indeterminate, 'DIMENSION); The result of NONDIMENSIONALIZE usually depends upon the value of the global variable %PURE, which is set to a list of none or more of the indeterminates ELECTRICPERMITTIVITYOFAVACUUM, BOLTZMANNSCONSTANT, SPEEDOFLIGHT, PLANCKSCONSTANT, GRAVITYCONSTANT, corresponding to the relation between charge and force, temperature and energy, length and time, length and momentum, and the inverse-square law of gravitation respectively. Each included relation is used to eliminate one of CHARGE, TEMPERATURE, LENGTH, TIME, or MASS from the dimensional basis. To avoid omission of a possibly relevant nondimensional grouping, either include the relevant constant in %PURE or in the argument of NONDIMENSIONALIZE if the corresponding physical effect is thought to be relevant to the problem. However, the inclusion of unnecessary constants, especially the latter three, tends to produce irrelevant or misleading dimensionless groupings, defeating the purpose of dimensional analysis. As an extreme example, if all five constants are included in %PURE, all physical quantities are already dimensionless. %PURE is initially set to '[ELECTRICPERMITTIVITYOFVACUUM, BOLTZMANNSCONSTANT], which is best for most engineering work. %PURE must not include any of the other 3 constants without also including these 2. Send problems and suggestions to STOUTE. REFERENCES: R. Kurth, "Dimensional Analysis and Group Theory in Astrophysics", Perggamon Press. H.L. Langhaar, "Dimensional Analysis and Theory of Models", John Wiley and Sons. D.R. Stoutemyer, "Automatic Dimensional Analysis, Using Computer Symbolic Mathematics", report, Electrical Engineering Department, University of Hawaii, Honolulu, Hawaii 96822. --- NEW FILE: facexp.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/facexp.dem was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: facexp.usg --- File: SHARE1;FACEXP USAGE The file SHARE1;FACEXP FASL contains several related functions that provide the user with the ability to structure expressions by controlled expansion. This capability is especially useful when the expression contains variables that have physical meaning, because it is often true that the most economical form of such an expression can be obtained by fully expanding the expression with respect to those variables, and then factoring their coefficients. While it is true that this procedure is not difficult to carry out using standard Macsyma functions, additional fine-tuning may also be desirable, and these finishing touches can be more difficult to apply. The function FACSUM and its related forms provide a convenient means for controlling the structure of expressions in this way. Another function, COLLECTTERMS, can be used to add two or more expressions that have already been simplified to this form, without resimplifying the whole expression again. This function is particularly useful when the expressions are large and address space or cpu time is in short supply. For a complete demonstration of the use of these functions, BATCH(FACEXP,DEMO,DSK,SHARE1)$ Brief descriptions of these functions and their switches follow. FACSUM(exp,arg1,arg2,...) returns a form of exp which depends on the argi. The argi can be any form suitable for RATVARS, or they can be lists of such forms. If the argi are not lists, then the form returned will be fully expanded with respect to the argi, and the coefficients of the argi will be factored. These coefficients will be free of the argi, except perhaps in a non-rational sense. If any of the argi are lists, then all such lists will be combined into a single list, and instead of calling FACTOR on the coefficients of the argi, FACSUM will call itself on these coefficients, using this newly constructed single list as the new argi for this recursive call. This process can be repeated to arbitrary depth by nesting the desired elements in lists. It is possible that one may wish to FACSUM with respect to more complicated subexpressions, such as LOG(X+Y). Such arguments are also permissible. With no variable specification, for example FACSUM(exp), the result returned is the same as that returned by RATSIMP(exp). Occasionally the user may wish to obtain any of the above forms for expressions which are specified only by their leading operators. For example, one may wish to FACSUM with respect to all LOG's. In this situation, one may include among the argi either the specific LOG's which are to be treated in this way, or alternatively, either the expression OPERATOR(LOG) or 'OPERATOR(LOG). If one wished to FACSUM the expression EXP with respect to the operators OP1, OP2, .. ..., OPn, one would evaluate FACSUM(EXP, OPERATOR(OP1, OP2, ...OPn)). The OPERATOR form may also appear inside list arguments. In addition, the setting of the switches FACSUM_COMBINE and NEXTLAYERFACTOR may affect the result of FACSUM as follows: NEXTLAYERFACTOR[FALSE] if TRUE will force the recursive calls of FACSUM to be applied to the factors of the factored form of the coefficients of the argi. If FALSE then FACSUM will be applied to each coefficient as a whole whenever recusive calls to FACSUM occur as described above. In addition, inclusion of the atom NEXTLAYERFACTOR in the argument list of FACSUM has the effect of NEXTLAYERFACTOR:TRUE, but for the next level of the expression ONLY. Since NEXTLAYERFACTOR is always bound to either TRUE or FALSE, it must be presented single-quoted whenever it is used in this way. FACSUM_COMBINE[TRUE] controls the form of the final result returned by FACSUM when its argument is a quotient of polynomials. If FACSUM_COMBINE is FALSE then the form will be returned as a fully expanded sum as described above, but if TRUE, then the formed returned is a ratio of polynomials, with each polynomial in the form described above. The TRUE setting of this switch is useful when one wants to FACSUM both the numerator and denominator of a rational expression, but does not want the denominator to be multiplied through the terms of the numerator. FACTORFACSUM(exp, arg1, arg2, ...argN) returns a form of exp which is obtained by calling FACSUM on the factors of exp with the argI as arguments. If any of the factors of exp is raised to a power, both the factor and the exponent will be processed in this way. COLLECTTERMS(arg1, arg2, ... argn) If several expressions have been simplified with FACSUM, FACTORFACSUM, FACTENEXPAND, FACEXPTEN or FACTORFACEXPTEN, and they are to be added together, it may be desirable to combine them using the function COLLECTERMS. COLLECTERMS can take as arguments all of the arguments that can be given to these other associated functions with the exception of NEXTLAYERFACTOR, which has no effect on COLLECTTERMS. The advantage of COLLECTTERMS is that it returns a form similar to FACSUM, but since it is adding forms that have already been processed by FACSUM, it does not need to repeat that effort. This capability is especially useful when the expressions to be summed are very large. Note: FACTENEXPAND, FACEXPTEN, and FACTORFACEXPTEN are available only after loading DIAGEVAL. They are special functions used for tensor manipulation. See the file ASB;MANUAL DOC. --- NEW FILE: hypgeo.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/hypgeo.dem was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: hypgeo.usg --- The Hypergeometric Special Functions Package HYPGEO is still under development. At the moment it will find the Laplace Transform or rather, the integral from 0 to INF of some special functions or combinations of them. The factor, EXP(-P*var) must be explicitly stated. The syntax is as follows: SPECINT(EXP(-P*var)*expr,var); where var is the variable of integration and expr may be any expression containing special functions (at your own risk). Special function notation follows: %J[index](expr) Bessel Funct 1st Kind %K[index](expr) " " 2nd Kind %I[ ]( ) Modified Bessel %HE[ ]( ) Hermite Poly %P[ ]( ) Legendre Funct %Q[ ]( ) Legendre of second kind HSTRUVE[ ]( ) Struve H Function LSTRUVE[ ]( ) " L Function %F[ ]([],[],expr) Hypergeometric Function GAMMA() GAMMAGREEK() GAMMAINCOMPLETE() SLOMMEL %M[]() Whittaker Funct 1st Kind %W[]() " " 2nd " For a better feeling for what it can do, do DEMO(HYPGEO,DEMO,SHARE1); . --- NEW FILE: inteqn.usg --- ****INTEGRAL EQUATION SOLVER*** A VERSION OF PROF. STOUTEMYER'S (STOUTE) INTEGRAL EQUATION SOLVER CAN BE OBTAINED BY DOING "BATCH(INTEQN,>,DSK,SHARE1)". IN ORDER TO FREE SOME STORAGE IT DOES A KILL(LABELS) AT THE END SO BE FORWARNED TO GIVE NAMES TO ANY EXPRESSIONS YOU WANT TO KEEP AROUND BEFORE BATCHING IN THE INTEGRAL EQUATION SOLVER. THE USAGE IS VERY SIMPLE. THE MAIN FUNCTION IS CALLED "IEQN" AND IT TAKES FIVE ARGUMENTS THOUGH ONLY THE FIRST TWO ARE REQUIRED. IF TRAILING ARGUMENTS ARE OMITTED THEY WILL DEFAULT TO PRESET VALUES WHICH WILL BE ANNOUNCED. TWO TYPES OF EQUATIONS ARE CONSIDERED. A SECOND-KIND AND A FIRST-KIND GIVEN BY D3 AND D4 BELOW, RESPECTIVELY. THE UNKNOWN FUNCTION IN THESE EQUATIONS IS P(X) WHILE Q,W,A, & B ARE GIVEN FUNCTIONS OF THE INDEPENDENT VARIABLE. (C3) P(X)=Q(X,P(X),'INTEGRATE(W(X,U,P(X),P(U)),U,A(X),B(X))); B(X) / [ (D3) P(X) = Q(X, P(X), I W(X, U, P(X), P(U)) dU) ] / A(X) (C4) 'INTEGRATE(W(X,U,P(U)),U,A(X),B(X))=F(X); B(X) / [ (D4) I W(X, U, P(U)) dU = F(X) ] / A(X) ALTHOUGH THESE ARE THE GENERAL FORMS, MOST OF THE SOLUTION TECHNIQUES REQUIRE PARTICULAR FORMS OF Q AND W. THE TECHNIQUES USED ARE: FOR SECONDKIND EQUATIONS -- FLFRNK2ND : FOR FIXED-LIMIT, FINITE-RANK INTEGRANDS. VLFRNK : FOR VARIABLE-LIMIT, FINITE-RANK INTEGRANDS. TRANSFORM : LAPLACE TRANSFORM FOR CONVOLUTION TYPES. FREDSERIES : FREDHOLM-CARLEMAN SERIES FOR LINEAR EQUATIONS. TAILOR : TAYLOR SERIES FOR QUASI-LINEAR VARIABLE-LIMIT EQUATIONS. NEUMANN : NEUMANN SERIES FOR QUASI-SECOND KIND EQUATIONS. COLLOCATE : COLLOCATION USING A POWER SERIES FORM FOR P(X) EVALUATED AT EQUALLY SPACED POINTS. FOR FIRSTKIND EQUATIONS -- FLFRNK1ST : FOR FIXED-LIMIT, FINITE-RANK INTEGRANDS. VLFRNK : FOR VARIABLE-LIMIT, FINITE-RANK INTEGRANDS. ABEL : FOR SINGULAR INTEGRANDS TRANSFORM :SEE ABOVE COLLOCATE :SEE ABOVE FIRSTKINDSERIES : ITERATION TECHNIQUE SIMILAR TO NEUMANN SERIES. ---------------------------------------- THE CALLING SEQUENCE IS ... IEQN(EQNS,UNKS,TECH,NAPPROX,GUESSES) WHERE: EQNS - THE INTEGRAL EQUATION OR A LIST OF THEM. UNKS - THE UNKNOWN FUNCTION OR A LIST OF THEM. TECH - THE TECHNIQUE TO BE TRIED SELECTED FROM THOSE GIVEN ABOVE, OR "FIRST" MEANING TO TRY THE FIRST TECHNIQUE WHICH FINDS A SOLUTION, OR "ALL" WHICH MEANS TO TRY ALL APPLICABLE TECHNIQUES. NAPPROX - THE MAXIMUM NUMBER OF TERMS TO TAKE FOR TAYLOR, NEUMANN, FIRSTKINDSERIES, OR FREDSERIES, OR THE NUMBER OF PARAMETERS TO USE IN COLLOCATION. GUESSES - THE INITIAL GUESS FOR NEUMANN OR FIRSTKINDSERIES OR A LIST OF THEM. DEFAULT VALUES FOR THE 3RD THRU 5TH PARAMETERS ARE: TECH - FIRST NAPPROX - 1 GUESSES - NONE. WHICH WILL CAUSE FIRSTKINDSERIES TO USE F(X) AND NEUMANN SERIES TO USE Q(X,0,0) AS INITIAL GUESS. ### THE VALUE RETURNED BY IEQN IS A LIST OF LABELS OF SOLUTION LISTS. THE SOLUTION LISTS ARE PRINTED AS THEY ARE FOUND UNLESS THE OPTION VARIABLE "IEQNPRINT" IS SET TO FALSE. THESE LISTS ARE OF THE FORM [SOLUTION, TECHNIQUE USED, NTERMS, FLAG] WHERE FLAG IS ABSENT IF THE SOLUTION IS EXACT. OTHERWISE IT IS THE WORD "APPROXIMATE" OR "INCOMPLETE" CORRESPONDING TO AN INEXACT OR NON-CLOSED FORM SOLUTION RESPECTIVELY. IF A SERIES METHOD WAS USED, NTERMS GIVES THE NUMBER OF TERMS TAKEN WHICH COULD BE LESS THAN THE N GIVEN TO IEQN IF AN ERROR WAS ENCOUNTERED PREVENTING GENERATION OF FURTHER TERMS. IF FLAG ****** FOR EXAMPLES DO BATCH(INTEXS,">",DSK,SHARE1) WHICH WILL LOAD AN ARRAY CALLED "EQ" WITH ABOUT 43 SAMPLE INTEGRAL EQUATIONS. THEN TRY IEQN(EQ[1],P(X)), IEQN(EQ[30],P(X),ALL), ETC. --- NEW FILE: intsce.usg --- INTSCE LISP contains a routine INTSCE(EXPR,VAR) which integrates EXPR w.r.t. VAR where EXPR is of the form: EXP(A*X+B)*COS(C*X)^N*SIN(C*X)^M EXPR may be any expression, but if it is not in the above form then the regular integration program will be invoked if the switch ERRINTSCE [FALSE] is FALSE. If it is TRUE then INTSCE will err out. The LISP file contains a main program $INTSCE which does the pattern matching and error checking and a subroutine $SCEINT which does the actual integration. The MACSYMA form of this routine can be found in SHARE;INTSCE > . Richard Bogen 7/22/74 --- NEW FILE: invert.usg --- BMT@... 09/29/77 21:32:25 The file SHARE;INVERT > finds the inverse of a matrix using the adjoint method. This allows a user to compute the inverse of a matrix with bfloat entries or polynomials with floating pt. coefficients without converting to cre-form. The DETERMINANT command is used to compute cofactors, so if RATMX is FALSE (the default) the inverse is computed without changing the representation of the elements. The functions ADJOINT and INVERT are provided. The current (temporary?) implementation is inefficient for matrices of high order. The DETOUT flag if true keeps the determinant factored out of the inverse. Note: the results are not automatically expanded. If the matrix originally had polynomial entries, better appearing output can be generated by EXPAND(INVERT(mat)),DETOUT. If it is desirable to then divide through by the determinant this can be accomplished by XTHRU(%) or alternatively from scratch by EXPAND(ADJOINT(mat))/EXPAND(DETERMINANT(mat)). INVERT(mat):=ADJOINT(mat)/DETERMINANT(mat). --- NEW FILE: kach.dem --- /* Demo for Purtilo's implementation of Hacijan's algorithm.... */ load("kach")$ load("haipart")$ fpprec:16$ If Status(Feature,ITS)=True Then Translate(hach,check,iterate,fx,fq,fh,getl,criterion,accelerate)$ /* dynamalloc:true$ */ /* now get a matrix a and a matrix b for the problem of solving a.x < b ... */ a:bfloat(matrix([1,0],[0,1],[-1,-1]))$ b:bfloat(matrix([-1],[-1],[4.00000001]))$ n:2$ m:3$ /* the user can find the L as in Hacijan's abstract by doing ... */ getl(a,b,n,m); /* and can find a solution vector x for x[1] < -1 x[2] < -1 x[1] + x[2] >= -4 by doing */ show%all:true$ hach(a,b,n,m,1); /* Note that this last example does not pay attention to the L of the theory. Doing so, taking a bit longer, is */ show%all:false$ hach(a,b,n,m,9.58); --- NEW FILE: nchrpl.dem --- /* First load the necessary file: */ LOAD('NCHRPL); /* Now generate a 6 x 6 matrix to work with: */ H[I,J]:=(I+J)/(I*J); TEST_MATRIX:60*GENMATRIX(H,6,6,1,1); /* We wish to compare the performance of NCHARPOLY and CHARPOLY. */ SHOWTIME:ALL$ NCHARPOLY(TEST_MATRIX,X); RATSIMP(CHARPOLY(TEST_MATRIX,X)); /* The RATSIMP was necessary to obtain a form as neat as the one returned by NCHARPOLY. Another function in this package is MATTRACE, which takes the trace of a matrix: */ MATTRACE(TEST_MATRIX); --- NEW FILE: qq.dem --- /* Numerical Integration Demo for QUANC8*/ /* SETUP_AUTOLOAD(QQ,QUANC8); */ /* show the time and the garbage collection time, and use the GC-demon to make better use of flonum space */ SHOWTIME:'ALL; IF STATUS(FEATURE,ITS)=TRUE THEN LOAD("GCDEMN"); /* Sample highly oscillatory problem */ G(X):=(MODE_DECLARE(X,FLOAT),SIN(1.0/X)); /* Inefficient way to use QUANC8 */ QUANC8(G(x),x,0.01,0.1); /* give G as the thing to compile, then a semi-colon to read it in */ compile(); /* see how much faster now, using fast version of QUANC8 */ QUANC8('G,0.01,0.1); /* note that romberg doesn't easily manage oscillatory behaviors */ ERRCATCH(ROMBERG('G,0.01,0.1)); /* a function with two large and narrow peaks, so hard to integrate with a fixed step-size */ H(X):=(MODE_DECLARE(X,FLOAT),1.0/(.0001+(X-.3)^2)+1.0/(.0025+(X-.9)^2)); /* give H and a semi-colon in response */ compile(); /* the exact answer for integrate(h(x),x) is 100.0*atan(100.0*x-30.0)+20.0*atan(20.0*x-18.0) */ /* the exact answer for integrate(h(x),x,0.0,1.0) is 361.847622. note that the romberg result is more accurate in this case */ ROMBERG('H,0.0,1.0); /* but at a large time cost compared to adaptive */ QUANC8('H,0.0,1.0); /* reduce QUANC8's relative error tolerance by a factor of 10 */ QUANC8_RELERR:1.0E-5; /* does not take much longer, and gives much better result now */ QUANC8('H,0.0,1.0); /* put QUANC8_RELERR back to 1.0e-4 */ QUANC8_RELERR:1.0E-4$ /* the exact answer for integrate(h(x),x,0.0,10.0) is 372.33606. while too tough for fixed step size method */ ERRCATCH(ROMBERG('H,0.0,10.0)); /*QUANC8 is still super speedy even now */ QUANC8('H,0.0,10.0); /* double integral of sample function exp(x)*sin(y) from x=0.0 to 2.0, y=0.0 to 1.0 */ /* you must be sure that the quanc8 source is loaded (load("qq");) before you translate a call to quanc8 in a function, else an error in macro expansion will occur when it can't figure out if the 3 or 4 arg version is used */ /* the exact answer for integrate(integrate(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0) is 2.93703434. integrate over y to get the final answer */ QUANC8(quanc8(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0); p():=QUANC8(lambda([y],mode_declare(y,float), quanc8(lambda([x],mode_declare(x,float), exp(x)*sin(y)), 0.0,2.0)),0.0,1.0); p(); translate(p); p(); /* give p and a semi-colon */ compile(); p(); --- NEW FILE: qq.usg --- The file SHARE1;QQ FASL contains a function QUANC8 which can take either 3 or 4 arguments. The 3 arg version computes the integral of the function specified as the first argument over the interval from lo to hi as in QUANC8('function name,lo,hi). The function name should be quoted. The 4 arg version will compute the integral of the function or expression (first arg) with respect to the variable (second arg) over the interval from lo to hi as in QUANC8(<f(x) or expression in x>,x,lo,hi). The method used is the Newton-Cotes 8th order polynomial quadrature, and the routine is adaptive. It will thus spend time dividing the interval only when necessary to achieve the error conditions specified by the global variables QUANC8_RELERR (default value=1.0e-4) and QUANC8_ABSERR (default value=1.0e-8) which give the relative error test |integral(function)-computed value|< quanc8_relerr*|integral(function)| and the absolute error test |integral(function)-computed value|<quanc8_abserr. The error from each subinterval is estimated and the contribution from a subinterval is accepted only when the integral over the subinterval satisfies the error test over the subinterval. The total estimated error of the integral is contained in the global variable QUANC8_ERREST (default value=0.0). The global variable QUANC8_FLAG (default value=0.0) will contain valuable information if the computation fails to satisfy the error conditions. The integer part will tell you how many subintervals failed to converge and the fractional part will tell you where the singular behavior is, as follows: singular point=lo+(1.-frac part)*(hi-lo). Thus quanc8(tan(x),x,1.57,1.6); gives frac=.97 so trouble is at 1.57+.03*.03= 1.5709 (=half pi). If QUANC8_FLAG is not 0.0, you should be cautious in using the return value, and should try ROMBERG or a Simpson method and see if the result checks. Analysis of possible singular behavior might be advisable. You may get QUANC8_FLAG=<integer>.0 and an error message (such as division by 0) when a singular point is hit in the interval. You will have to find the singularity and eliminate it before QUANC8 will get an answer. Functions which have very large derivatives may throw the error estimate way off and cause the wrong points to be used, and a wrong answer returned. Try romberg(exp(-.002*x^2)*cos(x)^2,x,0.,100.); with the default tolerance, and quanc8(exp(-.002*x^2)*cos(x)^2,x,0.,100.); with quanc8_relerr=1.e-7 and 1.e-8. The last result is consistent with romberg while the previous one is off by a factor of 2! This is due to the bad behavior of the derivatives near x=10.0 which cause the adaptive routine to have trouble. If you use quanc8('f,a,c)+quanc8('f,c,b) where a<c<b, you will do better in such cases. Typing a control-right-bracket () will give a printout of where the computation is at the moment, and how much time has been used. You can do batch("qq demo"); for some comparisons with the ROMBERG numerical integrator (which is not adaptive). Note that ROMBERG usually gives more accurate answers for comparable tolerances, while QUANC8 will get the same answer faster even with a smaller tolerance, because ROMBERG subdivides the whole interval if the total result is not within error tolerance, while QUANC8 improves only where needed, thus saving many function calls. ROMBERG will also fail to converge when oscillatory behavior is overwhelming, while QUANC8 will adapt in the regions as it sees fit. (The global variable ROMBERGMIN is designed to allow you a minimum number of function calls in such cases, so that exp(-x)*sin(12*x) can be integrated from 0 to 4*%pi without erroneously giving 0. from the first few function calls.) To make your macsyma user function callable in the fastest way, you must use MODE_DECLARE and then translate and compile the function. Read TRANSL;TRANSL NEWS for info on this, or ask me for more info. The speed of the computation may be increased by well over an order of magnitude when compilation is used. If you do multiple integrals, it is really necessary to compile the function in order to avoid the time spent on function calls. A sample use of QUANC8 for a double integral is in the batch("qq demo"); and compilation is nearly a hundred times faster in doing the work! Comments, bugs, etc. should be sent to LPH@... --- NEW FILE: rncomb.dem --- load("rncomb")$ /* COMBINE will not combine the terms in EXP1 because their denominators are not identical: */ EXP1:X/(2*(X+Y))+Z/(X+Y); COMBINE(EXP1); /* RNCOMBINE will: */ RNCOMBINE(%); /* RNCOMBINE will combine terms whenever their denominators are identical or differ by a numerical factor. */ EXP2:A/2+B/2+C/3+D/3; COMBINE(exp2); RNCOMBINE(exp2); --- NEW FILE: rncomb.usg --- RNCOMBINE(exp) transforms exp by combining all terms of exp that have identical denominators or denominators that differ from each other by numerical factors only. This is slightly different from the behavior of COMBINE, which collects terms that have identical denominators. Setting PFEFORMAT:TRUE and using COMBINE will achieve results similar to those that can be obtained with RNCOMBINE, but RNCOMBINE takes the additional step of cross-multiplying numerical denominator factors. This results in neater forms, and the possiblity of recognizing some cancellations. Bugs to ASB. --- NEW FILE: sets.html --- <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> <HTML><HEAD> <TITLE>404 Not Found</TITLE> </HEAD><BODY> <H1>Not Found</H1> The requested URL /crew/mike/maxima/html/macref/share1/sets.dmo was not found on this server.<P> <HR> <ADDRESS>Apache/1.3.27 Server at starship.python.net Port 80</ADDRESS> </BODY></HTML> --- NEW FILE: sets.usg --- A Fast Sets Package due to GJC ------------------------------ There is a fast sets package available by doing LOAD(SETS); . The set constructor is the "{". So X:{A,B,C,D,E}; creates a set. The usual primitives UNION, INTERSECTION, SETDIFF, SYMDIFF are defined. (SYMDIFF(A,B)=UNION(SETDIFF(A,B),SETDIFF(B,A)) = SETDIFF(UNION(A,B),INTERSECTION(A,B)) ). Predicates are: ELEMENTP(X,SET), EMPTYP(SET), and SUBSETP(SET1,SET2) CARDINAL(SET) returns the cardinality. There are two mapping-like functions which are provided for sets. PREDSET(PREDICATE,SET) returns the set of all elements of SET such that the PREDICATE returns TRUE. e.g. X:{1,2,3,4,5,6,7,8,9,10,11} PREDSET(LAMBDA([U],IS(ABS(U-6)<3)),X) returns {5,6,7} MAPSET(FUNCTION,SET) creates a set from the results of applying the FUNCTION to the elements of the SET. ELEMENTS(SET) returns a list of the elements. The sets are not represented as lists. The set-algebraic functions (UNION, INTERSECTION, SETDIFF, SYMDIFF, PREDSET, CARDINAL), all operate on the internal representation of sets and as such are fast. Things which have to be converted from the set representation to non-set are a bit slower, the things which make sets from raw elements are slower still, however, they are somewhat faster than CONS on the average. If there is interest, "INFINITE" sets could be implemented. There, of course, one has to be careful in converting between representations, obviously one can't print a list of the elements of an infinite set. -gjc |

Powered by

Apache Allura™