Найти реальные корни многочлена

24

Напишите автономную программу, которая при задании полинома и границы найдет все действительные корни этого полинома для абсолютной ошибки, не превышающей границы.

Ограничения

Я знаю, что Mathematica и, возможно, некоторые другие языки имеют односимвольное решение, и это скучно, поэтому вам следует придерживаться примитивных операций (сложение, вычитание, умножение, деление).

Существует определенная гибкость в форматах ввода и вывода. Вы можете принимать ввод через аргументы stdin или командной строки в любом приемлемом формате. Вы можете разрешить использование чисел с плавающей запятой или требовать использования некоторого представления рациональных чисел. Вы можете взять границу или обратную величину границы, и если вы используете плавающую точку, вы можете предположить, что граница будет не менее 2 ulp. Полином должен быть выражен в виде списка мономиальных коэффициентов, но он может быть старшим или младшим.

Вы должны быть в состоянии объяснить, почему ваша программа всегда будет работать (по модулю числовые проблемы), хотя нет необходимости предоставлять полные доказательства в строке.

Программа должна обрабатывать многочлены с повторяющимися корнями.

пример

x^2 - 2 = 0 (error bound 0.01)

Ввод может быть, например,

-2 0 1 0.01
100 1 0 -2
1/100 ; x^2-2

Вывод может быть, например,

-1.41 1.42

но нет

-1.40 1.40

поскольку это имеет абсолютные ошибки около 0,014 ...

Контрольные примеры

Просто:

x^2 - 2 = 0 (error bound 0.01)

x^4 + 0.81 x^2 - 0.47 x + 0.06 (error bound 10^-6)

Многократный корень:

x^4 - 8 x^3 + 18 x^2 - 27 (error bound 10^-6)

Полином Уилкинсона:

x^20 - 210 x^19 + 20615 x^18 - 1256850 x^17 + 53327946 x^16 -1672280820 x^15 +
    40171771630 x^14 - 756111184500 x^13 + 11310276995381 x^12 - 135585182899530 x^11 +
    1307535010540395 x^10 - 10142299865511450 x^9 + 63030812099294896 x^8 -
    311333643161390640 x^7 + 1206647803780373360 x^6 -3599979517947607200 x^5 +
    8037811822645051776 x^4 - 12870931245150988800 x^3 + 13803759753640704000 x^2 -
    8752948036761600000 x + 2432902008176640000  (error bound 2^-32)

NB. Этот вопрос был в Песочнице в течение приблизительно 3 месяцев. Если вы думаете, что это нужно улучшить перед публикацией, посетите Песочницу и прокомментируйте другие предложенные вопросы, прежде чем они будут опубликованы на Главном.

Питер Тейлор
источник
«Вы должны быть в состоянии объяснить, почему ваша программа всегда будет работать» На всякий случай, если кому-то понадобится помощь в этом,
доктор Белизариус
@belisarius, ??
Питер Тейлор
3
задумывался как шутка :(
Доктор Белизарий
Я знаю, что это старая проблема, поэтому не стесняйтесь отвечать, если вы не хотите снова ее открывать. (а) Можем ли мы написать функцию или только полную программу? (б) Если мы можем написать функцию, можем ли мы предположить, что для ввода используется какой-то удобный тип данных, например, Python fractions.Fraction(рациональный тип)? (c) Должны ли мы обращаться с полиномами степени <1? (d) Можем ли мы предположить, что ведущий коэффициент равен 1?
Ell
(e) Что касается многочленов с повторяющимися корнями, то стоит проводить различие между корнями нечетных и четных кратностей (в тестовых случаях есть только корни нечетных кратностей.) Хотя с корнями нечетной кратности не так уж сложно разобраться, я Я не уверен, насколько реально правильно численно обрабатывать корни четной кратности, тем более что вы указываете погрешность только для значений корней, а не для их существования. (...)
Ell

Ответы:

8

Математика, 223

r[p_,t_]:=Module[{l},l=Exponent[p[x],x];Re@Select[NestWhile[Table[#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}],{i,l}]&,Table[(0.9+0.1*I)^i,{i,l}],2*Max[Table[Abs[#1[[i]]-#2[[i]]],{i,l}]]>t&,2],Abs[Im[#]]<t^.5&]]

Это решение реализует метод Дюранда – Кернера для решения полиномов. Обратите внимание, что это не полное решение (как будет показано ниже), потому что я пока не могу обработать полином Уилкинсона с заданной точностью. Сначала объяснение того, что я делаю: код в формате Mathematica

#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}]&: Таким образом, функция вычисляет для каждого индекса iследующее приближение Дюранда-Кернера. Затем эта строка инкапсулируется в таблицу и применяется с помощью NestWhile к входным точкам, созданным с помощью Table[(0.9+0.1*I)^i,{i,l}]. Условием для NestWhile является то, что максимальное изменение (по всем терминам) от одной итерации к следующей больше, чем указанная точность. Когда все термины изменились меньше, чем NestWhile заканчивается и Re@Selectудаляет нули, которые не попадают на реальную линию.

Пример вывода:

> p[x_] := x^2 - 2
> r[p, .01]
{1.41421, -1.41421}

> p[x_] := x^4 - 8 x^3 + 18 x^2 - 27
> r[p, 10^-6]
{2.99999, 3., 3.00001, -1.}

> p[x_] := x^20 - 210 x^19 + ... + 2432902008176640000 (Wilkinson's)
> Sort[r[p, 2^-32]]
{1., 2., 3., 4., 5., 6., 7.00001, 7.99994, 9.00018, 10.0002, 11.0007, \
11.9809, 13.0043, 14.0227, 14.9878, 16.0158, 16.9959, 17.9992, \
19.0001, 20.}

Как вы, вероятно, можете видеть, когда степень возрастает, этот метод начинает отражаться вокруг правильных значений, никогда не возвращаясь полностью. Если я установлю условие остановки моего кода на что-то более строгое, чем «от одной итерации к следующей, догадки меняются не более, чем на эпсилон», тогда алгоритм никогда не останавливается. Я полагаю, я должен просто использовать Дюран-Кернер в качестве входных данных для метода Ньютона?

Kaya
источник
У Дюранда-Кернера также есть потенциальные проблемы с несколькими корнями. (Метод Ньютона тоже может не сильно помочь - полином Уилкинсона определенно выбран плохо обусловленным).
Питер Тейлор,
Вы совершенно правы: я отказался от этого курса действий после приближения к Уилкинсону около x = 17, это абсолютный беспорядок. Я волнуюсь, что мне придется пойти на символическое решение на основе Гребнера, чтобы получить гораздо большую точность.
Кая