Zjednodušenie racionálnej funkcie.

Tak jednoduchy pripad je ten, kdy mas polynomy jako algebraickou extenzi celych cisel (pripadne racionalnich, protoze je muzes snadno prevest na cela). V tom pripade aplikuj Eukliduv algoritmus https://en.wikipedia.org/wiki/Polynomial_greatest…

Pro obecne reseni s algebraickou extenzi nad "realnymi" cisly si musis pohrat s presnosti. Realna cisla na pocitaci nejdou reprezentovat, takze taky budes aplikovat Eukliduv algoritmus a divat se na to jako na aproximaci. Jen pouzij trik - vzdycky pred dalsim krokem si uprav oba polynomy do monickeho tvaru.

Nemoznost najit koreny polynomu vyssich stupnu znamena, ze je neumi najit ani pocitac, vsechno jsou jen aproximace. Tim padem muzes narazet na okrajove pripady, kdy algoritmus selze a s tim se neda nic delat.

Ďakujem za info, prečítal som si o Euklidovom algoritme a ten, že funguje pekne pri celých číslach, čo nie je môj prípad. Aj keby som to nejak obišiel, vždy tam bude tá numerická nepresnosť, a to bude aj pri hľadaní koreňov číže si nepomôžem.

Ostanem pri metóde hľadania koreňov ako doteraz, len budem musieť vylepšiť moju pôvodnú iteračnú metódu (laguerre). Keď som ju písal a testoval pred niekoľkými rokmi, vykazovala vynikajúcu presnosť pre celé čísla (koreňov) a celkom obstojne pre komplexné. Rovnako jej rýchlosť konvergenie je vynikajúca. Najhoršie pre túto metódu je ale ak vstupné koeficienty polynómu sú iracionálne čísla o veľkej dlážke digitov. Na tom to jej presnosť prudko padá, ba niekedy aj zlyhá. No a práve takýto prípad musí metóda vyriešiť aspoň na 4-5 miest ak chcem korene porovnávať v čitateli a menovateli. Kvôli náročnosti výpočtu a jeho presnosti som obmedzil výstupný polynóm na max. 6 rád. Nájsť presne korene väčšieho rádu, by bolo pomerne komplikované alebo by mohli mať chybu už na 2 digite, to človek vopred nevie (tým by som polynómy nemohol skrátiť).

Testoval som medzi časom rôzne algoritmy, a porovnával výsledky s Octave a môžem povedať, že ich interný algoritmus nie je úplne dokonalý čo ma nemilo prekvapilo. Tiež trpí pomerne veľkou chybou pri mojom zadaní polynómu. Ale našiel som aj lepší algoritmus skrytý v Matlabe. Neviem sa dočítať akú metódu používajú, ale ich algoritmus je fakt vynikajúci, lepší ako výsledok Laguerre a Newton a metódy založené na derivácií polynomov. Ich metóda nájde tieto korene s presnosťou na min na 4 miesta (niekedy na aj na 15), čo by mi postačovalo. Možno používajú Jenkins-Traub, ale kto vie. Fakt ale je, že za tie peniaze to majú fakt odladené a na výsledok sa môže človek spoľahnúť.

Vie mi prosím niekto skúsenejší poradiť akú metódu hľadania koreňov by som mohol aplikovať pre lepšie výsledky?

Matlab i Octave konstruují matici, jejíž charakteristický polynom se shoduje se zadaným polynomem a aplikují na ni QR rozklad. Z definice determinantu plyne, že diagonála matice R obsahuje hledané kořeny.

Trvalo mi asi 5 minut naprogramovat si Euklida - aplikace na ten příklad, co jsi posílal, dává společného dělitele \( [1.0, -0.785893111668692, 0.3678794411714336] \) s kořeny

• (0.392946555834346+0.46203078407110365j)

• (0.392946555834346-0.46203078407110365j)

Podle toho, co píšeš v obou příspěvcích asi nemáš moc dobré znalosti ani z algoritmizace, ani z matematiky, co? Jinak bys nemohl tvrdit, že Euklidův algoritmus bude mít horší přesnost než přímé hledání kořenů.

Takže takhle: mám dva polynomy stupně N

• Euklidův algoritmus bude potřebovat nejvýše N-1 kroků

• v každém kroku normalizujeme do monického tvaru, to je násobení skalárem

