Could you please suggest a program to find the zeros to these two polynomials

[Thomas_mp]

Could you please suggest a program to find the zeros to these two polynomials:

48 = ((0.3-x)*2)/((1+x)^2)  solution (using excel) is -0.78717

25 = ((0.1-x)*4)/(1.x)^2  sol is -0.65292

I have a list of 4,000 of these equations to solve. and the actual polynomial have 12 terms.  Perhaps if I know how to do these simple ones, I can do the more complex ones.

Thank you !!!

[Rick_SAS · Accepted solution]

First subtract the constant terms so that you have a polynomial of the form P(x)=0. Then find the roots of the polynomials.

Do you want only the real roots, or do you want complex roots as well?

The POLYROOT function finds all complex roots of a polynomial, but you would have to represent the functions as a12*x^12 + a11*x^11 + ... + a0. See the doc or the first example of "Finding the roots of a univariate function".

To find only the real roots, use the FROOT function. Or find all and throw out the complex roots. An example of the FROOT function is given in "A simple way to find the root of a function of one variable." For the FROOT function, you can keep the functions written as the are.

[Rick_SAS]

Um, I just reread your question. Those are not polynomials, they are rational functions.

[Thomas_mp]

Thank very much Rick.

How could use this to solve several thousands functions with same "functional form"?  Please, let me explain:

In your example  in the link http://blogs.sas.com/content/iml/2014/02/05/find-the-root-of-a-function.html

