Have a trouble with the function roots

Question

Hey guys I have multiple problems with using function 'roots'.

I Have to find zeros of 's^1000 + 1'. I made Y = zeros(1,1000) then manually changed the 1000th matrice to '1'. but then 'root' function does not work with it !

---

Another problem is that I am having trouble with matrix multiplication. The question is finding zeros(roots) of (s^6 + 6*s^5 + 15*s^4 + 20*s^3 + 15*s^2 + 6*s +1)*(s^6 + 6s^5 + 15*s^4 +15*s^2 +6*s +1) so i did:

a = [1 6 15 20 15 6 1]
b = [1 6 15 0 15 6 1]
y = a.*b;
roots(y)

but this gives me

-27.9355 + 0.0000i
  -8.2158 + 0.0000i
   0.1544 + 0.9880i
   0.1544 - 0.9880i
  -0.1217 + 0.0000i
  -0.0358 + 0.0000i

where I calculate the original equation with wolfram then I have made matrix as :

p = [1 12 66 200 375 492 524 492 375 200 66 12 1]
roots(p)

and this gives me :

-3.1629 + 2.5046i   
-3.1629 - 2.5046i
 0.3572 + 0.9340i
 0.3572 - 0.9340i   
-1.0051 + 0.0000i   
-1.0025 + 0.0044i   
-1.0025 - 0.0044i   
-0.9975 + 0.0044i   
-0.9975 - 0.0044i   
-0.9949 + 0.0000i   
-0.1943 + 0.1539i  
 -0.1943 - 0.1539i

and I think the second solution is right (that is what wolfram alpha gave me) How would you answer these two questions through matlab guys?

Answer 1

To multiply polynomials, you convolve their coefficients:

>> roots(conv(a,b))

ans =

  -3.1629 + 2.5046i
  -3.1629 - 2.5046i
   0.3572 + 0.9340i
   0.3572 - 0.9340i
  -1.0051          
  -1.0025 + 0.0044i
  -1.0025 - 0.0044i
  -0.9974 + 0.0044i
  -0.9974 - 0.0044i
  -0.9950          
  -0.1943 + 0.1539i
  -0.1943 - 0.1539i

Answer 2

Q1

Using roots to find the roots of s¹⁰⁰⁰ + 1 is a bit of an overkill. The solution is given by this code snippet (corrected the first version; may be deduced using De Moivre's formula):

n = 1000;
k = 0:n-1
u = (2*k + 1)*pi / n;
s = cos(u) + 1i*sin(u)

Also, this method is approx. 100000 times faster.

Q2

Multiplying two polynomials to find the roots of their product is a bit of an overkill. :-) The union of the two polynomials' root sets is the root set of the product polynomial:

s = [roots(a);roots(b)]

Also, this method is more accurate.

Answer 3

1) The vector describing s^1000 + 1 should end with a 1 as well.

2)

a = [1 6 15 20 15 6 1]
b = [1 6 15 0 15 6 1]
y = a.*b;

This is a DOT product, multiplication of polynomials do not multiply element-wise.

Answer 4

Question 1

You need to include the coefficient of x^0 in the vector of coefficients, so there are 1001 entries with the first and last being 1

coeffs=zeros(1001,1);
coeffs([1,1001])=1;
roots(coeffs)

Question 2

To multiply the coefficients of polynomials you need to use convolution:

roots(conv(a,b))