This is Info file calc.info, produced by Makeinfo-1.55 from the input file calc.texinfo. This file documents Calc, the GNU Emacs calculator. Copyright (C) 1990, 1991 Free Software Foundation, Inc. Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies. Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided also that the section entitled "GNU General Public License" is included exactly as in the original, and provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one. Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions, except that the section entitled "GNU General Public License" may be included in a translation approved by the author instead of in the original English. File: calc.info, Node: Customizing the Integrator, Next: Numerical Integration, Prev: Integration, Up: Calculus Customizing the Integrator -------------------------- Calc has two built-in rewrite rules called `IntegRules' and `IntegAfterRules' which you can edit to define new integration methods. *Note Rewrite Rules::. At each step of the integration process, Calc wraps the current integrand in a call to the fictitious function `integtry(EXPR,VAR)', where EXPR is the integrand and VAR is the integration variable. If your rules rewrite this to be a plain formula (not a call to `integtry'), then Calc will use this formula as the integral of EXPR. For example, the rule `integtry(mysin(x),x) := -mycos(x)' would define a rule to integrate a function `mysin' that acts like the sine function. Then, putting `4 mysin(2y+1)' on the stack and typing `a i y' will produce the integral `-2 mycos(2y+1)'. Note that Calc has automatically made various transformations on the integral to allow it to use your rule; integral tables generally give rules for `mysin(a x + b)', but you don't need to use this much generality in your `IntegRules'. As a more serious example, the expression `exp(x)/x' cannot be integrated in terms of the standard functions, so the "exponential integral" function `Ei(x)' was invented to describe it. We can get Calc to do this integral in terms of a made-up `Ei' function by adding the rule `[integtry(exp(x)/x, x) := Ei(x)]' to `IntegRules'. Now entering `exp(2x)/x' on the stack and typing `a i x' yields `Ei(2 x)'. This new rule will work with Calc's various built-in integration methods (such as integration by substitution) to solve a variety of other problems involving `Ei': For example, now Calc will also be able to integrate `exp(exp(x))' and `ln(ln(x))' (to get `Ei(exp(x))' and `x ln(ln(x)) - Ei(ln(x))', respectively). Your rule may do further integration by calling `integ'. For example, `integtry(twice(u),x) := twice(integ(u))' allows Calc to integrate `twice(sin(x))' to get `twice(-cos(x))'. Note that `integ' was called with only one argument. This notation is allowed only within `IntegRules'; it means "integrate this with respect to the same integration variable." If Calc is unable to integrate `u', the integration that invoked `IntegRules' also fails. Thus integrating `twice(f(x))' fails, returning the unevaluated integral `integ(twice(f(x)), x)'. It is still legal to call `integ' with two or more arguments, however; in this case, if `u' is not integrable, `twice' itself will still be integrated: If the above rule is changed to `... := twice(integ(u,x))', then integrating `twice(f(x))' will yield `twice(integ(f(x),x))'. If a rule instead produces the formula `integsubst(SEXPR, SVAR)', either replacing the top-level `integtry' call or nested anywhere inside the expression, then Calc will apply the substitution `U = SEXPR(SVAR)' to try to integrate the original EXPR. For example, the rule `sqrt(a) := integsubst(sqrt(x),x)' says that if Calc ever finds a square root in the integrand, it should attempt the substitution `u = sqrt(x)'. (This particular rule is unnecessary because Calc always tries "obvious" substitutions where SEXPR actually appears in the integrand.) The variable SVAR may be the same as the VAR that appeared in the call to `integtry', but it need not be. When integrating according to an `integsubst', Calc uses the equation solver to find the inverse of SEXPR (if the integrand refers to VAR anywhere except in subexpressions that exactly match SEXPR). It uses the differentiator to find the derivative of SEXPR and/or its inverse (it has two methods that use one derivative or the other). You can also specify these items by adding extra arguments to the `integsubst' your rules construct; the general form is `integsubst(SEXPR, SVAR, SINV, SPRIME)', where SINV is the inverse of SEXPR (still written as a function of SVAR), and SPRIME is the derivative of SEXPR with respect to SVAR. If you don't specify these things, and Calc is not able to work them out on its own with the information it knows, then your substitution rule will work only in very specific, simple cases. Calc applies `IntegRules' as if by `C-u 1 a r IntegRules'; in other words, Calc stops rewriting as soon as any rule in your rule set succeeds. (If it weren't for this, the `integsubst(sqrt(x),x)' example above would keep on adding layers of `integsubst' calls forever!) Another set of rules, stored in `IntegSimpRules', are applied every time the integrator uses `a s' to simplify an intermediate result. For example, putting the rule `twice(x) := 2 x' into `IntegSimpRules' would tell Calc to convert the `twice' function into a form it knows whenever integration is attempted. One more way to influence the integrator is to define a function with the `Z F' command (*note Algebraic Definitions::.). Calc's integrator automatically expands such functions according to their defining formulas, even if you originally asked for the function to be left unevaluated for symbolic arguments. (Certain other Calc systems, such as the differentiator and the equation solver, also do this.) Sometimes Calc is able to find a solution to your integral, but it expresses the result in a way that is unnecessarily complicated. If this happens, you can either use `integsubst' as described above to try to hint at a more direct path to the desired result, or you can use `IntegAfterRules'. This is an extra rule set that runs after the main integrator returns its result; basically, Calc does an `a r IntegAfterRules' on the result before showing it to you. (It also does an `a s', without `IntegSimpRules', after that to further simplify the result.) For example, Calc's integrator sometimes produces expressions of the form `ln(1+x) - ln(1-x)'; the default `IntegAfterRules' rewrite this into the more readable form `2 arctanh(x)'. Note that, unlike `IntegRules', `IntegSimpRules' and `IntegAfterRules' are applied any number of times until no further changes are possible. Rewriting by `IntegAfterRules' occurs only after the main integrator has finished, not at every step as for `IntegRules' and `IntegSimpRules'. File: calc.info, Node: Numerical Integration, Next: Taylor Series, Prev: Customizing the Integrator, Up: Calculus Numerical Integration --------------------- If you want a purely numerical answer to an integration problem, you can use the `a I' (`calc-num-integral') [`ninteg'] command. This command prompts for an integration variable, a lower limit, and an upper limit. Except for the integration variable, all other variables that appear in the integrand formula must have stored values. (A stored value, if any, for the integration variable itself is ignored.) Numerical integration works by evaluating your formula at many points in the specified interval. Calc uses an "open Romberg" method; this means that it does not evaluate the formula actually at the endpoints (so that it is safe to integrate `sin(x)/x' from zero, for example). Also, the Romberg method works especially well when the function being integrated is fairly smooth. If the function is not smooth, Calc will have to evaluate it at quite a few points before it can accurately determine the value of the integral. Integration is much faster when the current precision is small. It is best to set the precision to the smallest acceptable number of digits before you use `a I'. If Calc appears to be taking too long, press `C-g' to halt it and try a lower precision. If Calc still appears to need hundreds of evaluations, check to make sure your function is well-behaved in the specified interval. It is possible for the lower integration limit to be `-inf' (minus infinity). Likewise, the upper limit may be plus infinity. Calc internally transforms the integral into an equivalent one with finite limits. However, integration to or across singularities is not supported: The integral of `1/sqrt(x)' from 0 to 1 exists (it can be found by Calc's symbolic integrator, for example), but `a I' will fail because the integrand goes to infinity at one of the endpoints. File: calc.info, Node: Taylor Series, Prev: Numerical Integration, Up: Calculus Taylor Series ------------- The `a t' (`calc-taylor') [`taylor'] command computes a power series expansion or Taylor series of a function. You specify the variable and the desired number of terms. You may give an expression of the form `VAR = A' or `VAR - A' instead of just a variable to produce a Taylor expansion about the point A. You may specify the number of terms with a numeric prefix argument; otherwise the command will prompt you for the number of terms. Note that many series expansions have coefficients of zero for some terms, so you may appear to get fewer terms than you asked for. If the `a i' command is unable to find a symbolic integral for a function, you can get an approximation by integrating the function's Taylor series. File: calc.info, Node: Solving Equations, Next: Numerical Solutions, Prev: Calculus, Up: Algebra Solving Equations ================= The `a S' (`calc-solve-for') [`solve'] command rearranges an equation to solve for a specific variable. An equation is an expression of the form `L = R'. For example, the command `a S x' will rearrange `y = 3x + 6' to the form, `x = y/3 - 2'. If the input is not an equation, it is treated like an equation of the form `X = 0'. This command also works for inequalities, as in `y < 3x + 6'. Some inequalities cannot be solved where the analogous equation could be; for example, solving `a < b c' for `b' is impossible without knowing the sign of `c'. In this case, `a S' will produce the result `b != a/c' (using the not-equal-to operator) to signify that the direction of the inequality is now unknown. The inequality `a <= b c' is not even partially solved. *Note Declarations::, for a way to tell Calc that the signs of the variables in a formula are in fact known. Two useful commands for working with the result of `a S' are `a .' (*note Logical Operations::.), which converts `x = y/3 - 2' to `y/3 - 2', and `s l' (*note Let Command::.) which evaluates another formula with `x' set equal to `y/3 - 2'. * Menu: * Multiple Solutions:: * Solving Systems of Equations:: * Decomposing Polynomials:: File: calc.info, Node: Multiple Solutions, Next: Solving Systems of Equations, Prev: Solving Equations, Up: Solving Equations Multiple Solutions ------------------ Some equations have more than one solution. The Hyperbolic flag (`H a S') [`fsolve'] tells the solver to report the fully general family of solutions. It will invent variables `n1', `n2', ..., which represent independent arbitrary integers, and `s1', `s2', ..., which represent independent arbitrary signs (either +1 or -1). If you don't use the Hyperbolic flag, Calc will use zero in place of all arbitrary integers, and plus one in place of all arbitrary signs. Note that variables like `n1' and `s1' are not given any special interpretation in Calc except by the equation solver itself. As usual, you can use the `s l' (`calc-let') command to obtain solutions for various actual values of these variables. For example, `' x^2 = y RET H a S x RET' solves to get `x = s1 sqrt(y)', indicating that the two solutions to the equation are `sqrt(y)' and `-sqrt(y)'. Another way to think about it is that the square-root operation is really a two-valued function; since every Calc function must return a single result, `sqrt' chooses to return the positive result. Then `H a S' doctors this result using `s1' to indicate the full set of possible values of the mathematical square-root. There is a similar phenomenon going the other direction: Suppose we solve `sqrt(y) = x' for `y'. Calc squares both sides to get `y = x^2'. This is correct, except that it introduces some dubious solutions. Consider solving `sqrt(y) = -3': Calc will report `y = 9' as a valid solution, which is true in the mathematical sense of square-root, but false (there is no solution) for the actual Calc positive-valued `sqrt'. This happens for both `a S' and `H a S'. If you store a positive integer in the Calc variable `GenCount', then Calc will generate formulas of the form `as(N)' for arbitrary signs, and `an(N)' for arbitrary integers, where N represents successive values taken by incrementing `GenCount' by one. While the normal arbitrary sign and integer symbols start over at `s1' and `n1' with each new Calc command, the `GenCount' approach will give each arbitrary value a name that is unique throughout the entire Calc session. Also, the arbitrary values are function calls instead of variables, which is advantageous in some cases. For example, you can make a rewrite rule that recognizes all arbitrary signs using a pattern like `as(n)'. The `s l' command only works on variables, but you can use the `a b' (`calc-substitute') command to substitute actual values for function calls like `as(3)'. The `s G' (`calc-edit-GenCount') command is a convenient way to create or edit this variable. Press `M-# M-#' to finish. If you have not stored a value in `GenCount', or if the value in that variable is not a positive integer, the regular `s1'/`n1' notation is used. With the Inverse flag, `I a S' [`finv'] treats the expression on top of the stack as a function of the specified variable and solves to find the inverse function, written in terms of the same variable. For example, `I a S x' inverts `2x + 6' to `x/2 - 3'. You can use both Inverse and Hyperbolic [`ffinv'] to obtain a fully general inverse, as described above. Some equations, specifically polynomials, have a known, finite number of solutions. The `a P' (`calc-poly-roots') [`roots'] command uses `H a S' to solve an equation in general form, then, for all arbitrary-sign variables like `s1', and all arbitrary-integer variables like `n1' for which `n1' only usefully varies over a finite range, it expands these variables out to all their possible values. The results are collected into a vector, which is returned. For example, `roots(x^4 = 1, x)' returns the four solutions `[1, -1, (0, 1), (0, -1)]'. Generally an Nth degree polynomial will always have N roots on the complex plane. (If you have given a `real' declaration for the solution variable, then only the real-valued solutions, if any, will be reported; *note Declarations::..) Note that because `a P' uses `H a S', it is able to deliver symbolic solutions if the polynomial has symbolic coefficients. Also note that Calc's solver is not able to get exact symbolic solutions to all polynomials. Polynomials containing powers up to `x^4' can always be solved exactly; polynomials of higher degree sometimes can be: `x^6 + x^3 + 1' is converted to `(x^3)^2 + (x^3) + 1', which can be solved for `x^3' using the quadratic equation, and then for `x' by taking cube roots. But in many cases, like `x^6 + x + 1', Calc does not know how to rewrite the polynomial into a form it can solve. The `a P' command can still deliver a list of numerical roots, however, provided that symbolic mode (`m s') is not turned on. (If you work with symbolic mode on, recall that the `N' (`calc-eval-num') key is a handy way to reevaluate the formula on the stack with symbolic mode temporarily off.) Naturally, `a P' can only provide numerical roots if the polynomial coefficents are all numbers (real or complex). File: calc.info, Node: Solving Systems of Equations, Next: Decomposing Polynomials, Prev: Multiple Solutions, Up: Solving Equations Solving Systems of Equations ---------------------------- You can also use the commands described above to solve systems of simultaneous equations. Just create a vector of equations, then specify a vector of variables for which to solve. (You can omit the surrounding brackets when entering the vector of variables at the prompt.) For example, putting `[x + y = a, x - y = b]' on the stack and typing `a S x,y RET' produces the vector of solutions `[x = a - (a-b)/2, y = (a-b)/2]'. The result vector will have the same length as the variables vector, and the variables will be listed in the same order there. Note that the solutions are not always simplified as far as possible; the solution for `x' here could be improved by an application of the `a n' command. Calc's algorithm works by trying to eliminate one variable at a time by solving one of the equations for that variable and then substituting into the other equations. Calc will try all the possibilities, but you can speed things up by noting that Calc first tries to eliminate the first variable with the first equation, then the second variable with the second equation, and so on. It also helps to put the simpler (e.g., more linear) equations toward the front of the list. Calc's algorithm will solve any system of linear equations, and also many kinds of nonlinear systems. Normally there will be as many variables as equations. If you give fewer variables than equations (an "over-determined" system of equations), Calc will find a partial solution. For example, typing `a S y RET' with the above system of equations would produce `[y = a - x]'. There are now several ways to express this solution in terms of the original variables; Calc uses the first one that it finds. You can control the choice by adding variable specifiers of the form `elim(V)' to the variables list. This says that V should be eliminated from the equations; the variable will not appear at all in the solution. For example, typing `a S y,elim(x)' would yield `[y = a - (b+a)/2]'. If the variables list contains only `elim' specifiers, Calc simply eliminates those variables from the equations and then returns the resulting set of equations. For example, `a S elim(x)' produces `[a - 2 y = b]'. Every variable eliminated will reduce the number of equations in the system by one. Again, `a S' gives you one solution to the system of equations. If there are several solutions, you can use `H a S' to get a general family of solutions, or, if there is a finite number of solutions, you can use `a P' to get a list. (In the latter case, the result will take the form of a matrix where the rows are different solutions and the columns correspond to the variables you requested.) Another way to deal with certain kinds of overdetermined systems of equations is the `a F' command, which does least-squares fitting to satisfy the equations. *Note Curve Fitting::. File: calc.info, Node: Decomposing Polynomials, Prev: Solving Systems of Equations, Up: Solving Equations Decomposing Polynomials ----------------------- The `poly' function takes a polynomial and a variable as arguments, and returns a vector of polynomial coefficients (constant coefficient first). For example, `poly(x^3 + 2 x, x)' returns `[0, 2, 0, 1]'. If the input is not a polynomial in `x', the call to `poly' is left in symbolic form. If the input does not involve the variable `x', the input is returned in a list of length one, representing a polynomial with only a constant coefficient. The call `poly(x, x)' returns the vector `[0, 1]'. The last element of the returned vector is guaranteed to be nonzero; note that `poly(0, x)' returns the empty vector `[]'. Note also that `x' may actually be any formula; for example, `poly(sin(x)^2 - sin(x) + 3, sin(x))' returns `[3, -1, 1]'. To get the `x^k' coefficient of polynomial `p', use `poly(p, x)_(k+1)'. To get the degree of polynomial `p', use `vlen(poly(p, x)) - 1'. For example, `poly((x+1)^4, x)' returns `[1, 4, 6, 4, 1]', so `poly((x+1)^4, x)_(2+1)' gives the `x^2' coefficient of this polynomial, 6. One important feature of the solver is its ability to recognize formulas which are "essentially" polynomials. This ability is made available to the user through the `gpoly' function, which is used just like `poly': `gpoly(EXPR, VAR)'. If EXPR is a polynomial in some term which includes VAR, then this function will return a vector `[X, C, A]' where X is the term that depends on VAR, C is a vector of polynomial coefficients (like the one returned by `poly'), and A is a multiplier which is usually 1. Basically, `EXPR = A*(C_1 + C_2 X + C_3 X^2 + ...)'. The last element of C is guaranteed to be non-zero, and C will not equal `[1]' (i.e., the trivial decomposition EXPR = X is not considered a polynomial). One side effect is that `gpoly(x, x)' and `gpoly(6, x)', both of which might be expected to recognize their arguments as polynomials, will not because the decomposition is considered trivial. For example, `gpoly((x-2)^2, x)' returns `[x, [4, -4, 1], 1]', since the expanded form of this polynomial is `4 - 4 x + x^2'. The term X may itself be a polynomial in VAR. This is done to reduce the size of the C vector. For example, `gpoly(x^4 + x^2 - 1, x)' returns `[x^2, [-1, 1, 1], 1]', since a quadratic polynomial in `x^2' is easier to solve than a quartic polynomial in `x'. A few more examples of the kinds of polynomials `gpoly' can discover: sin(x) - 1 [sin(x), [-1, 1], 1] x + 1/x - 1 [x, [1, -1, 1], 1/x] x + 1/x [x^2, [1, 1], 1/x] x^3 + 2 x [x^2, [2, 1], x] x + x^2:3 + sqrt(x) [x^1:6, [1, 1, 0, 1], x^1:2] x^(2a) + 2 x^a + 5 [x^a, [5, 2, 1], 1] (exp(-x) + exp(x)) / 2 [e^(2 x), [0.5, 0.5], e^-x] The `poly' and `gpoly' functions accept a third integer argument which specifies the largest degree of polynomial that is acceptable. If this is `n', then only C vectors of length `n+1' or less will be returned. Otherwise, the `poly' or `gpoly' call will remain in symbolic form. For example, the equation solver can handle quartics and smaller polynomials, so it calls `gpoly(EXPR, VAR, 4)' to discover whether EXPR can be treated by its linear, quadratic, cubic, or quartic formulas. The `pdeg' function computes the degree of a polynomial; `pdeg(p,x)' is the highest power of `x' that appears in `p'. This is the same as `vlen(poly(p,x))-1', but is much more efficient. If `p' is constant with respect to `x', then `pdeg(p,x) = 0'. If `p' is not a polynomial in `x' (e.g., `pdeg(2 cos(x), x)', the function remains unevaluated. It is possible to omit the second argument `x', in which case `pdeg(p)' returns the highest total degree of any term of the polynomial, counting all variables that appear in `p'. Note that `pdeg(c) = pdeg(c,x) = 0' for any nonzero constant `c'; the degree of the constant zero is considered to be `-inf' (minus infinity). The `plead' function finds the leading term of a polynomial. Thus `plead(p,x)' is equivalent to `poly(p,x)_vlen(poly(p,x))', though again more efficient. In particular, `plead((2x+1)^10, x)' returns 1024 without expanding out the list of coefficients. The value of `plead(p,x)' will be zero only if `p = 0'. The `pcont' function finds the "content" of a polynomial. This is the greatest common divisor of all the coefficients of the polynomial. With two arguments, `pcont(p,x)' effectively uses `poly(p,x)' to get a list of coefficients, then uses `pgcd' (the polynomial GCD function) to combine these into an answer. For example, `pcont(4 x y^2 + 6 x^2 y, x)' is `2 y'. The content is basically the "biggest" polynomial that can be divided into `p' exactly. The sign of the content is the same as the sign of the leading coefficient. With only one argument, `pcont(p)' computes the numerical content of the polynomial, i.e., the `gcd' of the numerical coefficients of all the terms in the formula. Note that `gcd' is defined on rational numbers as well as integers; it computes the `gcd' of the numerators and the `lcm' of the denominators. Thus `pcont(4:3 x y^2 + 6 x^2 y)' returns 2:3. Dividing the polynomial by this number will clear all the denominators, as well as dividing by any common content in the numerators. The numerical content of a polynomial is negative only if all the coefficients in the polynomial are negative. The `pprim' function finds the "primitive part" of a polynomial, which is simply the polynomial divided (using `pdiv' if necessary) by its content. If the input polynomial has rational coefficients, the result will have integer coefficients in simplest terms. File: calc.info, Node: Numerical Solutions, Next: Curve Fitting, Prev: Solving Equations, Up: Algebra Numerical Solutions =================== Not all equations can be solved symbolically. The commands in this section use numerical algorithms that can find a solution to a specific instance of an equation to any desired accuracy. Note that the numerical commands are slower than their algebraic cousins; it is a good idea to try `a S' before resorting to these commands. (*Note Curve Fitting::, for some other, more specialized, operations on numerical data.) * Menu: * Root Finding:: * Minimization:: * Numerical Systems of Equations:: File: calc.info, Node: Root Finding, Next: Minimization, Prev: Numerical Solutions, Up: Numerical Solutions Root Finding ------------ The `a R' (`calc-find-root') [`root'] command finds a numerical solution (or "root") of an equation. (This command treats inequalities the same as equations. If the input is any other kind of formula, it is interpreted as an equation of the form `X = 0'.) The `a R' command requires an initial guess on the top of the stack, and a formula in the second-to-top position. It prompts for a solution variable, which must appear in the formula. All other variables that appear in the formula must have assigned values, i.e., when a value is assigned to the solution variable and the formula is evaluated with `=', it should evaluate to a number. Any assigned value for the solution variable itself is ignored and unaffected by this command. When the command completes, the initial guess is replaced on the stack by a vector of two numbers: The value of the solution variable that solves the equation, and the difference between the lefthand and righthand sides of the equation at that value. Ordinarily, the second number will be zero or very nearly zero. (Note that Calc uses a slightly higher precision while finding the root, and thus the second number may be slightly different from the value you would compute from the equation yourself.) The `v h' (`calc-head') command is a handy way to extract the first element of the result vector, discarding the error term. The initial guess can be a real number, in which case Calc searches for a real solution near that number, or a complex number, in which case Calc searches the whole complex plane near that number for a solution, or it can be an interval form which restricts the search to real numbers inside that interval. Calc tries to use `a d' to take the derivative of the equation. If this succeeds, it uses Newton's method. If the equation is not differentiable Calc uses a bisection method. (If Newton's method appears to be going astray, Calc switches over to bisection if it can, or otherwise gives up. In this case it may help to try again with a slightly different initial guess.) If the initial guess is a complex number, the function must be differentiable. If the formula (or the difference between the sides of an equation) is negative at one end of the interval you specify and positive at the other end, the root finder is guaranteed to find a root. Otherwise, Calc subdivides the interval into small parts looking for positive and negative values to bracket the root. When your guess is an interval, Calc will not look outside that interval for a root. The `H a R' [`wroot'] command is similar to `a R', except that if the initial guess is an interval for which the function has the same sign at both ends, then rather than subdividing the interval Calc attempts to widen it to enclose a root. Use this mode if you are not sure if the function has a root in your interval. If the function is not differentiable, and you give a simple number instead of an interval as your initial guess, Calc uses this widening process even if you did not type the Hyperbolic flag. (If the function *is* differentiable, Calc uses Newton's method which does not require a bounding interval in order to work.) If Calc leaves the `root' or `wroot' function in symbolic form on the stack, it will normally display an explanation for why no root was found. If you miss this explanation, press `w' (`calc-why') to get it back. File: calc.info, Node: Minimization, Next: Numerical Systems of Equations, Prev: Root Finding, Up: Numerical Solutions Minimization ------------ The `a N' (`calc-find-minimum') [`minimize'] command finds a minimum value for a formula. It is very similar in operation to `a R' (`calc-find-root'): You give the formula and an initial guess on the stack, and are prompted for the name of a variable. The guess may be either a number near the desired minimum, or an interval enclosing the desired minimum. The function returns a vector containing the value of the the variable which minimizes the formula's value, along with the minimum value itself. Note that this command looks for a *local* minimum. Many functions have more than one minimum; some, like `x sin(x)', have infinitely many. In fact, there is no easy way to define the "global" minimum of `x sin(x)' but Calc can still locate any particular local minimum for you. Calc basically goes downhill from the initial guess until it finds a point at which the function's value is greater both to the left and to the right. Calc does not use derivatives when minimizing a function. If your initial guess is an interval and it looks like the minimum occurs at one or the other endpoint of the interval, Calc will return that endpoint only if that endpoint is closed; thus, minimizing `17 x' over `[2..3]' will return `[2, 38]', but minimizing over `(2..3]' would report no minimum found. In general, you should use closed intervals to find literally the minimum value in that range of `x', or open intervals to find the local minimum, if any, that happens to lie in that range. Most functions are smooth and flat near their minimum values. Because of this flatness, if the current precision is, say, 12 digits, the variable can only be determined meaningfully to about six digits. Thus you should set the precision to twice as many digits as you need in your answer. The `H a N' [`wminimize'] command, analogously to `H a R', expands the guess interval to enclose a minimum rather than requiring that the minimum lie inside the interval you supply. The `a X' (`calc-find-maximum') [`maximize'] and `H a X' [`wmaximize'] commands effectively minimize the negative of the formula you supply. The formula must evaluate to a real number at all points inside the interval (or near the initial guess if the guess is a number). If the initial guess is a complex number the variable will be minimized over the complex numbers; if it is real or an interval it will be minimized over the reals. File: calc.info, Node: Numerical Systems of Equations, Prev: Minimization, Up: Numerical Solutions Systems of Equations -------------------- The `a R' command can also solve systems of equations. In this case, the equation should instead be a vector of equations, the guess should instead be a vector of numbers (intervals are not supported), and the variable should be a vector of variables. You can omit the brackets while entering the list of variables. Each equation must be differentiable by each variable for this mode to work. The result will be a vector of two vectors: The variable values that solved the system of equations, and the differences between the sides of the equations with those variable values. There must be the same number of equations as variables. Since only plain numbers are allowed as guesses, the Hyperbolic flag has no effect when solving a system of equations. It is also possible to minimize over many variables with `a N' (or maximize with `a X'). Once again the variable name should be replaced by a vector of variables, and the initial guess should be an equal-sized vector of initial guesses. But, unlike the case of multidimensional `a R', the formula being minimized should still be a single formula, *not* a vector. Beware that multidimensional minimization is currently *very* slow. File: calc.info, Node: Curve Fitting, Next: Summations, Prev: Numerical Solutions, Up: Algebra Curve Fitting ============= The `a F' command fits a set of data to a "model formula", such as `y = m x + b' where `m' and `b' are parameters to be determined. For a typical set of measured data there will be no single `m' and `b' that exactly fit the data; in this case, Calc chooses values of the parameters that provide the closest possible fit. * Menu: * Linear Fits:: * Polynomial and Multilinear Fits:: * Error Estimates for Fits:: * Standard Nonlinear Models:: * Curve Fitting Details:: * Interpolation:: File: calc.info, Node: Linear Fits, Next: Polynomial and Multilinear Fits, Prev: Curve Fitting, Up: Curve Fitting Linear Fits ----------- The `a F' (`calc-curve-fit') [`fit'] command attempts to fit a set of data (`x' and `y' vectors of numbers) to a straight line, polynomial, or other function of `x'. For the moment we will consider only the case of fitting to a line, and we will ignore the issue of whether or not the model was in fact a good fit for the data. In a standard linear least-squares fit, we have a set of `(x,y)' data points that we wish to fit to the model `y = m x + b' by adjusting the parameters `m' and `b' to make the `y' values calculated from the formula be as close as possible to the actual `y' values in the data set. (In a polynomial fit, the model is instead, say, `y = a x^3 + b x^2 + c x + d'. In a multilinear fit, we have data points of the form `(x_1,x_2,x_3,y)' and our model is `y = a x_1 + b x_2 + c x_3 + d'. These will be discussed later.) In the model formula, variables like `x' and `x_2' are called the "independent variables", and `y' is the "dependent variable". Variables like `m', `a', and `b' are called the "parameters" of the model. The `a F' command takes the data set to be fitted from the stack. By default, it expects the data in the form of a matrix. For example, for a linear or polynomial fit, this would be a 2xN matrix where the first row is a list of `x' values and the second row has the corresponding `y' values. For the multilinear fit shown above, the matrix would have four rows (`x_1', `x_2', `x_3', and `y', respectively). If you happen to have an Nx2 matrix instead of a 2xN matrix, just press `v t' first to transpose the matrix. After you type `a F', Calc prompts you to select a model. For a linear fit, press the digit `1'. Calc then prompts for you to name the variables. By default it chooses high letters like `x' and `y' for independent variables and low letters like `a' and `b' for parameters. (The dependent variable doesn't need a name.) The two kinds of variables are separated by a semicolon. Since you generally care more about the names of the independent variables than of the parameters, Calc also allows you to name only those and let the parameters use default names. For example, suppose the data matrix [ [ 1, 2, 3, 4, 5 ] [ 5, 7, 9, 11, 13 ] ] is on the stack and we wish to do a simple linear fit. Type `a F', then `1' for the model, then `RET' to use the default names. The result will be the formula `3 + 2 x' on the stack. Calc has created the model expression `a + b x', then found the optimal values of `a' and `b' to fit the data. (In this case, it was able to find an exact fit.) Calc then substituted those values for `a' and `b' in the model formula. The `a F' command puts two entries in the trail. One is, as always, a copy of the result that went to the stack; the other is a vector of the actual parameter values, written as equations: `[a = 3, b = 2]', in case you'd rather read them in a list than pick them out of the formula. (You can type `t y' to move this vector to the stack; *note Trail Commands::..) Specifying a different independent variable name will affect the resulting formula: `a F 1 k RET' produces `3 + 2 k'. Changing the parameter names (say, `a F 1 k;b,m RET') will affect the equations that go into the trail. To see what happens when the fit is not exact, we could change the number 13 in the data matrix to 14 and try the fit again. The result 2.6 + 2.2 x Evaluating this formula, say with `v x 5 RET TAB V M $ RET', shows a reasonably close match to the y-values in the data. [4.8, 7., 9.2, 11.4, 13.6] Since there is no line which passes through all the N data points, Calc has chosen a line that best approximates the data points using the method of least squares. The idea is to define the "chi-square" error measure chi^2 = sum((y_i - (a + b x_i))^2, i, 1, N) which is clearly zero if `a + b x' exactly fits all data points, and increases as various `a + b x_i' values fail to match the corresponding `y_i' values. There are several reasons why the summand is squared, one of them being to ensure that `chi^2 >= 0'. Least-squares fitting simply chooses the values of `a' and `b' for which the error `chi^2' is as small as possible. Other kinds of models do the same thing but with a different model formula in place of `a + b x_i'. A numeric prefix argument causes the `a F' command to take the data in some other form than one big matrix. A positive argument N will take N items from the stack, corresponding to the N rows of a data matrix. In the linear case, N must be 2 since there is always one independent variable and one dependent variable. A prefix of zero or plain `C-u' is a compromise; Calc takes two items from the stack, an N-row matrix of `x' values, and a vector of `y' values. If there is only one independent variable, the `x' values can be either a one-row matrix or a plain vector, in which case the `C-u' prefix is the same as a `C-u 2' prefix. File: calc.info, Node: Polynomial and Multilinear Fits, Next: Error Estimates for Fits, Prev: Linear Fits, Up: Curve Fitting Polynomial and Multilinear Fits ------------------------------- To fit the data to higher-order polynomials, just type one of the digits `2' through `9' when prompted for a model. For example, we could fit the original data matrix from the previous section (with 13, not 14) to a parabola instead of a line by typing `a F 2 RET'. 2.00000000001 x - 1.5e-12 x^2 + 2.99999999999 Note that since the constant and linear terms are enough to fit the data exactly, it's no surprise that Calc chose a tiny contribution for `x^2'. (The fact that it's not exactly zero is due only to roundoff error. Since our data are exact integers, we could get an exact answer by typing `m f' first to get fraction mode. Then the `x^2' term would vanish altogether. Usually, though, the data being fitted will be approximate floats so fraction mode won't help.) Doing the `a F 2' fit on the data set with 14 instead of 13 gives a much larger `x^2' contribution, as Calc bends the line slightly to improve the fit. 0.142857142855 x^2 + 1.34285714287 x + 3.59999999998 An important result from the theory of polynomial fitting is that it is always possible to fit N data points exactly using a polynomial of degree N-1, sometimes called an "interpolating polynomial". Using the modified (14) data matrix, a model number of 4 gives a polynomial that exactly matches all five data points: 0.04167 x^4 - 0.4167 x^3 + 1.458 x^2 - 0.08333 x + 4. The actual coefficients we get with a precision of 12, like `0.0416666663588', clearly suffer from loss of precision. It is a good idea to increase the working precision to several digits beyond what you need when you do a fitting operation. Or, if your data are exact, use fraction mode to get exact results. You can type `i' instead of a digit at the model prompt to fit the data exactly to a polynomial. This just counts the number of columns of the data matrix to choose the degree of the polynomial automatically. Fitting data "exactly" to high-degree polynomials is not always a good idea, though. High-degree polynomials have a tendency to wiggle uncontrollably in between the fitting data points. Also, if the exact-fit polynomial is going to be used to interpolate or extrapolate the data, it is numerically better to use the `a p' command described below. *Note Interpolation::. Another generalization of the linear model is to assume the `y' values are a sum of linear contributions from several `x' values. This is a "multilinear" fit, and it is also selected by the `1' digit key. (Calc decides whether the fit is linear or multilinear by counting the rows in the data matrix.) Given the data matrix, [ [ 1, 2, 3, 4, 5 ] [ 7, 2, 3, 5, 2 ] [ 14.5, 15, 18.5, 22.5, 24 ] ] the command `a F 1 RET' will call the first row `x' and the second row `y', and will fit the values in the third row to the model `a + b x + c 8. + 3. x + 0.5 y Calc can do multilinear fits with any number of independent variables (i.e., with any number of data rows). Yet another variation is "homogeneous" linear models, in which the constant term is known to be zero. In the linear case, this means the model formula is simply `a x'; in the multilinear case, the model might be `a x + b y + c z'; and in the polynomial case, the model could be `a x + b x^2 + c x^3'. You can get a homogeneous linear or multilinear model by pressing the letter `h' followed by a regular model key, like `1' or `2'. It is certainly possible to have other constrained linear models, like `2.3 + a x' or `a - 4 x'. While there is no single key to select models like these, a later section shows how to enter any desired model by hand. In the first case, for example, you would enter `a F ' 2.3 + a x'. Another class of models that will work but must be entered by hand are multinomial fits, e.g., `a + b x + c y + d x^2 + e y^2 + f x y'. File: calc.info, Node: Error Estimates for Fits, Next: Standard Nonlinear Models, Prev: Polynomial and Multilinear Fits, Up: Curve Fitting Error Estimates for Fits ------------------------ With the Hyperbolic flag, `H a F' [`efit'] performs the same fitting operation as `a F', but reports the coefficients as error forms instead of plain numbers. Fitting our two data matrices (first with 13, then with 14) to a line with `H a F' gives the results, 3. + 2. x 2.6 +/- 0.382970843103 + 2.2 +/- 0.115470053838 x In the first case the estimated errors are zero because the linear fit is perfect. In the second case, the errors are nonzero but moderately small, because the data are still very close to linear. It is also possible for the *input* to a fitting operation to contain error forms. The data values must either all include errors or all be plain numbers. Error forms can go anywhere but generally go on the numbers in the last row of the data matrix. If the last row contains error forms `y_i +/- sigma_i', then the `chi^2' statistic is chi^2 = sum(((y_i - (a + b x_i)) / sigma_i)^2, i, 1, N) so that data points with larger error estimates contribute less to the fitting operation. If there are error forms on other rows of the data matrix, all the errors for a given data point are combined; the square root of the sum of the squares of the errors forms the `sigma_i' used for the data point. Both `a F' and `H a F' can accept error forms in the input matrix, although if you are concerned about error analysis you will probably use `H a F' so that the output also contains error estimates. If the input contains error forms but all the `sigma_i' values are the same, it is easy to see that the resulting fitted model will be the same as if the input did not have error forms at all (`chi^2' is simply scaled uniformly by `1 / sigma^2', which doesn't affect where it has a minimum). But there *will* be a difference in the estimated errors of the coefficients reported by `H a F'. Consult any text on statistical modelling of data for a discussion of where these error estimates come from and how they should be interpreted. With the Inverse flag, `I a F' [`xfit'] produces even more information. The result is a vector of six items: 1. The model formula with error forms for its coefficients or parameters. This is the result that `H a F' would have produced. 2. A vector of "raw" parameter values for the model. These are the polynomial coefficients or other parameters as plain numbers, in the same order as the parameters appeared in the final prompt of the `I a F' command. For polynomials of degree `d', this vector will have length `M = d+1' with the constant term first. 3. The covariance matrix `C' computed from the fit. This is an MxM symmetric matrix; the diagonal elements `C_j_j' are the variances `sigma_j^2' of the parameters. The other elements are covariances `sigma_i_j^2' that describe the correlation between pairs of parameters. (A related set of numbers, the "linear correlation coefficients" `r_i_j', are defined as `sigma_i_j^2 / sigma_i sigma_j'.) 4. A vector of `M' "parameter filter" functions whose meanings are described below. If no filters are necessary this will instead be an empty vector; this is always the case for the polynomial and multilinear fits described so far. 5. The value of `chi^2' for the fit, calculated by the formulas shown above. This gives a measure of the quality of the fit; statisticians consider `chi^2 = N - M' to indicate a moderately good fit (where again `N' is the number of data points and `M' is the number of parameters). 6. A measure of goodness of fit expressed as a probability `Q'. This is computed from the `utpc' probability distribution function using `chi^2' with `N - M' degrees of freedom. A value of 0.5 implies a good fit; some texts recommend that often `Q = 0.1' or even 0.001 can signify an acceptable fit. In particular, `chi^2' statistics assume the errors in your inputs follow a normal (Gaussian) distribution; if they don't, you may have to accept smaller values of `Q'. The `Q' value is computed only if the input included error estimates. Otherwise, Calc will report the symbol `nan' for `Q'. The reason is that in this case the `chi^2' value has effectively been used to estimate the original errors in the input, and thus there is no redundant information left over to use for a confidence test.