roots(x) vs root= interp1(y,x,0)

I tried to use the root = interp1 command to find the root of my function and it does not find the correct root, however, the roots command does. what is the difference?

x= linspace(-100,100,100);

y=85.*x- 91.*(x.^2)+ 44.*(x.^3)- 8.*(x.^4)+(x.^5)-26;

plot(x,y)

grid on

p=[1 -8 44 -91 85 -26]

r=roots(p)

thanks

Answer 1

[the cyclist]

It's not clear, at least to me, what you mean by "root = interp1 command". But, my best guess is that if you were somehow trying to use interp1 to find roots, you would run into the problem that it does linear interpolation by default, and root-finding for your equation will be non-linear.

Can you post the code that you thought should work, but did not give the result you expected?

[Matt Thomas]

root=interp1(y,x,0) rather than roots(p)

the interp1 command yielded 0.9xxx rather than 0.5xxx

Answer 2

[Paul]

The query points in x are spaced too far apart for the linear interpolation used by interp1 to be accurate relative to the actual curve. Plotting the curve around the real root shows why the original data yields an apparent root at r = 0.9xxxx when the actual root is around r = 0.5xxx

p = [1 -8 44 -91 85 -26];
x1 = linspace(-100,100,100);
y1 = polyval(p,x1);
x2 = linspace(-100,100,1000);
y2 = polyval(p,x2);
plot(x1,y1,'-x',x2,y2,'-o'),grid
xlim([-1.1 2.1])

roots(p)

ans = 
   2.5907 + 4.4333i
   2.5907 - 4.4333i
   1.1308 + 0.7012i
   1.1308 - 0.7012i
   0.5570 + 0.0000i

[interp1(y1,x1,0) interp1(y2,x2,0)]

ans = 1×2

0.9713 0.5702

[Matt Thomas]

so the linspace needs to have a tighter tolerance then?

thanks

[Torsten]

You need more points with a smaller distance between the x-values.

The maximum distance between the x-values is an estimate of the maximum error in determining the root of the function.

x= linspace(-3,3,1000);
y=85.*x- 91.*(x.^2)+ 44.*(x.^3)- 8.*(x.^4)+(x.^5)-26;
x0 = interp1(y,x,0)

Answer 3

[John D'Errico]

You need to understand the mathematics. Lost you will be without that (to misquote someone named Yoda.)

For this polynomial:

syms x
y=85.*x- 91.*(x.^2)+ 44.*(x.^3)- 8.*(x.^4)+(x.^5)-26;
fplot(y,[-2,5])
grid on

There seems to be a root near zero. Thus a root is where the function has a value of y=0. But when it is flat like that, this is often an indication there may be multiple roots, or at least several roots that are very close to each other. Or you may have a complex root relatively close to that point too. And the scale on the y axis makes it difficult to be sure.

grid on

p=[1 -8 44 -91 85 -26];
roots(p)

ans = 
   2.5907 + 4.4333i
   2.5907 - 4.4333i
   1.1308 + 0.7012i
   1.1308 - 0.7012i
   0.5570 + 0.0000i

However, roots tells us there is only one real root of the polynomial. It lies at x = 0.5570. while the other roots are complex.

What happens when you tried to use interp1? What does interp1 do here? Actually, interp1 is a terribly bad idea. When you flip x and y like that, for a function that has a nearly zero slope, then flipping y and x produces a function with nearly infinite slope. And using an INTERPOLATION function to approximate something with a nearly infinite slope is a BAD IDEA.

You used

x = linspace(-100,100,100);

so 100 points, that varied from -100 to 200. That means a spacing of roughly 2 between each point. Do you see how much curvature you have in that function? Then using interp1, which by default uses LINEAR interpolation is just completely misunderstanding what is happening. The only points that are near zero in x at all were:

x(49:53)

ans = 1×5

-3.0303 -1.0101 1.0101 3.0303 5.0505

yfun = matlabFunction(y)

yfun = function_handle with value:
    @(x)x.*8.5e+1-x.^2.*9.1e+1+x.^3.*4.4e+1-x.^4.*8.0+x.^5-2.6e+1

plot(yfun(x(49:53)),x(49:53),'o-')
grid on

the problem is, if you use interp1 here, it uses only linear interpolation. And that means it tries to interpolate the connect-the-dots curve drawn there. Using interp1 here to solve for a root of this function is again, just a flat out bad idea. I know that people have told you to do that online. But they were wrong, at least to tell you to just swap x and y, and then to use interp1, if done without any understanding of what is happening.) That is certainly true if you use such a coarsely sampled grid in x.

CAN we do a better job, using interp1 more intelligently? Well, yes. By using a more finely sampled grid, we might do this:

xbetter = linspace(-10,10,100);
ybetter = yfun(xbetter);
xinterproot = interp1(ybetter,xbetter,0)

xinterproot = 0.5698

xbetter = linspace(-10,10,1000);
ybetter = yfun(xbetter);
xinterproot = interp1(ybetter,xbetter,0)

xinterproot = 0.5572

xbetter = linspace(-10,10,10000);
ybetter = yfun(xbetter);
xinterproot = interp1(ybetter,xbetter,0)

xinterproot = 0.5570

So you see that by use of a finely sampled grid, the LINEAR interpolation can do a better job. However, we could also have used a spline interpolant with interp1. This is still a dangerous thing to do with interp1 using a spline, because splines themselves have serious problems.

xbetter = linspace(-10,10,100);
ybetter = yfun(xbetter);
xinterproot = interp1(ybetter,xbetter,0,'spline')

xinterproot = 0.5549

xbetter = linspace(-10,10,1000);
ybetter = yfun(xbetter);
xinterproot = interp1(ybetter,xbetter,0,'spline')

xinterproot = 0.5570

which you can see that even on the moderately coarsely sample grid, did a decent job. I did not need to go to as finely sampled a grid to get 4 digits correct there, when using a spline.

Just use roots, which will produce the correct solution, here accurate to full double precision accuracy. Or, you could have used solve or vpasolve. Since this is a 5th degree polynomial, we cannot compute an algebraic form of the roots of a 5th degree polynomial anyway, so vpasolve is correct here.

vpasolve(y)

ans = 

[Matt Thomas]

thanks for that nice reply with all the examples. Yeah, i am very familar with the math. Im in my 13 semester of math at the moment. However, i have never used matlab or excell or any type of programing. Im old school and do iterations by hand, usually 4 of them.

Usually what i do is take the deriv of the function and use that to find my slopes of 0, concave and convex points, then just graph it. But had to use matlab for this one haha.

i have sent this thread to my class so they can see what is going on. most of them are really good at matlab, except me haha