Trouver les vraies racines d'un polynôme

Écrivez un programme autonome qui, lorsqu'il reçoit un polynôme et une borne, trouvera toutes les racines réelles de ce polynôme à une erreur absolue ne dépassant pas la borne.

Contraintes

Je sais que Mathematica et probablement d'autres langues ont une solution à un seul symbole, et c'est ennuyeux, vous devez donc vous en tenir aux opérations primitives (addition, soustraction, multiplication, division).

Il existe une certaine flexibilité sur les formats d'entrée et de sortie. Vous pouvez prendre des entrées via stdin ou des arguments de ligne de commande dans n'importe quel format raisonnable. Vous pouvez autoriser la virgule flottante ou exiger qu'une représentation des nombres rationnels soit utilisée. Vous pouvez prendre la limite ou l'inverse de la limite, et si vous utilisez une virgule flottante, vous pouvez supposer que la limite ne sera pas inférieure à 2 ulp. Le polynôme doit être exprimé sous la forme d'une liste de coefficients monétaires, mais il peut être de grande ou de petite taille.

Vous devez être en mesure de justifier pourquoi votre programme fonctionnera toujours (problèmes modulo numériques), bien qu'il ne soit pas nécessaire de fournir des preuves complètes en ligne.

Le programme doit gérer des polynômes avec des racines répétées.

Exemple

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

L'entrée pourrait être par exemple

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

La sortie pourrait être par exemple

-1.41 1.42

mais non

-1.40 1.40

car cela a des erreurs absolues d'environ 0,014 ...

Cas de test

Simple:

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

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

Racine multiple:

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

Polynôme de Wilkinson:

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 Cette question est restée dans le bac à sable pendant environ 3 mois. Si vous pensez que cela devait être amélioré avant de publier, visitez le Sandbox et commentez les autres questions proposées avant qu'elles ne soient publiées sur Main.

Mathematica, 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&]]

Cette solution met en œuvre la méthode Durand – Kerner pour la résolution de polynômes. Notez que ce n'est pas une solution complète (comme cela sera montré ci-dessous) car je ne peux pas encore gérer le polynôme de Wilkinson avec la précision spécifiée. D'abord une explication de ce que je fais:

#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}]&: Ainsi, la fonction calcule pour chaque indice ila prochaine approximation de Durand-Kerner. Cette ligne est ensuite encapsulée dans un tableau et appliquée à l'aide d'un NestWhile aux points d'entrée générés par Table[(0.9+0.1*I)^i,{i,l}]. La condition sur NestWhile est que la variation maximale (sur tous les termes) d'une itération à la suivante est supérieure à la précision spécifiée. Lorsque tous les termes ont moins changé, NestWhile se termine et Re@Selectsupprime les zéros qui ne tombent pas sur la ligne réelle.

Exemple de sortie:

> 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.}

Comme vous pouvez probablement le voir, lorsque le degré augmente, cette méthode commence à rebondir autour des valeurs correctes, sans jamais vraiment s'y retrouver complètement. Si je fixe la condition d'arrêt de mon code à quelque chose de plus strict que "d'une itération à l'autre, les suppositions n'ont changé que d'epsilon", alors l'algorithme ne s'arrête jamais. Je suppose que je devrais simplement utiliser Durand-Kerner comme entrée dans la méthode de Newton?