Encontre raízes reais de um polinômio

24

Escreva um programa independente que, ao receber um polinômio e um limite, encontre todas as raízes reais desse polinômio em um erro absoluto que não exceda o limite.

Restrições

Eu sei que o Mathematica e provavelmente algumas outras línguas têm uma solução de um símbolo, e isso é chato, então você deve se ater às operações primitivas (adição, subtração, multiplicação, divisão).

Há certa flexibilidade nos formatos de entrada e saída. Você pode receber informações via stdin ou argumentos da linha de comando em qualquer formato razoável. Você pode permitir ponto flutuante ou exigir que alguma representação de números racionais seja usada. Você pode pegar o limite ou o recíproco do limite e, se estiver usando ponto flutuante, pode assumir que o limite não será menor que 2 ulp. O polinômio deve ser expresso como uma lista de coeficientes monomiais, mas pode ser grande ou pequeno endiano.

Você deve ser capaz de justificar por que seu programa sempre funcionará (problemas numéricos de módulo), embora não seja necessário fornecer provas completas em linha.

O programa deve lidar com polinômios com raízes repetidas.

Exemplo

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

A entrada pode ser, por exemplo,

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

A saída pode ser, por exemplo,

-1.41 1.42

mas não

-1.40 1.40

como isso tem erros absolutos de cerca de 0,014 ...

Casos de teste

Simples:

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

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

Raiz múltipla:

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

Polinômio 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 Esta pergunta esteve na Sandbox por aproximadamente 3 meses. Se você acha que precisava melhorar antes de postar, visite a Sandbox e comente as outras perguntas propostas antes de serem publicadas na página Principal.

code-golf math polynomials Peter Taylor
fonte
"Você deve ser capaz de justificar por que seu programa irá funcionar sempre" Apenas no caso de alguém precisa de ajuda com isso
Dr. belisarius
@belisarius, ??
Peter Taylor
3
foi planejado como uma piada :(
Dr. belisarius
Sei que esse é um desafio antigo, por isso não se sinta obrigado a responder se não quiser reabri-lo. (a) Podemos escrever uma função ou apenas um programa completo? (b) Caso possamos escrever uma função, podemos assumir que a entrada usa algum tipo de dados conveniente, por exemplo, o Python fractions.Fraction(um tipo racional)? (c) Temos que lidar com polinômios de grau <1? (d) Podemos assumir que o coeficiente inicial é 1?
Ell
(e) Com relação aos polinômios com raízes repetidas, vale a pena fazer uma distinção entre raízes de multiplicidades ímpares e pares (os casos de teste têm apenas raízes de multiplicidades ímpares). Embora as raízes da multiplicidade ímpar não sejam muito difíceis de lidar, eu ' não tenho certeza de quão tangível é lidar corretamente com raízes de multiplicidade uniforme numericamente, especialmente porque você especifica apenas uma margem de erro para os valores das raízes, não para a sua existência. (...)
Ell 16/12

Respostas:

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

Esta solução implementa o método Durand-Kerner para resolver polinômios. Observe que essa não é uma solução completa (como será mostrado abaixo) porque ainda não consigo lidar com o polinômio de Wilkinson com a precisão especificada. Primeiro, uma explicação do que estou fazendo: código na formatação mathematica

#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}]&: Assim, a função calcula para cada índice ia próxima aproximação de Durand-Kerner. Em seguida, essa linha é encapsulada em uma tabela e aplicada usando um NestWhile aos pontos de entrada gerados por Table[(0.9+0.1*I)^i,{i,l}]. A condição no NestWhile é que a alteração máxima (em todos os termos) de uma iteração para a próxima seja maior que a precisão especificada. Quando todos os termos mudaram menos que isso, o NestWhile termina e Re@Selectremove os zeros que não caem na linha real.

Exemplo de saída:

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

Como você provavelmente pode ver, quando o grau aumenta, esse método começa a ricochetear nos valores corretos, nunca realmente chegando completamente. Se eu definir a condição de parada do meu código como algo mais rigoroso do que "de uma iteração para a próxima, as suposições mudaram em não mais do que epsilon", o algoritmo nunca para. Acho que devo usar Durand-Kerner como entrada para o método de Newton?

Kaya
fonte
Durand-Kerner também tem problemas em potencial com múltiplas raízes. (O método de Newton também pode não ajudar muito - o polinômio de Wilkinson é especificamente escolhido para ser mal condicionado).
Peter Taylor
Você está certo: abandonei esse curso de ação depois de ampliar o zoom de Wilkinson perto de x = 17, é uma bagunça absoluta. Preocupo-me por ter que buscar uma solução simbólica com base no Groebner para obter muito mais precisão.
Kaya