• a sčítáme oba polynomy, to je sčítání skalárů

Výsledné koeficienty se budou složené maximálně z 2N operací násobení a N operací sčítání. To je 3N operací - např. pro polynomy 10. stupně mám 30 operací, \(log_2{ 30} =5\) a ve standardu IEEE 754 dostaneš výsledky nejhůř na 11 desetinných míst.

Poznám vysokoškolskú základnú matematiku (bude to už hodne dlhá doma čo som ju skončil) ale nepoznám samozrejme všetky matematické algoritmy, nie som matematik ,preto píšem sem aby som sa poučil od skúsenejších. Som elektrotechnik a program ktorý píšem nie je v Matlabe ale v Delphi, som inžinier nie akademik. Výpočet polynómov je jedna funkcia z ďalších 50, čo som musel naprogramovať aj s grafmi, preto som sa pochopitelne nemohol venovať úplne všetkým root-finding algoritmom.

Toľko o mne, ale k veci... Euklida som nikdy neprogramoval, preto ho nepoznám čo od neho čakať. Budem rád ak mi ho pomôžete pochopiť. Rád sa niečo naučím na staré kolená. Ja sa to potom pokúsim naprogramovať. Ak to bude fungovať potom nemusím programovať novú presnejšiu metódu, čo bude úplne super. Môžeš mi prosím napísať tvoj kód, pretože čo som hľadal na internete všade sa píše že to funguje len s celými číslami. Dakujem.

Nikdo nemůže znát všechno, v tom nevidím žádný problém. Ale když používáš algoritmy, tak jim to chce rozumět, aby byly použitelné v praxi.

Skoro žádný z z numerických algoritmů, které najdeš na webu, se nepoužívá v té formě, jak je tam uvádí. Ještě za tím bývají triky, které řeší numerickou nestabilitu nebo se používají speciální instrukce procesoru, které zmenšují chybu.

Už je to hodně let, co jsem se v tom vrtal, ale Matlab používá MKL - numerickou knihovnu od Intelu, zatímco Octave používá OpenBLAS. Už jenom tohle může být důvodem přesnějších výsledků v Matlabu, i když budou používat stejný algoritmus.

Tady máš největšího společného dělitele (GCD) pro celá čísla.

https://gist.github.com/coells…

Je to jednoduchý algoritmus a doporučuju prostudovat si na wikipedii jeho popis.

GCD pro polynomy je prakticky stejný, akorát srovnávám koeficienty u stejných mocnin.

https://gist.github.com/coells…

Navíc potřebuju udělat dvě věci:

• normalizovat do monického tvaru (koeficient nejvyšší mocniny je 1) kvůli stabilitě

• odečítat polynomy a pak vyhazovat vedoucí nulové koeficienty s odpovídající přesností (vybral jsem 1e-12)

No a proč jsem se ptal na tu matematiku. Na polynomy se můžeme dívat jako na extenzi tělesa reálných/komplexních čísel, koeficienty jednotlivých mocnin definují prvek konečně-dimenzionálního vektorového prostoru a GCD bude vytvářet projekci do vlastního podprostoru. A právě díky té konečné dimenzi se GCD bude chovat dobře bez ohledu na to, že koeficienty nejsou celá čísla.

Když budeš mít koeficienty na stejné škále, bude to fungovat, ale možná bude nutné doladit pro specifické případy (dělení nulou, atd). Snad je Python ok.

Ďakujem za informácie o Euklidovom algoritme pri použití na polynómy s reálnymi číslami. Presne to rieši môj problém. Je vynikajúce, že nemusím tie korene hľadať iteračnou metódou a separovať ich. Ja som doteraz o tomto spôsobe nepočul. Klobúk dole nad tým ako to niekto takto vymyslel.

Preštudoval som si tvoj algoritmus. Chvíľu mi teda trvalo pochopiť samotný Python jazyk (nikdy som s ním nerobil, asi som už na tieto novodobé jazyky starý) ale za pomoci internetu a môjho známeho som to pochopil. Dnes som napísal testovací kód v mojom Delphi. Je asi 4x dlhší ale hlavne že funguje.

Ešte raz sa ti chcem Tomáš poďakovať za matematickú pomoc. Je vidno, že si v tom doma. A ja som rád že som sa niečo nové naučil. Prajem ti veľa zdravia a úspechov!