Numeric Library / Reference 
Omikron Basic on the Internet: http://www.berkhan.de 
General Information turn pageTable of Contents 
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 practiceoriented 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 1FN 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#(N1)*X#^(N1)+...+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#(N1)*X#^(N1)+...+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 polylogarithm 1<=N<=2 
X#  Floating point number 1<=X#<INF 
Calculates the value of the polylogarithm 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#^410*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#^410*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 computerrelated 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 GaussSeidel 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). Reallife 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 bandmatrices 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 nonequidistant 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 nonperiodical 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 reallife 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 nonperiodic 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 reallife 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,N1]R#(J)*COS(2*PI*I/N)I#(J)*SIN(2*PI*I/N) I#(I)=SUM[J=0,N1]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,N1]R#(J)*COS(2*PI*I/N)+I#(J)*SIN(2*PI*I/N))/N I#(I)=(SUM[J=0,N1]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 wellworth 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 reallife 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. 
X2#  Right interval limit. 
N  This parameter may be used to influence the accuracy. N can have values between 2 and 15. The higher the value for N, the more accurate is the result, but the slower is also the calculation time. 
The integration according to GaussLegendre is one of the fastest methods to calculate integrals numerically. It functions only well, however, if relatively harmless functions such as SIN(X#) or COS(X#) with a limited interval are involved. Calculating the integral 1/X# in the interval [1,1000] already leads to completely false results even for N=15. Therefore, it is better to use the function FN Romberg# discussed further below in such cases. 
Example: The integral of SIN(X#) in the interval [0,PI] is to be calculated: Num_Init PRINT "The integral has the value: "; PRINT FN Gauss_Leg#(&FN Sinus#(0,0),0,0,PI,5) PRINT:INPUT "End with [Return]";Dummy Num_Exit END DEF FN Sinus#(I,X#)=SIN(X#) The result of 2.00000011020611 is rather close to the theoretical result of 2, although only 5 has been chosen for N. Try out the calculation for yourself using other values for N and other interval limits. The folder DEMO contains the program "Integral.BAS," which compares different integration processes with one another. 
FN Romberg#(&FN Fun#(0,0),I,X1#,X2#)  
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. 
X2#  Right interval limit. 
This function calculates the integral according to the method by Romberg. The calculation speed
is significantly lower than with FN Gauss_Leg#,
in return, the results are rather accurate even for problematic functions. This is achieved with a recursive algorithm
called by the function itself until a relative error falls below 1D10. This requires the memory on the BASIC stack,
which you might have to enlarge using COMPILER "STACK X." In spite of the good convergence of FN Romberg#, you cannot just integrate passing over singularities, of course. You have to make sure that your function within the entire integration interval takes on only values within the depictable number range minus a safety interval. If the calculation does not converge, you will receive an error message. 
Please also note the sample program Romberg.BAS in the folder DEMO. It calculates the integral sine,
the integral cosine, the integral exponential function, and the integral logarithm. All 4 integrals cannot be indicated
in a closed form and the integrands of the last 3 also have singularities, which have to be replaced by finite
values. The values of 100 and/or 100 are the supposed integration limits for the second and third functions. Theoretically, of course, this would have to be infinite and/or minus infinite. However, since the integrand strives towards 0, one only receives rounding errors, which falsify the entire result if the selected interval is too large. 
FN Romberg_2#(&FN Fun1#(0,0),I1,&FN Fun2#(0,0),I2,X1#,X2#)  
Fun1#(I1,X#) Fun2#(I2,X#) 
These are the two functions, whose product is to be integrated. Of course, you have to define this
functions yourself and then pass its addresses. In I1 and I2, your functions will receive parameters, 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 functions have to return. 
I1,I2  An index, passed to FN Fun1# and/or FN Fun2#. 
X1#  Left interval limit. 
X2#  Right interval limit. 
Many tasks, for example, developing functions according to function systems, result in the problem
that one has to integrate not the individual function but rather the product of different functions. For that reason,
we have included FN Romberg_2# in the Numeric
Library. The integration takes place according to the Romberg method. The results are rather accurate even for
problematic functions. This is achieved with a recursive algorithm calling the function itself until a relative
error falls below 1D10. This requires the memory on the BASIC stack, which you might have to enlarge using COMPILER
"STACK X." In spite of the good convergence of FN Romberg_2#, you cannot just integrate passing over singularities, of course. You have to make sure that your function within the entire integration interval takes on only values within the depictable number range minus a safety interval. If the calculation does not converge, you will receive an error message. 
2.9 Solving Differential Equations 
Solving differential equations and differential equation systems is a very important aspect within numerical mathematics, because the solution functions cannot be found in any other way in many cases and often cannot be depicted in a closed form either. Therefore, we have integrated the most often used methods into the Numeric Library. 
Runge_Kutta &FN Fun#(0,0,0),I,Xs#,Xe#,R Y#,R S#,E#  
Fun#(I,X#,Y#)  This function defines the differential equation using the correlation Y#'=FN
Fun#(I,X#,Y#). 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. 
I  An index, passed to FN Fun#. 
Xs#  Initial value for X#. 
Xe#  Ending value for X#. 
Y#  Before the invocation, you have to pass the initial value for Y#. After the return, you will receive the ending value of Y# in this variable. 
S#  This parameter has to be used to pass a recommendation for the first interval (S#>=0.01). After the return, this contains the variable of the interval used last. 
E#  This serves to stipulate the maximum size of the local error (1D14 <= E# <= 1D4). The local error is to be understood as the processimplied error, which is still permissible for every integration step. This has the result that the value E# has a direct effect on the stepsize control. 
This procedure solves the initial value problem of a differential equation of the first order according to the RungeKutta procedure. The procedure takes over the complete interval control in this case. 
The folder DEMO contains the program "DiffEquation.BAS," which compares the RungeKutta process with the AdamBashforth process with the help of an example. 
FN Runge_Step#(&FN Fun#(0,0,0),I,X#,Y#,S#)  
Fun#(I,X#,Y#)  This function defines the differential equation using the correlation Y#'=FN
Fun#(I,X#,Y#). 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. 
I  An index, passed to FN Fun#. 
X#  X value for which the step is to be executed. 
Y#  Yvalue for which the step is to be executed. 
S#  This parameter has to be used to pass the stepsize. 
This function executes a single step according to the RungeKutta procedure. You may use this function for an individual interval control program or simply just use the procedure Runge_Kutta right away. 
Adam_Bash &FN Fun#(0,0,0),I,Xs#,Xe#,Y#,&X#(,),N  
Fun#(I,X#,Y#)  This function defines the differential equation using the correlation Y#'=FN
Fun#(I,X#,Y#). 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. 
I  An index, passed to FN Fun#. 
Xs#  Initial value for X#. 
Xe#  Ending value for X#. 
Y#  Before the invocation, you have to pass the initial value for Y#. 
X#(0:N,0:1)  In this field, you will receive the solution values of your differential equation at the equidistant 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.). 
This procedure solves the initial value problem of a differential equation of the first order according to the AdamsBashforth four step procedure. The procedure works a bit faster than Runge_Kutta, however, it functions only with a fixed interval. It is therefore especially suitable if you want to receive the solution of a differential equation between initial and ending value of equidistant supporting points. 
The folder DEMO contains the program "DiffEquation.BAS," which compares the RungeKutta process with the AdamBashforth process with the help of an example. 
FN Adam_Step#(&FN Fun#(0,0,0),I,X#,Y#,Y0#,Y1#,Y2#,Y3#,S#)  
Fun#(I,X#,Y#)  This function defines the differential equation using the correlation Y#'=FN
Fun#(I,X#,Y#). 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. 
I  An index, passed to FN Fun#. 
X#  X value for which the step is to be executed. 
Y#  Y value for which the step is to be executed. 
Y0#,Y1#, Y2#,Y3# 
Here you have to pass the solution values of your differential equation at the last four supporting points. 
S#  This parameter has to be used to pass the stepsize. 
This function executes a single step according to the AdamsBashforth four step process. You may use this function for an individual solution program or just use the procedure Adam_Bash right away. 
Runge_Step_Sys &FN Fun#(0,0,0,0),I,X#,&Y#(),N,S#  
Fun#(I,X#,Y,Yd)  This function defines the differential equation system. 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. X# is the independent variable. Y is the address of a floating point field, which represents the independent variables. Yd is the address of a floating point field, in which your function has to return the first derivative of the independent variable. 
I  An index, passed to FN Fun#. 
X#  X value for which the step is to be executed. 
Y#(0:N)  This field has to be used to pass the initial values of the dependent variables. 
N  Highest field index Y#() (number of differential equations minus one). 
S#  This parameter has to be used to pass the stepsize. 
This procedure executes one single step according to the RungeKutta procedure for differential equation systems. You may use this procedure for an individual solution program and then work with a constant interval or use an interval control adapted to the equation system. 
The folder DEMO contains the program "RungeKuttaSys.BAS," which uses the RungeKutta process for differential equation systems with the help of an example and also considers the impact of different stepsizes. 
De &FN Fun#(0,0,0,0),I,R Xs#,Xe#,&Y#(),N,R Re#,R Ae#,R F  
Fun#(I,X#,Y,Yd)  This function defines the differential equation system. 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. X# is the independent variable. Y is the address of a floating point field, which represents the independent variables. Yd is the address of a floating point field, in which your function has to return the first derivative of the independent variable. 
I  An index, passed to FN Fun#. 
Xs#,Xe#  These two parameters define the interval, that is supposed to be used during integration. Xs# may also be larger than Xe#, because the procedure can also integrate backwards. After the return, Xs# contains the last used value. This information may be important for you, in case the integration has to be terminated. 
Y#(0:N)  In this field, the initial values of the dependent variable have to be passed. After the return, this field contains the solution vector at the position Xe#. 
N  Highest field index Y#() (number of differential equations minus one). 
Re#,Ae#  This parameters are used to determine the permissible relative error (Re#) and the absolute error (Ae#). The procedure
uses the following inequation: ABS(Absolute_Error#)<=Re#*ABS(Y#(0:N))+Ae# By setting Re# and/or Ae# at zero, it is possible to obtain a pure absolute as well as a pure relative error criterion. However, this rarely has any practical applications, because most of the differential equations offer little information about the behavior of the solutions. This is the reason that it is usually useful to define a relative as well as an absolute error. If the solution then exhibits a rather large fluctuation concerning its absolute value, then Re# is mainly used for controlling the integrator; in case of a small solution concerning the absolute value, Ae# is used. This means you are always on the safe side. 
F  This variable is used by the procedure to communicate with your program. The detailed significance
of each variable from F (1 or 1 to 6) is as follows: F=1,1: This value has to be set by the user, if he or she wants to restart the integration, that is if an integration is to be started anew or if the differential equation possesses, for example, a location with discontinuity, which cannot be ignored when integrating. In the case of the discontinuity, it also makes sense to start with F=1. This commands the integrator not to integrate beyond Xe#, thus the discontinuity is not exceeded. For reasons of the optimal interval choice and in case of positive values of F, the algorithm does not try to reach Xe# exactly. Rather, the solution vector at position Xe# of the system is normally determined by interpolation. This, however, does not result in any loss of accuracy. During the integration processes, F can assume different values, which inform your program about the progress of the integration. You may not change these values. 
F  F=2: This value signals a successful conclusion of the integration. Xs# equals Xe#, so that you only have to set Xe# to a new value before the next call to continue the integration. For that, all internal fields of the integrator are already initialized. In Y#(0:N), the solution vector of the differential equation system at position Xe# is located. F=3: Re# and Ae# are too small for the existing calculation accuracy (16 decimal points). Xs# is set to the value, which was nearest to Xe# during the integration and Y#(0:N) contains the solution vector Here, Re# and Ae# were set from De to acceptable values, so that the integrator only has to be called again to continue the integration with these higher error limits. 
F  F=4: The integration was terminated, because more than 500 internal steps (that is 1000 evaluations of the differential equation system) are required to reach Xe#. Xs# is set to the value, which was closest to Xe# during the integration. Your program can now decide whether it wants to continue the integration in spite of the long computation time or if it wants to terminate the integration process. To continue the integration, all you need to do is call De again. F=5: This has the same significance as F=4; however, it was additionally determined that the differential equations seem to be stiff. This can lead to longer calculation times, even if the results are still rather good. To continue the integration, all you need to do is call De again. 
F  F=6: Every time De is called, a test checks whether the entered parameters are permitted. Accordingly, F=6 is set, if Xs#=Xe#, Re#, or Ae# is smaller than zero or if the absolute value of F is not between 1 and 5. 
This procedure is a very frequently used integrator to solve differential equation systems of the
first order. Originally, this was a FORTRAN program, which is here made available in the form of an Omikron Basic
variation. The integration is carried out according to the Adams procedure of variable orders and variable stepsizes. The procedure controls the order and stepsizes independently, while it tries to keep the processimplied error below a relative and absolute error defined by the user. The order can vary between 1 and 13, that is the algorithm can still exactly integrate differential equations, whose solutions are polynomials of a maximal 13th order (disregarding the rounding errors in this case). In addition, a stiffness of the differential equations can be recognized in most cases and the integration is then terminated. This procedure is slower than the ones discussed further up, however, it generally delivers better results. If you have to deal with the numerical solution of differential equations quite frequently, you should have a look at the book by Shampine & Gordon. They have discussed the derivation and use of the procedure De in the FORTRAN variation in detail. Note: Problems can occur during the integration process, if one of the differential equations will be exactly equal to zero during the course of the integration, because this then does not define the relative error. However, this rarely if ever occurs due to the internal calculation errors, except if the differential equations of the initial value are set to zero right from the start. In this case, you can manage the problem by setting the respective values to very low numbers (e.g., 1D10), instead of setting them to zero. This results in practically no additional calculation errors, and the situation is avoided. It may also serve a purpose for critical cases to set the control word COMPILER "DEBUG ON" to be notified of nonpermitted floating point operations (e.g., division by 0). 
We have used rather transparent differential equations for the following examples, which also have analytical solutions with which to compare the numerical results: 
Example 1: Num_Init 'Solution of Y'=Y. DIM Y#(0),Yd#(0) 'Initial values. Xs#=0:Xe#=10 Y#(0)=1 Re#=1D5:Ae#=0 F=1 De &FN Diff_Equation#(0,0,0,0),0,Xs#,Xe#,&Y#(),0,Re#,Ae#,F PRINT "The numerical solution is:.. "; PRINT "X=";Xs#,"Y=";Y#(0) PRINT "The theoretical solution is: "; PRINT "X=";Xs#,"Y=";EXP(10) PRINT:INPUT "End with [Return]";Dummy Num_Exit END DEF FN Diff_Equation#(I,X#,Y,Yd) *Yd#(0)=*Y#(0) RETURN 0 
Example 2: The circular motions of a body without mass around a mass center is described by the differential equations Y0''=Y0*R^3 and Y1''=Y1*R^3. This is a system of the second degree, which, as is well known, can always be converted into a system of the first degree. In this case, four differential equations of the first degree are produced: 
Num_Init 'Solving a celestrialmechanical problem. DIM Y#(3),Yd#(3) 'Initial values. Xs#=0:Xe#=2*PI:'Corresponds to a rotation. Y#(0)=1:Y#(1)=0:'Initial position. Y#(2)=0:Y#(3)=1:'Initial velocity. Re#=1D6:Ae#=1D6 F=1 REPEAT De &FN Orbit#(0,0,0,0),0,Xs#,Xe#,&Y#(),3,Re#,Ae#,F UNTIL F=2 or F=6 If F=6 THEN PRINT "Integration was terminated":GOTO Ende PRINT "Numerical solution:","Theoretical solution:" PRINT Y#(0),COS(Xe#) PRINT Y#(1),SIN(Xe#) PRINT Y#(2),SIN(Xe#) PRINT Y#(3),COS(Xe#) Ende:PRINT:INPUT "End with [Return]";Dummy Num_Exit END DEF FN Orbit#(I,X#,Y,Yd) LOCAL R#=(*Y#(0)^2+*Y#(1)^2)^1.5 'Clearance of the bodies. *Yd#(0)=*Y#(2):'Velocity is the *Yd#(1)=*Y#(3):'derivative of the location. *Yd#(2)=*Y#(0)/R#:'Acceleration is the *Yd#(3)=*Y#(1)/R#:'derivative of the velocity. RETURN 0 The folder DEMO contains the program "Riccati.BAS," which uses another example to illustrate the many application possibilities of this procedure. 
2.10 Graphical Output 
The graphical depiction of functions or measured values is most of the time much more expressive than a listing in a table. That is the reason the Numeric Library contains a few procedures, which draw the functions or individual values into a coordinate system. We have tried our best to keep these procedures as flexible as possible to retain their applicability for as many cases as possible. However, predefined routines always represent a compromise between ease of use and the highest level of applicability. 
Plot_Coordinates X,Y,W,H,X1#,X2#,Y1#,Y2#,F  
X  X value of the upper left corner of the coordinate system in pixels. 
Y  Y value of the upper left corner of the coordinate system in pixels. 
W  Width of the coordinate system in pixels. 
H  Height of the coordinate system in pixels. 
X1#  Minimal X# value. 
X2#  Maximal X# value. 
Y1#  Minimal Y# value. 
Y2#  Maximal Y# value. 
F  This parameter may be used to stipulate which element is to be drawn. Each bit of this integer number
represents a certain element, which is only drawn if the corresponding bit is set. Bits have the following significance: Bit 0 : Fill the area of the coordinate system. Bit 1 : Draw frame around coordinate system. Bit 2 : Draw coordinate crosshair. Bit 3 : Label coordinate axes. Bit 4 : Draw auxiliary lines. 
This procedure draws a coordinate crosshair with different attributes. The drawing elements use
the current settings for line width, line color, fill color, etc. This means that you have to make sure yourself
that these attributes are set suitably. But this also gives you the opportunity to influence the look of your coordinate
crosshair directly. If, for example, you would like to draw the coordinate axes in a thicker line width than the auxiliary lines, first set the line width to 2, then draw the coordinate system without auxiliary lines, adjust the line width to 1, and again call the procedure, while setting the flags so that only the auxiliary lines are being drawn. Caution: This procedure resets the clipping rectangle. If you would like to use a clipping rectangle yourself, you have to recreate it after the invocation. 
The folder DEMO contains a variety of programs, which use this procedure. 
Plot_Function &FN Fun#(0,0),I,X,Y,W,H,X1#,X2#,Y1#,Y2#,F[,S]  
Fun#(I,X#)  This is the function, whose graph is to be drawn. 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  X value of the upper left corner of the coordinate system in pixels. 
Y  Y value of the upper left corner of the coordinate system in pixels. 
W  Width of the coordinate system in pixels. 
H  Height of the coordinate system in pixels. 
X1#  Minimal X value. 
X2#  Maximal X value. 
Y1#  Minimal Y value. 
Y2#  Maximal Y value. 
F  This parameter may be used to stipulate which element is to be drawn. Each bit of this integer number
represents a certain element, which is only drawn if the corresponding bit is set. Bits have the following significance: Bit 0 : Fill the area of the coordinate system. Bit 1 : Draw frame around coordinate system. Bit 2 : Draw coordinate crosshair. Bit 3 : Label coordinate axes. Bit 4 : Draw auxiliary lines. 
S  This parameter may be used to determine at what interval of how many pixels in the direction of
X a function value is being calculated. If you omit this parameter or pass a value of <=1, S=1 will be set, which means that one function value per pixel is being
calculated. In those cases where the calculation of the function values is rather laborious, it might make sense to select a higher value for S. For example, see the program Riccati.BAS in the DEMO folder. 
This procedure draws the graph of a function. A coordinate crosshair may be drawn simultaneously.
The drawing elements use the current settings for line width, line color, fill color, etc. This means that you
have to make sure yourself that these attributes are set suitably. If you would like your graphical depiction to be a bit fancier, first draw only the coordinate system using Plot_Coordinates , then draw with a different line color, for example, the actual graph, while the attributes are turned off (F=0). Caution: This procedure resets the clipping rectangle. If you would like to use a clipping rectangle yourself, you have to recreate it after the invocation. 
The folder DEMO contains a variety of programs, which use this procedure. 
Plot_Array &X#(0,0),N,X,Y,W,H,X1#,X2#,Y1#,Y2#,F  
X#(0:N,0:1)  In this field, you pass the points to be drawn. X#(0:N,0) contains the X values and X#(0:N,1) the Y values. 
N  Highest field index X#(,) (number of points minus one). 
X  X value of the upper left corner of the coordinate system in pixels. 
Y  Y value of the upper left corner of the coordinate system in pixels. 
W  Width of the coordinate system in pixels. 
H  Height of the coordinate system in pixels. 
X1#  Minimal X value. 
X2#  Maximal X value. 
Y1#  Minimal Y value. 
Y2#  Maximal Y value. 
F  This parameter may be used to stipulate which element is to be drawn. Each bit of this integer number
represents a certain element, which is only drawn if the corresponding bit is set. Bits have the following significance: Bit 0 : Fill the area of the coordinate system. Bit 1 : Draw frame around coordinate system. Bit 2 : Draw coordinate crosshair. Bit 3 : Label coordinate axes. Bit 4 : Draw auxiliary lines. Bits 8  10 determine whether the individual points are connected by lines or are to be drawn as individual, geometrical objects. Bit 10 9 8 0 0 0 : Connect points with lines. 0 0 1 : Xshaped crosses. 0 1 0 : Normal crosses. 0 1 1 : Triangles located on base line. 1 0 0 : Inverted triangles (resting on points). 1 0 1 : Squares resting on points. 1 1 0 : Squares resting on the side lines. 1 1 1 : Circles. The option to define different types of depictions for points is an advantage, especially if the coordinate cross is to contain several curves but a black and white print out, for example, does not have any other colors available to distinguish between the different lines. It is then possible to make that distinction by using differentlyshaped points. Tip: Filled objects may also be generated simply by enlarging the line width. 
This procedure draws discrete points (e.g., measurement values), which are either connected with
lines or depicted as geometrical objects. A coordinate crosshair can be drawn simultaneously. The drawing elements
use the current settings for line width, line color, fill color, etc. This means that you have to make sure yourself
that these attributes are set suitably. If you would like your graphical depiction to be a bit fancier, first draw only the coordinate system using Plot_Coordinates , then draw with a different line color, for example, the actual graph, while the attributes are turned off (F=0). Caution: This procedure resets the clipping rectangle. If you would like to use a clipping rectangle yourself, you have to recreate it after the invocation. 
The folder DEMO contains a variety of programs, which use this procedure. 
General Information turn pageTable of Contents 
Omikron Basic on the Internet: http://www.berkhan.de 
Copyright 19981999 by BerkhanSoftware