Temukan akar polinomial yang sebenarnya

24

Tulislah program mandiri yang ketika diberi polinomial dan terikat akan menemukan semua akar nyata polinomial itu menjadi kesalahan absolut tidak melebihi batas.

Kendala

Saya tahu bahwa Mathematica dan mungkin beberapa bahasa lain memiliki solusi satu-simbol, dan itu membosankan, jadi Anda harus tetap berpegang pada operasi primitif (penambahan, pengurangan, perkalian, pembagian).

Ada fleksibilitas tertentu pada format input dan output. Anda dapat mengambil input melalui argumen stdin atau baris perintah dalam format apa pun yang masuk akal. Anda dapat mengizinkan titik apung atau mengharuskan beberapa representasi bilangan rasional digunakan. Anda dapat mengambil terikat atau kebalikan dari terikat, dan jika Anda menggunakan floating point Anda dapat mengasumsikan bahwa terikat tidak akan kurang dari 2 ulp. Polinomial harus dinyatakan sebagai daftar koefisien monomial, tetapi mungkin besar atau sedikit-endian.

Anda harus dapat membenarkan mengapa program Anda akan selalu bekerja (masalah numerik modulo), meskipun itu tidak perlu untuk menyediakan bukti lengkap sebaris.

Program harus menangani polinomial dengan akar berulang.

Contoh

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

Masukan bisa berupa misalnya

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

Keluaran bisa misalnya

-1.41 1.42

tapi tidak

-1.40 1.40

karena memiliki kesalahan absolut sekitar 0,014 ...

Uji kasus

Sederhana:

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

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

Banyak root:

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

Polinomial 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 Pertanyaan ini ada di Sandbox selama kurang lebih 3 bulan. Jika Anda merasa perlu ditingkatkan sebelum memposting, kunjungi Sandbox dan komentari pertanyaan lain yang diajukan sebelum diposkan di Main.

code-golf math polynomials Peter Taylor
sumber
"Anda harus dapat membenarkan mengapa program Anda akan selalu bekerja" Untuk berjaga-jaga jika seseorang membutuhkan bantuan dengan itu
Dr. belisarius
@belisarius, ??
Peter Taylor
3
dimaksudkan sebagai lelucon :(
Dr. belisarius
Saya tahu ini adalah tantangan lama, jadi jangan merasa wajib untuk menjawab jika Anda tidak ingin membukanya lagi. (a) Bisakah kita menulis fungsi, atau hanya program penuh? (B) Dalam hal kita dapat menulis fungsi, dapatkah kita mengasumsikan bahwa input menggunakan beberapa tipe data yang nyaman, misalnya, Python fractions.Fraction(tipe rasional)? (c) Apakah kita harus menangani polinomial derajat <1? (d) Bisakah kita mengasumsikan bahwa koefisien yang memimpin adalah 1?
Ell
(e) Berkenaan dengan polinomial dengan akar berulang, ada baiknya membuat perbedaan antara akar ganjil dan genap genap (kasus uji hanya memiliki akar gandakan multiplikasi). Sementara akar gandakan tidak terlalu sulit untuk ditangani, saya ' Saya tidak yakin seberapa nyata menangani akar dengan benar bahkan multiplisitas secara numerik, terutama karena Anda hanya menentukan margin kesalahan untuk nilai-nilai root, bukan untuk keberadaannya. (...)
Ell

Jawaban:

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

Solusi ini mengimplementasikan metode Durand-Kerner untuk menyelesaikan polinomial. Perhatikan bahwa ini bukan solusi lengkap (seperti yang akan ditunjukkan di bawah) karena saya belum dapat menangani Polinomial Wilkinson dengan presisi yang ditentukan. Pertama, penjelasan tentang apa yang saya lakukan: kode dalam format Mathematica

#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}]&: Dengan demikian fungsi menghitung untuk setiap indeks iperkiraan Durand-Kerner berikutnya. Kemudian baris ini dienkapsulasi dalam Tabel dan diterapkan menggunakan NestWhile ke titik input yang dihasilkan oleh Table[(0.9+0.1*I)^i,{i,l}]. Kondisi pada NestWhile adalah bahwa perubahan maksimum (atas semua persyaratan) dari satu iterasi ke yang berikutnya lebih besar dari presisi yang ditentukan. Ketika semua istilah telah berubah kurang dari ini, NestWhile berakhir dan Re@Selectmenghapus nol yang tidak jatuh pada garis nyata.

Contoh output:

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

Seperti yang mungkin Anda lihat, ketika tingkat tumbuh lebih tinggi metode ini mulai bangkit di sekitar nilai-nilai yang benar, tidak pernah benar-benar pulang sepenuhnya. Jika saya mengatur kondisi berhenti kode saya menjadi sesuatu yang lebih ketat daripada "dari satu iterasi ke berikutnya tebakan diubah oleh tidak lebih dari epsilon" maka algoritma tidak pernah berhenti. Saya kira saya harus menggunakan Durand-Kerner sebagai input untuk metode Newton?

Kaya
sumber
Durand-Kerner juga memiliki masalah potensial dengan beberapa akar. (Metode Newton mungkin tidak banyak membantu juga - polinomial Wilkinson secara khusus dipilih untuk dikondisikan buruk).
Peter Taylor
Anda benar: Saya mengabaikan tindakan setelah memperbesar pada Wilkinson's dekat x = 17, ini benar-benar berantakan. Saya khawatir bahwa saya harus mencari solusi simbolis dengan basis Groebner untuk mendapatkan akurasi lebih.
Kaya