find roots of function

[WeiChen]

I confused why i get 3.3 for root of function for froot call. I know function have root at x=1 because 1-4+2+1 = 0. So why froot call not give me root>?

proc iml;
start func(x);
return x**3 - 4*x**2 + 2*x + 1;
finish;

root = froot("func", {-4 4});
print root;
y = func(root);
print y;

[Rick_SAS · Accepted solution]

The FROOT function returns one root within the interval that you specify. The value x=3.3 is a root, just not the one that you want!

So how can you find the other roots? Read the article, "Finding roots: Automating the search for an initial guess." It suggests graphing the function on a domain that is wide enough to display all the roots. FOr your function, I converted the formula to use the '#' multiplication operator so that I can pass in a vector of x values with one call:

proc iml;
/* vectorize computation so that x can be a vector */
start func(x);
   return x##3 - 4*x##2 + 2#x + 1;
finish;

/* Graph function on wide domain. See
   https://blogs.sas.com/content/iml/2015/06/22/root-guess.html
*/
x = do(-2, 4, 0.1);
y = func(x);
call scatter(x, y) grid={x y};

 

 From the graph, it looks like you can subdivide the x axis into nonoverlapping intervals so that there is one root in each interval. The FROOT function will then compute the three roots of this function:

intervals = {-4 0,
              0 2,
              2 4};
roots = froot("func", intervals);    
y = func(roots);
print roots y;

roots	y
-0.302776	0
1	0
3.3027756	1.776E-15

[WeiChen]

What if i dont know how many roots are in function?

[Rick_SAS]

Read the article. It discusses that situation.

[Ksharp]

Rick,

It is about to find the roots of the polynomial .Shouldn't use POLYROOT()?

proc iml;
p={1 -4 2 1};
r=polyroot(p);
print r;
quit;

[Rick_SAS]

Yes, POLYROOT would work for this example. It wasn't clear to me whether the OP was actually interested in a polynomial or was interested in a harder problem and is using the polynomial as an example.