Это ошибка в NSolve в Mathematica?

Question

Можно было бы ожидать и надеяться, что если вы попросите Mathematica найти корни многочлена, он должен дать одинаковые (приблизительные) ответы, делаете ли вы это символически, а затем найти числовые приближения к этим точным ответам или сделать это численно. Вот пример, который (в Mathematica 7, работает на OS X), где это терпит неудачу:

poly = -112 + 1/q^28 + 1/q^26 - 1/q^24 - 6/q^22 - 14/q^20 - 25/q^18 -
  38/q^16 - 52/q^14 - 67/q^12 - 81/q^10 - 93/q^8 - 102/q^6 - 108/
  q^4 - 111/q^2 - 111 q^2 - 108 q^4 - 102 q^6 - 93 q^8 - 81 q^10 -
  67 q^12 - 52 q^14 - 38 q^16 - 25 q^18 - 14 q^20 - 6 q^22 - q^24 +
  q^26 + q^28;

Total[q^4 /. NSolve[poly == 0, q]] - Total[q^4 /. N[Solve[poly == 0, q]]]

(Примечание: это на самом деле многочлен Лорана, и если вы умножите его на большую степень q, проблема исчезнет.)

Последняя строка здесь - просто демонстрация того, что найденные решения очень разные; фактически это количество, которое мы пытались вычислить в задаче, над которой работали.

Если вы внимательно посмотрите на вывод NSolve[poly == 0, q] и N[Solve[poly == 0, q], вы увидите, что NSolve дает только корни 54 вместо ожидаемых 56. Дело не в том, что он просто пропустил повторный корень или что-то еще; в нем отсутствуют два самых больших корня по величине (приблизительно +/- 1.59)

Это ошибка в Mathematica? У кого-нибудь есть объяснение, почему это происходит?

Answer 1

К сожалению, ваши ожидания необоснованны.

Как вы говорите, Solve[] дает точное решение, а N[] вносит небольшую ошибку, но только один раз, когда вы ее оцениваете. NSolve[] с другой стороны использует числовые приближения с самого начала, и, таким образом, может накапливаться значительная ошибка округления.

Вы также ограничены точностью вычисления по умолчанию, и это может привести к полному отказу численного метода, например, отсутствию корней (см. Полином Уилкинсона ). Вы можете противодействовать этому, сказав NSolve[], чтобы использовать большую точность, например так:

In[1]  := Total[q^4 /. NSolve[poly == 0, q, WorkingPrecision -> 50]] - 
          Total[q^4 /. N[Solve[poly == 0, q]]]
Out[1] := 0. - 3.66374*10^-15 I

При использовании численных методов всегда важно помнить об ошибках. Поскольку это справедливо для большого разнообразия численных анализов задач от решения длинных полиномов до диагонализации больших матриц, интеграции странных функций и т. Д., Единого правильного подхода не существует, и Mathematica нужно сказать Например, повысить Рабочую точность или применить другую числовую технику.