Skip to main content
Fixed code
Source Link
Michael E2
  • 244.7k
  • 18
  • 351
  • 774

The roots returned by Solve are nearly as accurate as possible within rounding error at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be significantly more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, dx$, where $dx$ is the smallest possible change in $x$. The smallest possible change in $x$ is about $x\cdot\epsilon$$|a_k x^k| \, \epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $4\cdot10^{16}$$1.6\cdot10^{11}$.

Here are the rounding-error bounds of the largest monomial terms and the residual errors in "y" for each root, which are all about the same size:

Abs[List @@ eqn[l]*l*$MachineEpsilon]eqn[l]*$MachineEpsilon] /. sol // Map@Max
(*  {31.8073*10^166153*10^11, 43.7861*10^4102*10^-98, 51.9168*10^1990*10^-1410, 43.7861*10^4102*10^-98, 31.8073*10^166153*10^11}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a significantly different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will need at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 indicates that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.

The roots returned by Solve are as accurate as possible at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, dx$, where $dx$ is the smallest possible change in $x$. The smallest possible change in $x$ is about $x\cdot\epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $4\cdot10^{16}$.

Here are the rounding-error bounds of the largest monomial terms and the errors in "y" for each root:

Abs[List @@ eqn[l]*l*$MachineEpsilon] /. sol // Map@Max
(*  {3.8073*10^16, 4.7861*10^-9, 5.9168*10^-14, 4.7861*10^-9, 3.8073*10^16}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will need at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 indicates that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.

The roots returned by Solve are nearly as accurate as possible within rounding error at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be significantly more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, \epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $1.6\cdot10^{11}$.

Here are the rounding-error bounds of the largest monomial terms and the residual errors in "y" for each root, which are all about the same size:

Abs[List @@ eqn[l]*$MachineEpsilon] /. sol // Map@Max
(*  {1.6153*10^11, 3.4102*10^-8, 1.1990*10^-10, 3.4102*10^-8, 1.6153*10^11}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a significantly different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will need at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 indicates that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.

Added missing braces to TeX expression, emphasized last paragraph
Source Link
MarcoB
  • 67.7k
  • 18
  • 96
  • 198

The roots returned by Solve are as accurate as possible at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, dx$, where $dx$ is the smallest possible change in $x$. The smallest possible change in $x$ is about $x\cdot\epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $4\cdot10^16$$4\cdot10^{16}$.

Here are the rounding-error bounds of the largest monomial terms and the errors in "y" for each root:

Abs[List @@ eqn[l]*l*$MachineEpsilon] /. sol // Map@Max
(*  {3.8073*10^16, 4.7861*10^-9, 5.9168*10^-14, 4.7861*10^-9, 3.8073*10^16}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will need at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 means it turns outindicates that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.

The roots returned by Solve are as accurate as possible at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, dx$, where $dx$ is the smallest possible change in $x$. The smallest possible change in $x$ is about $x\cdot\epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $4\cdot10^16$.

Here are the rounding-error bounds of the largest monomial terms and the errors in "y" for each root:

Abs[List @@ eqn[l]*l*$MachineEpsilon] /. sol // Map@Max
(*  {3.8073*10^16, 4.7861*10^-9, 5.9168*10^-14, 4.7861*10^-9, 3.8073*10^16}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 means it turns out that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.

The roots returned by Solve are as accurate as possible at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, dx$, where $dx$ is the smallest possible change in $x$. The smallest possible change in $x$ is about $x\cdot\epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $4\cdot10^{16}$.

Here are the rounding-error bounds of the largest monomial terms and the errors in "y" for each root:

Abs[List @@ eqn[l]*l*$MachineEpsilon] /. sol // Map@Max
(*  {3.8073*10^16, 4.7861*10^-9, 5.9168*10^-14, 4.7861*10^-9, 3.8073*10^16}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will need at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 indicates that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.

Source Link
Michael E2
  • 244.7k
  • 18
  • 351
  • 774

The roots returned by Solve are as accurate as possible at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, dx$, where $dx$ is the smallest possible change in $x$. The smallest possible change in $x$ is about $x\cdot\epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $4\cdot10^16$.

Here are the rounding-error bounds of the largest monomial terms and the errors in "y" for each root:

Abs[List @@ eqn[l]*l*$MachineEpsilon] /. sol // Map@Max
(*  {3.8073*10^16, 4.7861*10^-9, 5.9168*10^-14, 4.7861*10^-9, 3.8073*10^16}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 means it turns out that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.