topic find roots of function in SAS/IML Software and Matrix Computations

find roots of function

WeiChen — Sat, 05 Feb 2022 11:06:49 GMT

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;

Re: find roots of function

Rick_SAS — Sat, 05 Feb 2022 11:16:21 GMT

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

Re: find roots of function

WeiChen — Sat, 05 Feb 2022 11:18:31 GMT

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

Re: find roots of function

Rick_SAS — Sat, 05 Feb 2022 11:19:43 GMT

Read the article. It discusses that situation.

Re: find roots of function

Ksharp — Sat, 05 Feb 2022 14:13:06 GMT

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;

Re: find roots of function

Rick_SAS — Sat, 05 Feb 2022 20:59:53 GMT

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.