return( exp(-x##2) - x##3 + 5#x +1 );

let me call -2,-3, 5,1 "parameters"

I have a file with 4,000 different rows with these parameters; thus each row defines a rational function; each has (at least) one root in the interval (0.05,0.30)

I tried to read the file with the 4,000 rows (I have it in a csv format) ,

infile  'H;file.csv' ; input a b c d;

The parameters in your example are  (a,b,c,d), so your example is, in a generic form: ;

( exp(-x##a) - bx##c + c#x +d );

I run this, but I get the message

  z = froot("Func", intervals);

ERROR: (execution) Matrix has not been set to a value.

Why?

As you can tell my knowledge of SAS is very limited..

Thank you for your help

[PGStats]

Here is the method to find the function roots with FCMP. Note however that the SOLVE function in FCMP requires a fairly good guess as initial value to find the root.

/*  Equation has the form p1 = (p2-x)*p3/(p4+x)**p5   */

data parms;

eqNo + 1;

input p1-p5;

datalines;

48 0.3 2 1 2

25 0.1 4 1 2

;

/* Define FCMP routine to find function roots */

proc fcmp outlib=sasuser.fcmp.test;

function fn(p2, p3, p4, p5, x);

    if fuzz(p4+x) = 0 then return (.);

    else return ((p2-x)*p3/(p4+x)**p5);

endsub;

function findZero(p1, p2, p3, p4, p5, init);

    array solvopts[1] initial (0);

    initial = init;

    return (solve("fn", solvopts, p1, p2, p3, p4, p5, .));

endsub;

run;

options cmplib=sasuser.fcmp;

data roots;

set parms;

root = findZero(of p1-p5, -0.6); /* -0.6 is the initial guess  */

run;

proc print data=roots noobs; run;

PG

[Rick_SAS]

I think this is a difficult problem in general because root-finding works best when you have knowledge of the function and the interval on which you are searching.  As PG says, a good estimate for the root is important.   This problem is compounded because these rational functions are not continuous and have singularities. However, in your recent response you claim that function "has (at least) one root in the interval (0.05,0.30)."

Is it always true that the function is continuous on that interval and f(0.05) and f(0.3) have different signs? If so, then you can use the FROOT function in SAS/IML and pass in that interval.  The function will return the first root that it finds in that interval.

The following example reads the parameters from a SAS data set and calls the FROOT function for each row.

data parms;
eqNo + 1;
input p1-p5;
datalines;
48 0.3 2 1 2
25 0.1 4 1 2
;

proc iml;
/* Define objective function */
start fn(x)  global( p );
  return( (p[2]-x)*p[3]/(p[4]+x)##p[5] - p[1] );
finish;

varNames="p1":"p5";
use parms;
read all var varNames into m;
close parms;

root = j(nrow(m), 1, .);
do i = 1 to nrow(m);
   p = m[i,];           /* store i_th row into global variable */
   /* Optional: Call function that computes interval that has root */
   interval = {-0.9 2};
   root = froot("fn", interval);  /* find a root in the interval [0,1] */
end;

print root;

Message was edited by: Rick Wicklin

[Ksharp]

Rick.

Not sure . If I am wrong. Your root is not right for the first function .

I get -48.16637008  not approximate 0 .

[Rick_SAS]

Thanks! I was reading all numeric variables, so I was getting the EqNo variable and using that for the first parameter.  I have corrected the program to read only the parameters.

[Thomas_mp]

Thank you again !!

This was my first experience in this forum and I am very grateful.

The 3 methods work for me.  I will be using Rick's because I understand a little bit better how it works.  For your information, I will be using this to solve for the implied cost of capital as explained in the paper :

Toward an Implied Cost of Capital

Author(s): William R. Gebhardt, Charles M. C. Lee and Bhaskaran Swaminathan

Source: Journal of Accounting Research, Vol. 39, No. 1 (Jun., 2001), pp. 135-176

Take care

[Ksharp]

By the way, I found a better ROOT than you for the first funtion via GA. Of course for the second function .

1 OPTIONS NONOTES NOSTIMER NOSOURCE NOSYNTAXCHECK;

57 

58 

59 data _null_;

60 my_root=-1.254501;

61 your_root=-0.787165;

62 my_value=(0.3-my_root)*2/(1+my_root)**2 - 48;

63 your_value=(0.3-your_root)*2/(1+your_root)**2 - 48;

64 put my_root= my_value= / your_root= your_value=;

65 run;

my_root=-1.254501 my_value=0.000085964

your_root=-0.787165 your_value=-0.00020722

NOTE: DATA statement used (Total process time):

  real time 0.00 seconds

  cpu time 0.00 seconds

Xia Keshan

[PGStats]

it is not a better root. The function has two roots! :smileysilly:  - PG

[Ksharp]

Really ?  Why ?  But -1.254501 is more approximate 0 . See my LOG .

Maybe it is the problem about accuracy .

[Rick_SAS]

As I said, there is a singularity, and the function has a root to the left of the singularity and another branch to the right.

Here is a graph of the function, which has a singularity at x = -p4.:

data rational;
p1=48; p2=0.3; p3=2; p4=1; p5=2;
/* p1=25; p2=0.1; p3=4; p4=1; p5=2; */
delta = 0.05;
do x = -3 to -p4-delta by delta;
   Branch = "Neg";
   y =  (p2-x)*p3/(p4+x)**p5 - p1;
   output;
end;
do x = -p4+delta to 1 by delta;
   Branch = "Pos";
   y =  (p2-x)*p3/(p4+x)**p5 - p1;
   output;
end;
run;

proc sgplot data=rational;
series x=x y=y / group=Branch;
refline 0 /axis=y;
refline -1 /axis=x;
yaxis max =100;
run;

[Ksharp]

Thank you . PG .

[Thomas_mp]

Hello again Rick,

Thank you gain,

I thought that your program worked for my generic problem , but I am having some trouble .  In my problem one of the terms is of the form:

x/(x*(1+x)##12) 

For these type of problems, when x appears also multiplying in the denominator,  it seems that the program cannot find the solution.  For instance, in your example, below, slightly modified (I changed the first number from 48 to 3.3058, and added the x multiplying in the denominator), the solution, a zero, is   x= 0.1, but the program seems to have problems finding it.

:

data parms;

eqNo + 1;

input p1-p5;

datalines;

3.3058 0.3 2 1 2

;

proc iml;
/* Define objective function */
start fn(x)  global( p );
  return( (p[2]-x)*p[3]/x*(p[4]+x)##p[5] - p[1] );
finish;

varNames="p1":"p5";
use parms;
read all var varNames into m;
close parms;

root = j(nrow(m), 1, .);
do i = 1 to nrow(m);
   p = m[i,];           /* store i_th row into global variable */
   /* Optional: Call function that computes interval that has root */
   interval = {-0.9 2};
   root = froot("fn", interval);  /* find a root in the interval [0,1] */
end;

print root;