Jak mohu získat řešení složité rovnice?

Úvod

aktualizace 6/2014

Původně jsem našel tři řešení; nyní je to šest. Myslel jsem, že vím, že jsou tři, protože racionalizace rovnice vede k polynomu stupně 24 a NSolve najde 24 kořenů a poté odstraní kořeny, které nebyly nulami původní rovnice. Ukázalo se, že rovnice byla problematická z hlediska numeriky a chyběly mi tři.

další modifikace:

  • zdá se přirozenější použít Rationalize k převodu koeficientů na přesná čísla v tomto případě než SetPrecision.

  • podobně jsem nahradil First @ eq Subtract @@ eq.

Racionalizovat rovnice.

Jeden by měl mít na paměti, že podle dokumentace na NSolve

NSolve pojednává především o lineárních a polynomiálních rovnic.

rovnice OP je algebraická rovnice, kterou lze proměnit v polynomiální rovnici, kterou NSolve pak lze snadno vyřešit. Protože NSolve to nedělá samo o sobě, musíme racionalizovat ” ručně.”Racionalizace má tendenci zavádět cizí řešení, takže musíme zkontrolovat odpovědi vrácené NSolve. Poněkud ošklivé koeficienty eq dělat to těžší, protože v rationalizating eq, round-off error vede k některé z náměstí-kořeny neruší, jak by měly. Můžeme to opravit nastavením přesnosti eq na Infinity.

můžeme najít na náměstí-kořeny ve výrazu expr s

DeleteDuplicates @ Cases, Infinity]

Pak pro každý square-root, můžeme vynásobit expr výraz s square-root nahrazen jeho negativní a zjednodušit.

níže je uveden výsledek použití metody na funkci získanou nahrazením Equal v eq Subtract Apply (@@). Můžeme vyřešit rovnici přesně.

eqExact = Rationalize;rationalized = Simplify @ Fold &, eqExact, DeleteDuplicates@Cases, Infinity]];rootsRatExact = Solve;rootsRat = NSolve(* 24 roots {{ζ -> -24559.6 + 24559.6 I}, <<22>>, {ζ -> 24559.6 - 24559.7 I}}*)

mezi oběma řešeními jsou některé významné rozdíly:

(ζ /. rootsRat) - (ζ /. rootsRatExact) // Abs // Max(* 3.15226*10^-10*)

výběr kořenů dané rovnice

aktualizace: Zde jsem vynechal některé nuly.

použil jsem Pick k výběru kořenů, u kterých má původní rovnice eq malou hodnotu, menší než nějaká malá prahová hodnota, řekněme 10^-1. Později můžeme zkontrolovat, že každý je kořen, ale nekontroloval jsem jiné kořeny. Původní kořeny byly:

rootsEq = Pick < 10^-1 /. rootsRat](* {{ζ -> -3.78042 - 5.50655 I}, {ζ -> -3.20562 + 5.39914 I}, {ζ -> 6.98478 + 0.493405 I}}*)

Tyto odpovídají kořeny 7, 9 a 20 v rootsRat:

Position < 10^-1 &)](* {{7}, {9}, {20}}*)

Když se podívám na přesnou rovnici pro nuly s PossibleZeroQ, mám víc kořenů:

Position(* {{1}, {7}, {9}, {14}, {20}, {21}}*)

trvá celou věčnost na tyto potenciální kořeny, nebo alespoň déle, než jsem byl ochoten čekat.]

je Zajímavé, že tyto další kořeny jsou spojeny s rozdíly mezi výsledky vrácené NSolve a Solve:

Position, _?Positive](* {{1}, {2}, {3}, {14}, {21}, {22}, {23}, {24}}*)

nechme naše nové kořeny být následující a můžeme podívejte se na ně níže:

rootsEqNew = Pick];rootsEqNew // N(* {{ζ -> -24559.6 + 24559.6 I}, {ζ -> -3.78042 - 5.50655 I}, {ζ -> -3.20562 + 5.39914 I}, {ζ -> -0.0832786 - 722.827 I}, {ζ -> 6.98478 + 0.493405 I}, {ζ -> 722.642 - 0.100823 I}}*)

Kontroly

nejprve se na to nevypadá moc dobře pro nové přírůstky:

eqExact /. rootsEqNew // N(* { 0. + 9.7614*10^12 I, (* {1} *) -1.49012*10^-8 + 3.91155*10^-8 I, (* {7} *) -1.39698*10^-8 - 3.63216*10^-8 I, (* {9} *) -2. + 1536. I, (* {14} *) 1.49012*10^-8 + 1.11759*10^-8 I, (* {20} *) 0. + 1024. I} (* {21} *)*)

Problém je v hodnocení rovnice s přesnosti strojů.

N(* N::meprec warning about $MaxExtraPrecision being reached *)(* {0``69.62973568978663 + 0``69.69302870899077 I, 0``90.5174054423328 + 0``90.55817837498498 I, 0``90.50250822096415 + 0``90.54414468499085 I, 0``80.1824915073549 + 0``79.76650578675965 I, 0``90.47483650216002 + 0``90.49782363232914 I, 0``80.17292602755023 + 0``79.76710897249409 I}*)

varování se nezdá být významné. Při nárazu do $MaxExtraPrecision na 2000 se stále vytváří varování a nuly s přesností kolem 20202040. Pokud jsou nuly, N pravděpodobně bude vždy produkovat varování, protože nemůže mít Precision (jiné než MachinePrecision).

Napsat komentář

Vaše e-mailová adresa nebude zveřejněna.