PO 21 Customized Base SAS implemented solvers

Solvers are software tools or subroutines, designed for solving a problem in the area of numerical analysis. The use of efficient reusable mathematical software as the building blocks for complex programs is beneficial for creation of viable applications. Presently, there are numerous scientific software libraries created on the basis of known high-level programming languages. SAS also has a large number of solvers available in SAS/OR and SAS/IML modules, most of which are not available to SAS Base only users. Thus, the creation of Base SAS implemented scientific software repository of efficient customized solvers is of great importance for SAS community. The goal of this paper is a presentation of three customized Base SAS implemented solvers, based on well known mathematical algorithms. Particularly, the first solver finds a solution x to the equation 0 ) ( = x f using the bisection method as a root-finding algorithm. The second solver is intended for calculating of the numerical value of a one-dimensional definite integral dx x f b a ∫ ) ( , using the modified adaptive Simpson’s method based on polynomial of order 2 with the adaptive increment. The third solver finds the roots of a system of linear equations{ j j ij b x a = , where n j i ,..., 1 , = using Cramer’s rule and Leibniz formula. We hope that these solvers might be considered as prototypes for future enhancements, and spur other SAS users and developers to contribute to the future Base SAS solvers repository. INTRODUCTION There is a wide range of mathematical routines, solvers, subroutines, designed for solving a problem in the area of numerical analysis in applied mathematics. The use of efficient reusable mathematical software as the building blocks for complex programs is beneficial for creation of viable applications. Presently, there are numerous scientific software libraries created on the basis of known high-level programming languages (see, e.g., [1]). SAS also has a large number of solvers available in SAS/OR and SAS/IML modules, however, most of which are not available to SAS Base only users. Thus, the creation of Base SAS implemented scientific software repository of efficient customized solvers is of great importance for SAS community. The goal of this paper is a demonstration of using SAS language capabilities for creation of three customized Base SAS implemented solvers, based on well known mathematical algorithms. Particularly, the first solver finds a solution x to the equation 0 ) ( = x f using the bisection method as a root-finding algorithm. The second solver is intended for calculating of the numerical value of a one-dimensional definite integral dx x f b a ∫ ) ( , using the modified adaptive Simpson’s method based on polynomial of order 2 with the adaptive increment. The third solver finds the roots of a system of linear equations { j j ij b x a = , where n j i ,..., 1 , = using Cramer’s rule and Leibniz formula. We hope that these solvers might be considered as prototypes for future enhancements, and spur other SAS users and developers to contribute to the future Base SAS solvers repository. THE SOLVER FOR FINDING SOLUTION X TO THE EQUATION 0 ) ( = x f The first solver presented to your attention handles one of the most important topics in numerical analysis, namely finding a solution x to the equation 0 ) ( = x f . From many root-finding algorithms developed for these methods (e.g., bisection method, Newton’s method, the secant method, etc.) we will present SAS implemented code for the bisection method (see, e.g., [2]), which is based on dividing an interval in half and then selecting the subinterval in which the root exists. It works when the function is a continuous function, and it requires previous knowledge of 2 initial guesses, a and b , such that ) (a f and ) (b f have opposite signs. The bisection method is one of the bracketing methods for finding roots of equation since it keeps a 2 bounded interval where there is at least one root at each step. Given a function ) (x f and an interval which might contain a root, a predetermined number of iterations are carried out. A simple bisection procedure for iteratively converging on a solution which is known to lie inside some interval [ ] b a, proceeds by evaluating the function in question at the midpoint of the original interval ( ) 2 b a x + = and testing to see in which of the subintervals ( ) [ ] 2 , b a a + or ( ) [ ] b b a , 2 + the solution lies. The procedure then repeated with the new interval as often as needed to locate the solution to desired accuracy. It is interesting to note that convergence is guaranteed. The advantages and limitations of this method are well documented in the existing mathematical literature (see, e.g., [3] and references therein). Here are some explanations how the macro %Bisect(eqtn,x1,x2,eps) works. Eqnt variable represents the equation under consideration, e.g., x x 3 4 3 = + or 0 4 3 3 = + − x x . 1 x is the starting point, and 2 x is the ending point of the interval, where the root will be searched, i.e., 2 1 x root x ≤ ≤ . These values are chosen in such a way that the condition 0 ) ( ) ( 2 1 < x f x f is satisfied. The method terminates the interval bracketing the root has the length eps root f ≤ ) ( , where eps defines the accuracy of the calculation. The comments are presented in the code below. THE CODE FOR FINDING THE ROOTS USINF BISECTION METHOD: %macro Bisect(eqtn,x1,x2,eps); %global root; * standardizing the initial equation, i.e., f1(x)=f2(x) is presented as f(x)=0 ; %let eqtn2=%scan(&eqtn,1,"="); %let eqtn3=%scan(&eqtn,2,"="); %let eqt=%str(y=&eqtn2-&eqtn3;); data one(keep=root); flag=0; x=&x1; &eqt; y1=y; calculating ( ) 1 x f and saving it as 1 y ; x=&x2; &eqt; y2=y; calculating ( ) 2 x f and saving it as 2 y ; x1=&x1; x2=&x2; *if the root is not between the initial values. In this case, one need to change their values. if y1*y2 > 0 then do; put 'No root within these borders!'; end; *the initial range is divided in 2 sub ranges and one of them, f(x1)*f(x2)<0 is chosen ; else do until ( flag=1); x3=0.5*(x1+x2); x=x3; &eqt; y3=y; calculating ( ) 3 x f and saving as 3 y ; if y1*y3<0 then x2=x3; * sets the middle point x3 as the end point of the new interval , while the starting point (x1) remains unchanged ; else x1=x3; * sets the middle point x3 as the staring point of the new interval, while the ending point (x2) remains unchanged ; *the process repeats until the the required accuracy is achieved. If the middle point is less than the defined accuracy then the corresponding value of x3 is considered as a root; if abs(y3)< &eps then do; flag=1; root=x3; call symput('root',root); end; * the root is presented in the form of a global macro variable ; end; run;