Trouver les vraies racines d'un polynôme

24

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

Peter Taylor
la source
"Vous devez être en mesure de justifier pourquoi votre programme fonctionnera toujours" Juste au cas où quelqu'un aurait besoin d'aide pour cela
Dr. belisarius
@belisarius, ??
Peter Taylor
3
a été conçu comme une blague :(
Dr belisarius
Je sais que c'est un vieux défi, alors ne vous sentez pas obligé de répondre si vous n'avez pas envie de le rouvrir. (a) Peut-on écrire une fonction, ou seulement un programme complet? (b) Dans le cas où nous pouvons écrire une fonction, pouvons-nous supposer que l'entrée utilise un type de données pratique, par exemple, Python fractions.Fraction(un type rationnel)? (c) Devons-nous gérer des polynômes de degré <1? (d) Peut-on supposer que le coefficient dominant est 1?
Ell
(e) En ce qui concerne les polynômes à racines répétées, il convient de faire une distinction entre les racines de multiplicités impaires et paires (les cas de test n'ont que des racines de multiplicités impaires.) Alors que les racines de multiplicité impaire ne sont pas trop difficiles à traiter, je '' Je ne sais pas à quel point il est tangible de gérer correctement des racines de multiplicité paire, d'autant plus que vous ne spécifiez qu'une marge d'erreur pour les valeurs des racines, pas pour leur existence. (...)
Ell

Réponses:

8

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: code au format mathématique

#[[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?

Kaya
la source
Durand-Kerner a également des problèmes potentiels avec plusieurs racines. (La méthode de Newton pourrait ne pas aider beaucoup non plus - le polynôme de Wilkinson est spécifiquement choisi pour être mal conditionné).
Peter Taylor
Vous avez tout à fait raison: j'ai abandonné cette ligne de conduite après avoir zoomé sur Wilkinson près de x = 17, c'est un gâchis absolu. Je crains de devoir opter pour une solution symbolique basée sur Groebner pour obtenir beaucoup plus de précision.
Kaya