Statistical programming, matrix languages, and more

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

Accepted Solution Solved
Reply
Occasional Contributor
Posts: 14
Accepted Solution

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

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 !!!


Accepted Solutions
Solution
‎06-11-2015 01:51 PM
SAS Super FREQ
Posts: 3,413

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

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.

View solution in original post


All Replies
Solution
‎06-11-2015 01:51 PM
SAS Super FREQ
Posts: 3,413

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

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.

SAS Super FREQ
Posts: 3,413

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

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

Occasional Contributor
Posts: 14

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

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

Respected Advisor
Posts: 4,608

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

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

PG
SAS Super FREQ
Posts: 3,413

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

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

Grand Advisor
Posts: 9,584

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

Rick.

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

I get -48.16637008  not approximate 0 .

SAS Super FREQ
Posts: 3,413

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

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.

Occasional Contributor
Posts: 14

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

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

Grand Advisor
Posts: 9,584

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

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

Respected Advisor
Posts: 4,608

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

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

PG
Grand Advisor
Posts: 9,584

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

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

Maybe it is the problem about accuracy .

SAS Super FREQ
Posts: 3,413

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

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;

Grand Advisor
Posts: 9,584

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

Thank you . PG .

Occasional Contributor
Posts: 14

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

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;

☑ This topic is SOLVED.

Need further help from the community? Please ask a new question.

Discussion stats
  • 25 replies
  • 868 views
  • 7 likes
  • 5 in conversation