Numeric Library / Reference
Omikron Basic on the Internet: http://www.berkhan.de
General Information -turn page-Table of Contents
|
Numeric Library / Reference |
|
Omikron Basic on the Internet: http://www.berkhan.de |
|
|
2. Numeric Library Reference
| This section will explain the procedures and functions of the Numeric Library. The theoretical background, however, cannot be explained at this point. To employ the Library in a purposeful manner, we thus recommend the use of additional references. We especially recommend the book "Numerical Recipes," which explains many aspects of numerical mathematics in practice-oriented detail. |
| 2.1 Library Logging On and Off | ||||
|
||||
| Num_Exit | |
| Call this procedure once at the end of your program. You cannot use the Numeric Library afterwards. |
| Numeric | |
| Only a copyright message of the Numeric Library is returned. |
| Num_Mode M | |
| M | Operating mode: M=0 : Display error messages. M=1 : Do not display error messages. |
| With this procedure, you can set the operating mode of the Numeric Library. If an error occurs, it will usually be displayed in an alert box. You can suppress the display with M=1, for example, to perform an error handling routine yourself. | |
| FN Num_Error | |
| Returns the number of the last error that occurred. If you have turned off the automatic error display
using Num_Mode 1, you may use this function to determine which error has occurred. The following error numbers are possible: 0 : No error. 1 : FN Li#(N,X#) works only with N=1 or N=2. 2 : FN Bessel#(N,X#) works only with N>=0. 3 : FN Weber#(N,X#) works only with N>=0. 4 : FN Mod_Bessel#(N,X#) works only with N>=0. 5 : FN Macdo#(N,X#) works only with N>=0. 6 : FN Zero# could not find a zero point. 7 : FN Minimum# could not find a minimum. 8 : FN Maximum# could not find a maximum. 9 : PROC Eigen works only with symmetrical matrices. 10 : Eigenvalue calculation in PROC Eigen does not converge. 11 : The equation system in PROC Gauss has no solution, because the matrix is singular. 12 : PROC Gauss_Seidel does not converge. 13 : The number of equations in PROC Band does not serve a purpose. 14 : PROC Newton_Sys does not converge. 15 : Too few supporting points in PROC Spline_Int. 16 : PROC Gauss_Fit does not converge. 17 : Error in PROC Fft. Use only even, positive numbers for N! 18 : PROC Mean_Approx does not converge. 19 : Too few data in PROC Derive (N>=7). 20 : Too few data in FN Newton_Cotes (N>=5). 21 : FN Romberg does not converge. 22 : FN Romberg_2 does not converge. 23 : Invalid parameter in FN Gauss_Leg# (2<=N<=15). |
| 2.2 Special Functions |
| Attention: Not all functions can be exactly calculated.
Many functions have singularities in a few locations and have to be calculated recursively or by using polynomial
interpolations. Since the accuracy of floating point numbers is limited, it can come to a built up of errors, causing
completely false results. For example, the modified Bessel function returns for N=100 in the area of X#=0 completely
inaccurate results. However, the functions are quite acceptable to use for normal applications, where they return
rather accurate function values without problems. Tip: In the test phase, set the control word COMPILER "DEBUG ON," which has the result of notifying you of problematic areas, because all processor exceptions are reported. |
| FN Atn2#(X#,Y#) | |
| X# | Floating point number -INF<X#<INF |
| Y# | Floating point number -INF<Y#<INF AND Y#<>0 |
| Calculates the arc tangent with two arguments (ARCTAN(X#/Y#)). Depending on the operational sign of X#
and Y#, the quadrant is chosen in such a way as
to place the result into the interval [-PI,PI].
This function is thus very suitable to convert Cartesian coordinates to polar coordinates, for example. Note: The BASIC function ARCTAN(X#/Y#) only returns results in the interval [-PI/2,PI/2] and thus only corresponds with positive X# with FN Atn2#. |
|
| FN Hypot#(X#,Y#) | |
| X# | Floating point number -INF<X#<INF |
| Y# | Floating point number -INF<Y#<INF |
| Calculates SQR(X#^2+Y#^2). The calculation is performed directly, that is without first creating the sum of the squares. This prevents some of the overflows that otherwise occur in the case of a direct evaluation of the formula. | |
| FN Compound#(X#,Y#) | |
| X# | Floating point number 0<X#<INF |
| Y# | Floating point number 0<Y#<INF |
| Calculates (1+X#)^Y#. For small X#, this function returns more accurate values than a direct evaluation of the equation. This function is especially important for financial applications. | |
| FN Annuity#(X#,Y#) | |
| X# | Floating point number 0<X#<INF |
| Y# | Floating point number 0<Y#<INF |
| Calculates (1-(1+X#)^(-Y#)/X#. For small X#, this function returns more accurate values than a direct evaluation of the equation. This function is especially important for financial applications. | |
| FN Exp2#(X#) | |
| X# | Floating point number -INF<X#<INF |
| Calculates 2^X#. | |
| FN Expm1#(X#) | |
| X# | Floating point number -INF<X#<INF |
| Calculates EXP(X#)-1. For small X#, this function returns more accurate values than a direct evaluation of the formula. | |
| FN Lnp1#(X#) | |
| X# | Floating point number -1<X#<INF |
| Calculates LN(X#+1). For small X#, this function returns more accurate values than a direct evaluation of the formula. | |
| FN Erf#(X#) | |
| X# | Floating point number -INF<X#<INF |
| Calculates (2/PI)*INTEGRAL[0,X#]EXP((-t)^2)dt. This integral is also known as the error function and is very significant for statistics and probability calculus. | |
| FN Erfc#(X#) | |
| X# | Floating point number -INF<X#<INF |
| Calculates 1-FN Erf#(X#). This integral is also known as a complementary error function and is very significant for statistics and probability calculus. This function returns more accurate results for large positive number (starting at 10) than a direct evaluation of the formula. | |
| FN Gamma#(X#) | |
| X# | Floating point number -INF<X#<INF AND X#<>0,-1,-2,-3 ... |
| Calculates INTEGRAL[0,INF]EXP(-t)*t^(X#-1)dt. This integral is also known as a gamma function. This function has singularities at zero and all negative whole numbers. The factorial function of Omikron Basic yields the correlation: FN Gamma#(X#)=FACT(X#-1). | |
| FN Lgamma#(X#) | |
| X# | Floating point number -INF<X#<INF AND X#<>0,-1,-2,-3 ... |
| Calculates LN(ABS(FN Gamma#(X#)). Because the gamma function reaches overflow very fast as X# approaches to greater numbers, it is often better to work with the logarithm of the gamma function to avoid any overflows. | |
| FN Hermite#(N,X#) | |
| N | Degree of polynomial 0<=N<INF |
| X# | Floating point number -INF<X#<INF |
| Calculates the value of the Hermite polynomial of degree N at position X#. | |
| FN Cheby#(N,X#) | |
| N | Degree of polynomial 0<=N<INF |
| X# | Floating point number -INF<X#<INF |
| Calculates the value of the Tschebyscheff polynomial of degree N at position X#. Do not confuse this function with the Tschebyscheff approximation, which is discussed further below. | |
| FN Legendre#(N,X#) | |
| N | Degree of polynomial 0<=N<INF |
| X# | Floating point number -INF<X#<INF |
| Calculates the value of the Legendre polynomial of degree N at position X#. | |
| FN Legendre_D#(N,X#) | |
| N | Degree of polynomial 0<=N<INF |
| X# | Floating point number -INF<X#<INF AND X#<>-1 AND X#<>1 |
| Calculates the value of the Legendre polynomial derivation of degree N at position X#. This is also called a Gegenbauer polynomial. | |
| FN Horner#(&A#(),N,X#) | |
| A#(0:N) | Coefficients of the whole rational function. |
| N | Degree of function 0<=N<INF |
| X# | Floating point number -INF<X#<INF |
| Calculates the value of the entirely rational function A#(N)*X#^N+A#(N-1)*X#^(N-1)+...+A#(1)*X#+A#(0) at position X#. |
|
| FN Horner_D#(&A#(),N,X#) | |
| A#(0:N) | Coefficients of the whole rational function. |
| N | Degree of function 0<=N<INF |
| X# | Floating point number -INF<X#<INF |
| Calculates the value of the derivation of the whole rational function A#(N)*X#^N+A#(N-1)*X#^(N-1)+...+A#(1)*X#+A#(0)at position X#. |
|
| FN Beta#(X#,Y#) | |
| X# | Floating point number -INF<X#<INF AND X#<>0,-1,-2,-3 ... |
| Y# | Floating point number -INF<Y#<INF AND X#<>0,-1,-2,-3 ... |
| Calculates the beta function. The following correlation applies: FN Beta#(X#,Y#)=FN Gamma#(X#)*FN Gamma#(Y#)/FN Gamma#(X#+Y#) |
|
| FN Zeta#(X#) | |
| X# | Floating point number -INF<X#<INF |
| Calculates the Riemann zeta function. | |
| FN Kappa#(X#) | |
| X# | Floating point number -INF<X#<INF |
| Calculates the kappa function. | |
| FN Eta#(X#) | |
| X# | Floating point number -INF<X#<INF |
| Calculates the eta function. | |
| FN Li#(N,X#) | |
| N | Order of the poly-logarithm 1<=N<=2 |
| X# | Floating point number -1<=X#<INF |
| Calculates the value of the poly-logarithm of the order N at position X#. | |
| FN Chi#(N,X#) | |
| N | Order of the chi function 1<=N<=2 |
| X# | Floating point number -1<=X#<INF |
| Calculates the value of the chi function of the order N at position X#. | |
| FN Bessel#(N,X#) | |
| N | Order of Bessel function 0<=N |
| X# | Floating point number -INF<=X#<INF |
| Calculates the value of the Bessel function of the order N at position X#. | |
| The folder DEMO contains the program "Bessel.BAS," with which you can depict all Bessel functions graphically. | |
| FN Weber#(N,X#) | |
| N | Order of Weber function 0<=N |
| X# | Floating point number 0<X#<INF |
| Calculates the value of the Weber function of the order N at position X#. Pertinent literature frequently calls the Weber function a Bessel function of the 2nd type. | |
| The folder DEMO contains the program "Bessel.BAS," with which you can depict all Bessel functions graphically. | |
| FN Mod_Bessel#(N,X#) | |
| N | Order of the modified Bessel function 0<=N |
| X# | Floating point number -INF<=X#<INF |
| Calculates the value of the modified Bessel function of the order N at position X#. | |
| The folder DEMO contains the program "Bessel.BAS," with which you can depict all Bessel functions graphically. | |
| FN Macdo#(N,X#) | |
| N | Order of MacDonald function 0<=N |
| X# | Floating point number 0<X#<INF |
| Calculates the value of the MacDonald function of the order N at position X#. Pertinent literature frequently calls the MacDonald function a modified Bessel function of the 2nd type. | |
| The folder DEMO contains the program "Bessel.BAS," with which you can depict all Bessel functions graphically. | |
| 2.3 Zero Point and Extreme Value Calculations |
| Zero points and extreme values are usually very difficult to evaluate in analytical terms. That is the reason we have implemented a few functions into the Numeric Library, which solve such problems using iterative methods. |
| FN Zero#(&FN Fun#(0,0),I,X0#,X1#) | |
| Fun#(I,X#) | This is the function with the zero point to be calculated. Of course, you have to define this function
yourself and then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| I | An index, passed to FN Fun#. |
| X0# | Left limit of interval, in which the zero point is located. |
| X1# | Right limit of interval, in which the zero point is located. |
| This function searches for a zero point of any desired function using the secant method. The function values of X0# and X1# are used as starting values. If no zero point was found after 100 steps, FN Zero# generates an error message and terminates. In this case, the function returns a NAN (Not A Number). | |
| Example: The position at which SIN(3*X#)=COS(5*X#) is supposed to be calculated. We first convert this equation to the form SIN(3*X#)-COS(5*X#)=0. The boundaries of the search interval are X0#=0 and X1#=PI. That results in the following program: Num_Init PRINT FN Zero#(&FN Trig#(0,0),0,0,PI) INPUT "End with [Return]";Dummy Num_Exit END DEF FN Trig#(I,X#)=SIN(3*X#)-COS(5*X#) If you done everything correctly, you should receive the value 2.356... . |
| FN Minimum#(&FN Fun#(0,0),I,E#,S#) | |
| Fun#(I,X#) | This is the function, whose minimum is supposed to be calculated. Of course, you have to define
this function yourself and then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| I | An index, passed to FN Fun#. |
| E# | Estimated value for the position of the minimum. |
| S# | Interval, with which to start the search. A typical value is S#=5, for example. |
| This function searches for a minimum according to the method of parabolic inverse interpolation. Starting point is the passed estimated value E# and the interval S#. If no minimum is found after 200 steps, FN Minimum# generates an error message and terminates. In this case, NAN (Not a Number) is returned as the function value. | |
| Example: The minimum of the function X#^4-10*X#^2+5*X# in the proximity of 2 is supposed to be found: Num_Init PRINT FN Minimum#(&FN Poly#(0,0),0,2,1) INPUT "End with [Return]";Dummy Num_Exit END DEF FN Poly#(I,X#)=X#^4-10*X#^2+5*X# If you done everything correctly, you should receive the value 2.098... . |
| FN Maximum#(&FN Fun#(0,0),I,E#,S#) | |
| Fun#(I,X#) | This is the function, whose maximum is supposed to be calculated. Of course, you have to define
this function yourself and then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| I | An index, passed to FN Fun#. |
| E# | Estimated value for the position of the maximum. |
| S# | Interval, with which to begin the search. A typical value is, for example, S#=5. |
| This function searches for a maximum according to the method of parabolic inverse interpolation. Starting point is the passed estimated value E# and the interval S#. If no minimum is found after 200 steps, FN Maximum# generates an error message and terminates. In this case, NAN (Not A Number) is returned as the function value. | |
| 2.4 Calculation of Eigenvalues and Solving Equation Systems. |
| Equations and equation systems are often difficult to solve in analytical terms. It is either altogether impossible to find an analytical solution or the occurring equations are very complex and hard to handle. We have therefore integrated a few functions into the Numeric Library, which will solve such equations at least approximately. |
| Eigen &M#(,),N,&E#() | |
| M#(0:N,0:N) | The elements of the matrix have to be passed in this field. After the return, this field contains the eigenvectors. If you still need your matrix later on in your program, you have to copy it first into another field using the BASIC command MAT. |
| N | Highest matrix index (dimension minus one). |
| E#(0:N) | After the return, this field contains the eigenvalues. |
| This procedure calculates the eigenvalues of a real symmetrical matrix. If the matrix is not symmetrical or if the algorithm does not converge, you will receive an error message. | |
| Example: The eigenvalues of the matrix are to be calculated, whose elements are indicated in the DATA line: Num_Init -Matrix DATA 2,-1,2,-1,2,-2,2,-2,5 DIM M#(2,2),E#(2) RESTORE Matrix PRINT "Matrix:" FOR J=0 TO 2 FOR I=0 TO 2 READ M#(J,I):PRINT M#(J,I), NEXT I:PRINT NEXT J Eigen &M#(,),2,&E#() PRINT:PRINT "Eigenvalues:" FOR I=0 TO 2:PRINT E#(I):NEXT I INPUT "End with [Return]";Dummy Num_Exit END |
| Gauss &M#(,),&C#(),N,&X#(),R D# | |
| M#(0:N,0:N) | In this field, the coefficients of the equation system have to be passed in form of a matrix. If you still need your matrix later on in your program, you have to copy it first into another field using the BASIC command MAT. |
| C#(0:N) | This field has to contain the constant values on the right side of the equations system. If you still need your field later on in your program, you have to copy it first into another field using the BASIC command MAT. |
| N | Highest field index (number of equations minus one). |
| X#(0:N) | After the return, this field contains the solution vector. |
| D# | The determinant is returned in this variable. |
| This procedure calculates the solution vector X#(), which meets the equation system A#(,)*X#()=C#(). The Gauss process with pivotization is used to keep the computer-related rounding errors to a minimum. | |
| Example: The equation system 3*X0#+4*X1#=11 2*X0#-X1#=0 is supposed to be solved: |
| Num_Init -Equation DATA 3,4,11,2,-1,0 DIM M#(1,1),C#(1),X#(1) RESTORE Equation FOR J=0 TO 1 FOR I=0 TO 1 READ M#(J,I) NEXT I READ C#(J) NEXT J Gauss &M#(,),&C#(),1,&X#(),D# PRINT:PRINT "The solution is:" PRINT "X0#=";X#(0) PRINT "X1#=";X#(1) PRINT "Determinant D#=";D# INPUT "End with [Return]";Dummy Num_Exit END |
| Gauss_Seidel &M#(,),&C#(),N,&X#() | |
| M#(0:N,0:N) | The coefficients of the equation systems have to be passed in this field in the form of matrix. |
| C#(0:N) | This field has to contain the constant values on the right side of the equation system. |
| N | Highest field index (number of equations minus one). |
| X#(0:N) | In this field, you have to pass an approximate solution. After the return, this field contains the improved solution vector. |
| This procedure calculates the solution vector X#(), which meets the equation system A#(,)*X#()=C#(). In case of very large equation systems (100 or more variables), the Gauss algorithm requires a too long
calculation time and the rounding error becomes too great. In these cases, it is better to use an iteration procedure
such as the one by Gauss and Seidel, especially if some approximate values for the solution vector are known. Of course, it is also possible to implement this procedure after the Gauss algorithm and thus reduce rounding errors. Caution: Since the Gauss-Seidel procedure functions iteratively, it is possible that it does not converge. This effect occurs especially if no approximated solution vector is passed. If the algorithm does not converge, you will receive an error message. |
|
| Band &M#(,),N,K,&X#() | |
| M#(0:N,0:2*K+1) | In this field, the coefficients of the equation systems have to be passed. Please note that the constant values from the right side of the equation system have to be passed in the last column. If you still need this field later on in your program, you have to copy it first into another field using the BASIC command MAT. |
| N | Highest field index (number of equations minus one). |
| K | Distance of the element not equal to zero, which is farthest apart from the main diagonal. |
| X#(0:N) | After the return, this field contains the solution vector. |
| This procedure calculates the solution vector X#(), which meets the equation system A#(0:N,0:N)*X#(0:N)=A#(0:N,N+1). Real-life applications often feature equation systems, which have a coefficient matrix that only has elements different of zero in the proximity of the main diagonal. Such matrices are called band-matrices with the width of K. Such an equation system is solved much faster with the band procedure than with the Gauss procedure, because superfluous calculation operations with the zeros of the matrix are avoided. The larger the matrix and the narrower the band in relationship to the dimension of the matrix, the larger the advantage. |
|
| Example: The following equation system is supposed to be solved: 2*X0#+9*X1#=1 2*X0#-3*X1#+12*x2#=2 8*X1#+3*X2#-5*X3#=3 6*X2#+4*X3#+X4#=4 X3#-2*X4#=-5 As one can easily see, only the elements on the main diagonal, as well as the elements on the left and right of it, are not equal to zero. The coefficient matrix can thus be reduced to a band matrix with 3 columns. The constant elements from the right side of the equation system are then also entered into the fourth column. Thus, the following program results with K=1: Num_Init -Band DATA 0,2,9,1 DATA 2,-3,12,2 DATA 8,3,-5,3 DATA 6,4,1,4 DATA 1,-2,0,-5 DIM M#(4,3),X#(4) RESTORE Band FOR J=0 TO 4 FOR I=0 TO 3 READ M#(J,I) NEXT I NEXT J Band &M#(,),4,1,&X#() PRINT:PRINT "The solution is:" FOR I=0 TO 4 PRINT "X"+MID$(STR$(I),2)+"#=";X#(I) NEXT I INPUT "End with [Return]";Dummy Num_Exit END |
| Newton_Sys &FN Fun#(0),&X#(),N | |
| Fun#(I) | This function is used to define the equation system to be solved. Of course, you have to define
this function yourself and then pass its address. In I, your function will receive a parameter, which identifies the number of the equation. The equation has to be transformed in such a way that the right side contains a 0. Your function then has to return the function value from the left side. |
| X#(0:N) | You should use this field to pass an approximated solution, if known. After the return, this field contains the solution vector. |
| N | Highest field index (number of equations minus one). |
| This procedure solves a nonlinear equation system according to the Newton procedure. If no solution was found after 100 steps, Newton_Sys generates an error message and terminates. | |
| Example: The nonlinear equation system X0#*X1#-5*X1#-5=0 X0#-2+X1#*COS(X2#)=0 SIN(X0#+3*X2#)=0 is supposed to be solved: |
| Num_Init DIM X#(2) Newton_Sys &FN Equation#(0),&X#(),2 PRINT:PRINT "The solution is:" FOR I=0 TO 2 PRINT "X"+MID$(STR$(I),2)+"#=";X#(I) NEXT I PRINT:INPUT "End with [Return]";Dummy Num_Exit END DEF FN Equation#(I):LOCAL X# SELECT I CASE 0:X#=X#(0)*X#(1)-5*X#(1)-5 CASE 1:X#=X#(0)-2+X#(1)*COS(X#(2)) CASE 2:X#=SIN(X#(0)+3*X#(2)) END_SELECT RETURN X# |
| 2.5 Interpolation |
| It is often the case that only tabulated values are known of required functions. In such cases,
the obvious thing to do would be to determine the intermediate values by interpolating with a known function. The
Numeric Library offers a variety of procedures and functions, which are suitable for almost all interpolation problems. However, always remember please that an interpolation does not necessarily return the correct intermediate values. For example, functions that cannot be differentiated or functions that are non continuous, cannot even be depicted correctly at all. |
| Lagrange_Int &X#(,),N,&C#() | |
| X#(0:N,0:1) | Use this field to pass the supporting points. X#(0:N,0) contains the X values and X#(0:N,1) the Y values. |
| N | Highest field index (number of supporting points minus one). |
| C#(0:N) | After the return, this field contains the coefficient of the interpolation polynomial. This field has to be passed to FN Lagrange# to calculate the interpolation function. |
| The Lagrange interpolation is a very good procedure, because - unlike others - it also permits non-equidistant values. If very many supporting points exist, the interpolation polynomial begins to oscillate at the interval ends, which causes the interpolation quality to worsen. | |
| The folder DEMO contains the program "Interpolation.BAS," which you may use to interpolate supporting points you have selected yourself by applying different procedures. |
| FN Lagrange# X#,&X#(,),N,&C#() | |
| X# | Here, the function value of the Lagrange interpolation is calculated. Of course, X# has to be located within the interval, in which the function values are known; otherwise, you will obtain nonsensical results. |
| X#(0:N,0:1) | Here, you have to pass the same field with the supporting points as with the procedure Lagrange_Int. |
| N | Highest field index (number of supporting points minus one). |
| C#(0:N) | This field has to contain the coefficients of the interpolation polynomial. It is the exact same field you have received from Lagrange_Int. |
| This function serves to calculate the values of the interpolation polynomial. In practical applications, you would then call the procedure Lagrange_Int first to determine the coefficients of the interpolation polynomial and then calculate the intermediate values using FN Lagrange#. | |
| Spline_Int &X#(,),N,&C#(,) | |
| X#(0:N,0:1) | Use this field to pass the supporting points. X#(0:N,0)contain the X values and X#(0:N,1) the Y values. |
| N | Highest field index (number of supporting points minus one). |
| C#(0:N,0:3) | After the return, this field contains the coefficients of the interpolation polynomial (4 per supporting point).To calculate the interpolation function, this field has to be passed to FN Spline#. |
| This procedure calculates a cubic spline interpolation. By also considering the derivative at the
end of each interpolation interval, the spline interpolation delivers a very smooth interpolation process. The
spline interpolation also permits different intervals between the supporting points. Unlike the Lagrange interpolation, the use of polynomials of the third degree, which are assembled at the supporting points, makes it possible to avoid the strong oscillation of the polynomials when a high number of supporting points exists. The spline interpolation has the characteristic that the tension in the curve is minimized; the curve thus corresponds to an elastic spline, which crosses the supporting points. |
|
| The folder DEMO contains the program "Interpolation.BAS," which you may use to interpolate supporting points you have selected yourself by applying different procedures. |
| FN Spline# X#,&X#(,),N,&C#() | |
| X# | Here, the function value of the spline interpolation is calculated. Of course, X# has to be located within the interval, in which the function values are known; otherwise, you might obtain nonsensical results. |
| X#(0:N,0:1) | Here you have to pass the same field with the supporting points as for the procedure Spline_Int. |
| N | Highest field index (number of supporting points minus one). |
| C#(0:N,0:3) | This field has to contain the coefficients of the interpolation polynomial (4 per supporting point). It is the exact same field, which you have received from Spline_Int. |
| This function serves to calculate the values of the interpolation. In practical applications, you would first call the procedure Spline_Int to determine the coefficients of the interpolation polynomial and then calculate the intermediate values using FN Spline#. | |
| Gauss_Int &FN Fun#(0,0),I,&X#(,),N,&P#(),M | |
| Fun#(I,X#) | This is a function with M freely
selectable parameters, which are in the field P#( ).
Of course, you have to define this function yourself and then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| I | An index, passed to FN Fun#. |
| X#(0:N,0:2) | Use this field to pass the supporting points. X#(0:N,0) contain the X values and X#(0:N,1) the
Y values. In X#(0:N,2), you can still indicate an error tolerance for the individual supporting points. Generally, X#(0:N,2)=1 can be applied. |
| N | Highest field index X#(,) (number of supporting points minus one). |
| P#(0:M) | After the return, this field contains the parameters of your function. |
| M | Highest field index P#() (number of freely selectable parameters of the function FN Fun minus one). |
| This procedure adjusts a function defined by you with M+1 freely selectable parameters to N+1 supporting
points. The method of the smallest error square according to Gauss is used, which means that the sum from the squares
of the errors is minimized at the supporting points. In addition, the Gauss interpolation permits different intervals
between the supporting points. Note: If your function has fewer freely selectable parameters than supporting points are passed, then the function usually does not run through the supporting points. |
|
| The folder DEMO contains the program "Interpolation.BAS," which you may use to interpolate supporting points you have selected yourself by applying different procedures. |
| Fourier_Int &X#(,),N,X1#,X2#,&C#(),M | |
| X#(0:N,0:1) | Use this field to pass the supporting points. X#(0:N,0) contains the X values and X#(0:N,1) the Y values. |
| N | Highest field index X#(,) (number of supporting points minus one). |
| X1# | Left interval limit. |
| X2# | Right interval limit. |
| C#(0:M+M) | After the return, this field contains the coefficients of the trigonometric interpolation function. To calculate the function values, this field has to be passed to FN Fourier#. |
| M | The number of trigonometric terms. |
| This procedure calculates the coefficient of an interpolation function, which consists of cosine
and sine functions. This interpolation type is thus especially suitable for the interpolation of values, which
are periodically repeated, for example, sounds. The interpolation function takes the following form: C#(0)+C#(1)*COS(X#)+C#(2)*COS(2*X#)+ ... +C#(M)*COS(M*X#)+ ... +C#(M+1)*SIN(X#)+C#(M+2)*SIN(2*X#)+C#(M+M)*SIN(M*X#) with X#=2*PI*X#/(X2#-X1#)-PI Caution: The supporting points have to be equidistant. In the case of periodical functions, the interval has to extend over at least one or several whole periods. |
|
| The folder DEMO contains the program "FourierInterpolation.BAS," which demonstrates the use of this interpolation method. |
| FN Fourier# X#,X1#,X2#,&C#(),M | |
| X# | Here, the function value of the Fourier interpolation is being calculated. If a non-periodical series of points has to be interpolated, then X# has to be located within the interval, in which the function values are known; otherwise, you will possibly receive faulty results. |
| X1# | Left interval limit. |
| X2# | Right interval limit. |
| C#(0:M+M) | This field has to contain the coefficients of the trigonometric interpolation function. It is exactly the field, which Fourier_Int has returned to you. |
| M | The number of the trigonometric terms. |
| This function serves to calculate the values of the interpolation. For real-life applications, you would first call the procedure Fourier_Int to determine the coefficients of the trigonometric term and then calculate the intermediate values using FN Fourier#. For periodic sequences of measurements, the interpolation also applies outside of the interpolated interval. | |
| Fourier_Int_Sin &X#(,),N,X1#,X2#,&C#(),M | |
| X#(0:N,0:1) | You use this field to pass the supporting points. X#(0:N,0) contains the X values and X#(0:N,1) the Y values. |
| N | Highest field index X#(,) (number of supporting points minus one). |
| X1# | Left interval limit. |
| X2# | Right interval limit. |
| C#(0:M+M) | After the return, this field contains the coefficients of the trigonometric interpolation function. To calculate the function value, this field has to be passed to FN Fourier_Sin#. |
| M | The number of the sine functions. |
| This procedure calculates the coefficients of an interpolation function, which consists of nothing
but sine functions and a constant term C#(0). This
interpolation type is thus very suitable for the interpolation of values, which repeat themselves periodically.
The interpolation function takes the following form: C#(0)+C#(1)*SIN(X#+C#(M+1))+C#(2)*SIN(2*X#+C#(M+2))+ ... +C#(M)*SIN(M*X#+C#(M+M) with X#=2*PI*X#/(X2#-X1#)-PI Unlike Fourier_Int, the interpolation function contains only half the amount of trigonometric functions; calculations should thus be carried out much speedier. Caution: The supporting points have to be equidistant. In the case of periodical functions, the interval has to extend over at least one or several whole periods. |
|
| The folder DEMO contains the program "FourierInterpolation.BAS," which demonstrates the use of this interpolation method. |
| FN Fourier_Sin# X#,X1#,X2#,&C#(),M | |
| X# | Here, the function value of the Fourier sine interpolation is being calculated. If a non-periodic series of points is to be interpolated, X# has to be located within the interval, in which the function values are known or you might receive false results. |
| X1# | Left interval limit. |
| X2# | Right interval limit. |
| C#(0:M+M) | This field has to contain the coefficients of the trigonometric interpolation function. It is exactly the field, which Fourier_Int_Sin has returned to you. |
| M | The number of the sine functions. |
| This function serves to calculate the values of the interpolation. In real-life applications, you would thus first call the procedure Fourier_Int_Sin to determine the coefficients and phases of the sine functions and then calculate the intermediate values with FN Fourier_Sin#. For periodic sequences of measurements, the interpolation also applies outside of the interpolated interval. | |
| Fft &R#(),&I#(),N,Flag | |
| R#(0:N) | Real parts of the complex values. After the return, this field contains the real parts of the transformed values. If you still need the original values later on in your program, you have to copy them first into another field using the BASIC command MAT. |
| I#(0:N) | Imaginary parts of the complex values. After the return, this field contains the imaginary parts of the transformed values. If you still need the original values in this field later on in your program, you have to copy them first into another field using the BASIC command MAT. |
| N | Number of values. This has to be an even number. Note: The procedure works especially fast if N is a power of two, but also numbers, which can be factored well into 2, 3, 4, 5 result in short calculation times. |
| Flag | Flag = 1 : The normal Fourier transformation is carried
out. Flag = -1 : The inverse Fourier transformation is carried out. |
| This procedure conducts a discrete Fourier transformation. The following sums are basically created: Normal Fourier Transformation: R#(I)=SUM[J=0,N-1]R#(J)*COS(2*PI*I/N)-I#(J)*SIN(2*PI*I/N) I#(I)=SUM[J=0,N-1]I#(J)*COS(2*PI*I/N)+R#(J)*SIN(2*PI*I/N) R#(N)=R#(0):I#(N)=I#(0) Inverse Fourier Transformation: R#(I)=(SUM[J=0,N-1]R#(J)*COS(2*PI*I/N)+I#(J)*SIN(2*PI*I/N))/N I#(I)=(SUM[J=0,N-1]I#(J)*COS(2*PI*I/N)-R#(J)*SIN(2*PI*I/N))/N R#(N)=R#(0):I#(N)=I#(0) However, this procedure also works according to the method of the fast Fourier transform, which shortens the calculation time considerably as if one were to calculate the above sums directly. |
|
| The folder DEMO contains the program "FFT.BAS," which demonstrates the use of this procedure. |
| 2.6 Approximation of Functions Through Function Systems |
| One often has to face the problem of a function being known, but its evaluation being rather difficult. It is then advantageous to approximate the result using other, simpler functions. The most often applied functions of this type are included in the Numeric Library. |
| Cheby_App &FN Fun#(0,0),I,X1#,X2#,&C#(),N | |
| Fun#(I,X#) | This is the function to be approximated. Of course, you have to define this function yourself and then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| I | An index, passed to FN Fun#. |
| X1# | Left interval limit. |
| X2# | Right interval limit. |
| C#(0:N) | After the return, this field contains the coefficients for the Chebyshev approximation. To calculate the function values, this field has to be passed to FN Cheby#. |
| N | Highest field index C#() (number of Chebyshev polynomials minus one). |
| This procedure adjusts an analytical function FN Fun# defined by you to the orthogonal function system from the Chebyshev polynomials. The advantage of the Chebyshev approximation consists of normally needing only few coefficients to achieve a rather close approximation and that a very efficient method exists to calculate the Chebyshev polynomials rather rapidly. For this, you may use the function FN Cheby#. Note: This procedure contrasts with the procedure Mean_App, discussed further below, because no integrals have to be evaluated to calculate the coefficients, which means that larger numbers are also possible for the number of polynomials (N > 100). The calculation is indeed a bit longer, but the results are well-worth the waiting time. |
| The folder DEMO contains the program "Approximation.BAS," which demonstrates the use of this procedure. |
| FN Cheby# X#,X1#,X2#,&C#(),N | |
| X# | Here, the function value of the Chebyshev approximation is being calculated. Of course, X# has to be located within the interval, in which the function values are known or you might receive false results. |
| X1# | Left interval limit. |
| X2# | Right interval limit. |
| C#(0:N) | This field has to contain the coefficients of the Chebyshev approximation function. It is exactly the field, which Cheby_App has returned to you. |
| N | Highest field index C#() (number of Chebyshev polynomials minus one). |
| This function serves to calculate the values of the approximation. For real-life applications, you
first call the procedure Cheby_App to determine
the coefficients of the Chebyshev polynomials and then calculate the approximation values using FN Cheby#. Do not confuse this function with the one discussed above (FN Cheby#(I,X#)), which only serves to calculate the function value of a very specific Chebyshev polynomial. |
|
| The folder DEMO contains the program "Approximation.BAS," which demonstrates the use of this procedure. |
| Mean_App &FN Fun#(0,0),I,X1#,X2#,&FN Sys#(0,0),&C#(),N | |
| Fun#(I,X#) | This is the function to be approximated. Of course, you have to define this function yourself and then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| I | An index, passed to FN Fun#. |
| X1# | Left interval limit. |
| X2# | Right interval limit. |
| Sys#(I,X#) | This is the function systems, with which to approximate FN
Fun#. It has to consist of N+1
elements and use the coefficients C#(0:N). Of course, you have to define this function system yourself and then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| C#(0:N) | After the return, this field contains the coefficients for your function systems. |
| N | Highest field index C#() (number of functions defined by you minus one). |
| This procedure adapts an analytical function FN Fun# defined by you to the function system FN Sys# also
defined by you in such a way that the error square integral is minimized. The adjusted function calculation is
as follows: C#(0)*FN Sys#(0,X#)+C#(1)*FN Sys#(1,X#)+ ... +C#(N)*FN Sys#(N,X#) Note: Since many integrals have to be evaluated to calculate the coefficients, the computation time might be quite long in the case of larger function systems (N >> 10) . In addition, number overflows may occur, which falsify the result completely. It is therefore recommended to use the control word COMPILER "DEBUG ON" - at least during the test phase, even if the calculation time might increase a bit. |
| Example: To approximate the Euler function (EXP(X#)) using a power series: |
| Num_Init N=5 DIM C#(N) Mean_App &FN Euler#(0,0),0,-1,1,&FN Potenz_Sys#(0,0),&C#(),N PRINT:PRINT " X#";TAB(7);"Initial function"; PRINT TAB(27);"Approximation" FOR X#=-1 TO 1 STEP 0.1 PRINT USING "##.##";X#;TAB(6);USING "";EXP(X#); PRINT TAB(25);FN Potenz_App#(X#) NEXT X# PRINT:INPUT "End with [Return]";Dummy Num_Exit END DEF FN Euler#(I,X#)=EXP(X#) DEF FN Potenz_Sys#(I,X#)=X#^I DEF FN Potenz_App#(X#):LOCAL I,Y#=0 FOR I=0 TO N Y#+=C#(I)*X#^I NEXT I RETURN Y# |
| As you can try out easily with this program, N=5 yields better approximation values than N=10.
At N=15, already false results occur, because rounding
errors accumulate during the calculation of the coefficients and at N=20,
an underflow is already occurring. The folder DEMO contains the program "Approximation.BAS," which also shows the use of this procedure. |
| 2.7 Numerical Differentiation | ||||||||||||
| This chapter deals with the derivative of functions. The numerical differentiation is very sensitive when it comes to rounding errors. The results are thus always only approximation values. | ||||||||||||
|
||||||||||||
| The folder DEMO contains the program "Derivative.BAS," which shows the use of this function. |
| Derivative &X#(,),N | |
| X#(0:N,0:3) | In this field, you pass the function values. X#(0:N,0) contains the X values and X#(0:N,1) the
Y values. After the return, X#(0:N,1) contains the 1st derivative, X#(0:N,2) the 2nd derivative, and X#(0:N,3) the 3rd derivative at the position X#(0:N,0). |
| N | Highest field index X#(,) (number of the function values minus one). |
| This procedure calculates the first three derivatives of a tabulated function. Please note that the quality of the results very strongly depends on the quality of the data entered in X#(0:N,0:1). Here, the old adage "Garbage in, garbage out" applies. | |
| Example: The first three derivatives of the function X#^3 are to be calculated: Num_Init N=20:I=0 DIM X#(N,4) FOR X#=-1 TO 1 STEP 0.1 X#(I,0)=X# X#(I,1)=X#^3 I+=1 NEXT X# Derivative &X#(,),N PRINT:PRINT " X#";TAB(7);"function value"; PRINT TAB(27);"1st derivative";TAB(47);"2nd derivative"; PRINT TAB(67);"3rd derivative" I=0 FOR X#=-1 TO 1 STEP 0.1 PRINT USING "##.##";X#;TAB(7);USING "##.#########";X#(I,1); PRINT TAB(27);USING "##.#########";X#(I,2);TAB(47); PRINT TAB(27);USING "##.#########";X#(I,3);TAB(67);X#(I,4) I+=1 NEXT X# PRINT:INPUT "End with [Return]";Dummy Num_Exit END |
| Since the function values stem from a smooth analytical function, the results are in this case rather good as well as you may see easily by calculating yourself. |
| 2.8. Numerical Integration |
| Unlike differentiation, the integration of functions does not yield any generally valid algorithms
any more, which always permit to indicate the value of the integrals in the form of a closed expression. This means
that already functions such as SIN(X#)/X#, which
appear relatively easy, cannot be integrated in an elementary way. This increases the importance of having numerical methods at one's disposal, which can be used to calculate such integrals with a good level of accuracy. |
| FN Newton_Cotes#(&X#(,),N) | |
| X#(0:N,0:1) | In this field, you pass the function values. X#(0:N,0) contains the X values and X#(0:N,1) the
Y values. Caution: The X values have to be equidistant. |
| N | Highest field index X#(,) (number of function values minus one). N has to be larger or equal 5; otherwise you will receive an error message. |
| This function calculates the integral of the discrete function values passed in the field X#(,). Note: If your values are not equidistant, you can initially also use one of the previously discussed interpolation functions to obtain equidistant values and then pass these to FN Newton_Cotes#. An alternative is to use one of the integration functions discussed further below, to which you then pass the address of the interpolation function. |
|
| Example: A measurement resulted in the measurement values indicated in the DATA line. Now, the integral, that is the area under the measurement curve is to be calculated. Since the X values are not equidistant, we will first interpolate: |
| Num_Init -Values DATA 0,1,0.5,1.649,1,2.718,2,1.4,4,0.5,5,2 N=5 DIM X#(N,1),Y#(N,1),C#(N) RESTORE Values FOR I=0 TO 5:READ X#(I,0),X#(I,1):NEXT I 'First, generate equidistant measurement 'values through interpolation. Lagrange_Int &X#(,),N,&C#() FOR I=0 TO 5 Y#(I,0)=I Y#(I,1)=FN Lagrange#(I,&X#(,),N,&C#()) NEXT I PRINT "The integral has the value: "; PRINT FN Newton_Cotes#(&Y#(,),N) PRINT:INPUT "End with [Return]";Dummy Num_Exit END |
| FN Gauss_Leg#(&FN Fun#(0,0),I,X1#,X2#,N) | |
| Fun#(I,X#) | This is the function to be integrated. Of course, you have to define this function yourself and
then pass its address. In I, your function will receive a parameter, which you may use for your own purposes, for example, to identify different variations of a function. In X#, the argument is passed, whose function value your function has to return. |
| I | An index, passed to FN Fun#. |
| X1# | Left interval limit. |