BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
Thomas_mp
Obsidian | Level 7

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

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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

25 REPLIES 25
Rick_SAS
SAS Super FREQ

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
SAS Super FREQ

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

Thomas_mp
Obsidian | Level 7

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
Opal | Level 21

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
Rick_SAS
SAS Super FREQ

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
Super User

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
SAS Super FREQ

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
Obsidian | Level 7

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
Super User

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
Opal | Level 21

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

PG
Ksharp
Super User

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

Maybe it is the problem about accuracy .

Rick_SAS
SAS Super FREQ

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;

Thomas_mp
Obsidian | Level 7

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;

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 25 replies
  • 2107 views
  • 7 likes
  • 5 in conversation