root of a polynomial fails to satisfy the polynomial

Please see the last loop. L must be zero. Because DE is polynomial whose roots are given by h. When we substitute the roots in the DE, we must get zero. But its not the case here.

I am not able to detect the mistake. please help

function sukh1
eps=0.;
beta_T=2e-4;
lambda=77.6;mu=38.6;C_e=383;k=400; rho=8.920;
T_0=318;tau_0=1e-12; tau_1=2*tau_0;
bulk=lambda+2*mu;
for m=10:10:10;
    omg=2*pi*(2^m);
alpha_T=(3*lambda+2*mu)*beta_T;
A=bulk*k;
B=rho*C_e*(bulk)*(tau_0+1i/omg)+rho*k+(alpha_T^2)*T_0*(1-1i*omg(1-eps)*tau_1)*(eps*tau_0+1i/omg);
C=(rho*rho*C_e*(tau_0+i/omg));
%Instead of alpha, alpha^2 is used everywhere. so we will directly find
%alpha^2 and denote it by alf2.
alf2=roots([C -B A]);
%dialational wave with velocity alpha1(alpha2)is called qP(q(T))
nu1=(alpha_T*T_0*omg*omg*(eps*tau_0+i/omg))/(rho*C_e*alf2(1)*(tau_0+i/omg)-k);
nu2=(alpha_T*T_0*omg*omg*(eps*tau_0+i/omg))/(rho*C_e*alf2(2)*(tau_0+i/omg)-k);
eta=nu1/nu2;
%gm stands for gamma;
gm1=bulk/rho*alf2(1);
gm2=bulk/rho*alf2(2);
a=(gm1-eta*gm2); b=1-eta;
%beta=sqrt(mu/rho), bets=beta^2
bets=(mu/rho);
e1=bets/alf2(1);e2=bets/alf2(2);
d=a*b-(e1+eta*eta*e2); 
g=a*a+b*b+4*d; es=e1+e2; ep=e1*e2;
beta=sqrt(bets);
C0=1024*eta*(d+eta*es);
C1=256*(d*d-eta*(g+4*d))-1024*eta*eta*(ep+2*es);
C2=128*(2*eta*a*(a+b)-(d-2*eta)*g)+1024*(eta^2)*(2*ep+es);
C3=16*((g^2)+8*a*(a+b)*(d-2*eta)-4*eta*(a^2))-1024*(eta^2)*ep;
C4=-32*a*(a*(d-2*eta)+(a+b)*g);
C5=8*(a^2)*(3*g+4*a*b-8*d);
C6=-8*(a^3)*(a+b);
C7=a^4;  
h=roots([C6 C7]);
for j=1:1;
    hj=h(j)
    DE=C6*(hj)+C7;
    L=abs(DE)
    end
end

Answer 1

[Walter Roberson]

Your C6 and C7 are in the range of 10^39. eps() of those is around 10^24, so calculations involving C6 and C7 cannot be expected to be more accurate than +/- 10^24. Your L is on the order of 10^23, so you are getting as much accuracy out of the calculation as is possible.

If you change your calculations over to symbolic, then L comes out 0, but C6 and C7 come out quite different. That might be because the order of roots() for symbolic expressions is not certain to be the same as the order of roots() for numeric expressions. In turn, the order of complex roots is not defined -- so your nu1 and nu2 should not be counting on a particular ordering like they are now.

[Sukhversha Luthra]

so what should i do for getting my polynomial satisfied by the roots

[Walter Roberson]

You know that C6*x + C7 has a single root:

C6*x + C7 = 0
C6*x = 0 - C7
x = (0 - C7) / C6
x = - C7 / C6

This is going to give you as an accurate a solution as is possible considering the finite floating point values.

You can convert your code to use the Symbolic Toolbox so that you are working with algebraic roots in closed form (in which case you get an exact 0 for L).

You can use the Multiprecision Toolbox from the File Exchange, https://www.mathworks.com/matlabcentral/fileexchange/6446-multiple-precision-toolbox-for-matlab . However since the results are going to be finite length, no matter how many specific finite digits you tell it to use, you are not going to get an exact 0 out.

The situation is similar to saying that you have to express 1/3 in decimal to some finite number of digits, and then you want to multiply the result by 3 and get an exact 1 out.

0.33333333333333333333
* 3 =
0.99999999999999999999

no matter how many more digits you add to the list of 3s, after the multiplication you are not going to recover an exact 1. Here I used 20 decimal digits, which would be between 66 and 67 binary digits of precision, but it is not enough -- no finite number of digits is enough in this situation.

Because no finite number of digits is enough, if you do choose to work with a finite number of digits, you need to compare to 0 within a tolerance. 20 decimal digits, for example, would not be able to distinguish between numbers that were closer together than 1 part in 10^20, so two calculations would have to be considered to give the same result if the results were no more than 1 part in 10^20 apart for a 20 decimal digit scheme.

The corresponding fraction for binary floating point is 1 part in 2^52. MATLAB gives a special name to this quantity: eps. 1 part in 2^52 times the 10^39 range you are working with is somewhere around 10^23, so instead of looking for an exact 0 you should be checking for something on the order of +/- 10^24: such a result cannot be told apart from 0 considering your range of values.