Artifact c593dd14cb3092bd50ce2737020e2226d6deff5528064118b2384d48de5589d6:
- Executable file
r37/packages/numeric/numeric.txt
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 15636) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/numeric/numeric.txt
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 15636) [annotate] [blame] [check-ins using]
NUMERIC Herbert Melenk Konrad--Zuse--Zentrum fuer Informationstechnik Berlin E--mail: Melenk@sc.zib--berlin.de The NUMERIC package implements some numerical (approximative) algorithms for REDUCE, based on the REDUCE rounded mode arithmetic. These algorithms are implemented for standard cases. They should not be called for ill-conditioned problems; please use standard mathematical libraries for these. 1 Syntax 1.1 Intervals, Starting Points Intervals are generally coded as lower bound and upper bound connected by the operator `..', usually associated to a variable in an equation. E.g. x= (2.5 .. 3.5) means that the variable x is taken in the range from 2.5 up to 3.5. Note, that the bounds can be algebraic expressions, which, however, must evaluate to numeric results. In cases where an interval is returned as the result, the lower and upper bounds can be extracted by the PART operator as the first and second part respectively. A starting point is specified by an equation with a numeric righthand side, e.g. x=3.0 If for multivariate applications several coordinates must be specified by intervals or as a starting point, these specifications 1 2 MINIMA 2 can be collected in one parameter (which is then a list) or they can be given as separate parameters alternatively. The list form is more appropriate when the parameters are built from other REDUCE calculations in an automatical style, while the flat form is more convenient for direct interactive input. 1.2 Accuracy Control The keyword parameters accuracy=a and iterations=i, where a and i must be positive integer numbers, control the iterative algorithms: -a the iteration is continued until the local error is below 10 ; if that is impossible within i steps, the iteration is terminated with an error message. The values reached so far are then returned as the result. 1.3 tracing Normally the algorithms produce only a minimum of printed output during their operation. In cases of an unsuccessful or unexpected long operation a trace of the iteration can be printed by setting on trnumeric; 2 Minima The Fletcher Reeves version of the steepest descent algorithms is used to find the minimum of a function of one or more variables. The function must have continuous partial derivatives with respect to all variables. The starting point of the search can be specified; if not, random values are taken instead. The steepest descent algorithms in general find only local minima. Syntax: NUMMIN (exp,var [=val ][,var [=val ]...] 1 1 2 2 [,accuracy=a][,iterations=i]) or 3 ROOTS OF FUNCTIONS/ SOLUTIONS OF EQUATIONS 3 NUMMIN (exp,{var [=val ][,var [=val ]...]} 1 1 2 2 [,accuracy=a][,iterations=i]) where exp is a function expression, var ,var ,... are the variables in exp and val ,val ,... are 1 2 1 2 the (optional) start values. MIN tries to find the next local minimum along the descending path starting at the given point. The result is a list with the minimum function value as first element followed by a list of equations, where the variables are equated to the coordinates of the result point. Examples: num_min(sin(x)+x/5, x); {4.9489585606,{X=29.643767785}} num_min(sin(x)+x/5, x=0); { - 1.3342267466,{X= - 1.7721582671}} % Rosenbrock function (well known as hard to minimize). fktn := 100*(x1**2-x2)**2 + (1-x1)**2; num_min(fktn, x1=-1.2, x2=1, iterations=200); {0.00000021870228295,{X1=0.99953284494,X2=0.99906807238}} 3 Roots of Functions/ Solutions of Equations An adaptively damped Newton iteration is used to find an approximative zero of a function, a function vector or the solution of an equation or an equation system. Equations are internally converted to a difference of lhs and rhs such that the Newton method (=zero detection) can be applied. The expressions must have continuous derivatives for all variables. A starting point for the iteration can be given. If not given, random values are taken 3 ROOTS OF FUNCTIONS/ SOLUTIONS OF EQUATIONS 4 instead. If the number of forms is not equal to the number of variables, the Newton method cannot be applied. Then the minimum of the sum of absolute squares is located instead. With ON COMPLEX solutions with imaginary parts can be found, if either the expression(s) or the starting point contain a nonzero imaginary part. Syntax: NUM_SOLVE (exp ,var [=val ][,accuracy=a][,iterations=i]) or 1 1 1 NUM_SOLVE ({exp ,... ,exp },var [=val ],... ,var [=val ] 1 n 1 1 1 n [,accuracy=a][,iterations=i]) or NUM_ SOLVE ({exp ,... ,exp },{var [=val ],... ,var [=val ]} 1 n 1 1 1 n [,accuracy=a][,iterations=i]) where exp ,... ,exp are function expressions, 1 n var ,... ,var are the variables, 1 n val ,... ,val are optional start values. 1 n SOLVE tries to find a zero/solution of the expression(s). Result is a list of equations, where the variables are equated to the coordinates of the result point. The Jacobian matrix is stored as side effect the shared variable JACOBIAN. Example: num_solve({sin x=cos y, x + y = 1},{x=1,y=2}); {X= - 1.8561957251,Y=2.856195584} jacobian; [COS(X) SIN(Y)] [ ] 4 INTEGRALS 5 [ 1 1 ] 4 Integrals For the numerical evaluation of univariate integrals over a finite interval the following strategy is used: 1. If the function has an antiderivative in close form which is bounded in the integration interval, this is used. 2. Otherwise a Chebyshev approximation is computed, starting with order 20, eventually up to order 80. If that is recognized as sufficiently convergent it is used for computing the integral by directly integrating the coefficient sequence. 3. If none of these methods is successful, an adaptive multilevel quadrature algorithm is used. For multivariate integrals only the adaptive quadrature is used. This algorithm tolerates isolated singularities. The value iterations here limits the number of local interval intersection levels. Accuracy is a measure for the relative total discretization error (comparison of order 1 and order 2 approximations). Syntax: NUMINT (exp,var =(l ..u )[,var =(l ..u )... ] 1 1 1 2 2 2 [,accuracy=a][,iterations=i]) where exp is the function to be integrated, var ,var ,... are the integration variables, 1 2 l ,l ,... are the lower bounds, 1 2 u ,u ,... are the upper bounds. 1 2 Result is the value of the integral. Example: num_int(sin x,x=(0 .. pi)); 2.0000010334 5 ORDINARY DIFFERENTIAL EQUATIONS 6 5 Ordinary Differential Equations A Runge-Kutta method of order 3 finds an approximate graph for the solution of a ordinary differential equation real initial value problem. Syntax: NUMODESOLVE (exp,depvar=dv,indepvar=(from..to) [,accuracy=a][,iterations=i]) where exp is the differential expression/equation, depvar is an identifier representing the dependent variable (function to be found), indepvar is an identifier representing the independent variable, exp is an equation (or an expression implicitly set to zero) which contains the first derivative of depvar wrt indepvar, from is the starting point of integration, to is the endpoint of integration (allowed to be below from), dv is the initial value of depvar in the point indepvar=from. The ODE exp is converted into an explicit form, which then is used for a Runge Kutta iteration over the given range. The number of steps is controlled by the value of i (default: 20). If the steps are too coarse to reach the desired accuracy in the neighborhood of the starting point, the number is increased automatically. Result is a list of pairs, each representing a point of the approximate solution of the ODE problem. Example: num_odesolve(df(y,x)=y,y=1,x=(0 .. 1), iterations=5); {{0.0,1.0},{0.2,1.2214},{0.4,1.49181796},{0.6,1.8221064563}, 6 BOUNDS OF A FUNCTION 7 {0.8,2.2255208258},{1.0,2.7182511366}} Remarks: -- If in exp the differential is not isolated on the lefthand side, please ensure that the dependent variable is explicitly declared using a DEPEND statement, e.g. depend y,x; otherwise the formal derivative will be computed to zero by REDUCE. -- The REDUCE package SOLVE is used to convert the form into an explicit ODE. If that process fails or has no unique result, the evaluation is stopped with an error message. 6 Bounds of a Function Upper and lower bounds of a real valued function over an interval or a rectangular multivariate domain are computed by the operator BOUNDS. The algorithmic basis is the computation with inequalities: starting from the interval(s) of the variables, the bounds are propagated in the expression using the rules for inequality computation. Some knowledge about the behavior of special functions like ABS, SIN, COS, EXP, LOG, fractional exponentials etc. is integrated and can be evaluated if the operator BOUNDS is called with rounded mode on (otherwise only algebraic evaluation rules are available). If BOUNDS finds a singularity within an interval, the evaluation is stopped with an error message indicating the problem part of the expression. Syntax: BOUNDS (exp,var =(l ..u )[,var =(l ..u )...]) 1 1 1 2 2 2 BOUNDS (exp,{var =(l ..u )[,var =(l ..u )...]}) 1 1 1 2 2 2 where exp is the function to be investigated, 7 CHEBYSHEV CURVE FITTING 8 var ,var ,... are the variables of exp, 1 2 l ,l ,... and u ,u ,... specify the area (intervals). 1 2 1 2 BOUNDS computes upper and lower bounds for the expression in the given area. An interval is returned. Example: bounds(sin x,x=(1 .. 2)); {-1,1} on rounded; bounds(sin x,x=(1 .. 2)); 0.84147098481 .. 1 bounds(x**2+x,x=(-0.5 .. 0.5)); - 0.25 .. 0.75 7 Chebyshev Curve Fitting The operator family Chebyshev-... implements approximation and (a,b) evaluation of functions by the Chebyshev method. Let T (x) be the Chebyshev polynomial of order n transformed to the ninterval (a,b). Then a function f(x) can be approximated in (a,b) by a series N (a,b) f(x)?si=0 c T (x) i i The operator Chebyshev-fit computes this approximation and returns a list, which has as first element the sum expressed as a polynomial 7 CHEBYSHEV CURVE FITTING 9 and as second element the sequence of Chebyshev coefficients c . i Chebyshevdf and Chebyshevint transform a Chebyshev coefficient list into the coefficients of the corresponding derivative or integral respectively. For evaluating a Chebyshev approximation at a given point in the basic interval the operator Chebysheveval can be used. Note that Chebyshev-eval is based on a recurrence relation which is in general more stable than a direct evaluation of the complete polynomial. CHEBYSHEVFIT (fcn,var=(lo..hi),n) CHEBYSHEVEVAL (coeffs,var=(lo..hi),var=pt) CHEBYSHEVDF (coeffs,var=(lo..hi)) CHEBYSHEVINT (coeffs,var=(lo..hi)) where fcn is an algebraic expression (the function to be fitted), var is the variable of fcn, lo and hi are numerical real values which describe an interval (lo<hi), n is the approximation order,an integer >0, set to 20 if missing, pt is a numerical value in the interval and coeffs is a series of Chebyshev coefficients, computed by one of CHEBYSHEVCOEFF, -DF or -INT. Example: on rounded; w:=chebyshev_fit(sin x/x,x=(1 .. 3),5); 3 2 w := {0.03824*x - 0.2398*x + 0.06514*x + 0.9778, {0.8991,-0.4066,-0.005198,0.009464,-0.00009511}} chebyshev_eval(second w, x=(1 .. 3), x=2.1); 0.4111 8 GENERAL CURVE FITTING 10 8 General Curve Fitting The operator NUMFIT finds for a set of points the linear combination of a given set of functions (function basis) which approximates the points best under the objective of the least squares criterion (minimum of the sum of the squares of the deviation). The solution is found as zero of the gradient vector of the sum of squared errors. Syntax: NUMFIT (vals,basis,var=pts) where vals is a list of numeric values, var is a variable used for the approximation, pts is a list of coordinate values which correspond to var, basis is a set of functions varying in var which is used for the approximation. The result is a list containing as first element the function which approximates the given values, and as second element a list of coefficients which were used to build this function from the basis. Example: % approximate a set of factorials by a polynomial pts:=for i:=1 step 1 until 5 collect i$ vals:=for i:=1 step 1 until 5 collect for j:=1:i product j$ num_fit(vals,{1,x,x**2},x=pts); 2 {14.571428571*X - 61.428571429*X + 54.6,{54.6, - 61.428571429,14.571428571}} num_fit(vals,{1,x,x**2,x**3,x**4},x=pts); 4 3 8 GENERAL CURVE FITTING 11 {2.2083333234*X - 20.249999879*X 2 + 67.791666154*X - 93.749999133*X + 44.999999525, {44.999999525, - 93.749999133,67.791666154, - 20.249999879,2.2083333234}}