Hur kan jag få lösningen av komplicerad ekvation?
introduktion
Uppdatering 6/2014
ursprungligen hittade jag tre lösningar; nu är det sex. Jag trodde att jag visste att det fanns tre eftersom rationalisering av ekvationen resulterar i ett polynom av grad 24 och NSolve
hittar 24 rötter och eliminerar sedan rötter som inte var nollor i den ursprungliga ekvationen. Det visar sig att ekvationen var besvärlig ur numerikens synvinkel och jag missade tre.
andra ändringar:
-
det verkar mer naturligt att använda
Rationalize
för att konvertera koefficienterna till exakta siffror i detta fall änSetPrecision
. -
på samma sätt ersatte jag
First @ eq
medSubtract @@ eq
.
rationalisera ekvationen
man bör komma ihåg att enligt dokumentationen om NSolve
NSolve
behandlar främst linjära och polynomekvationer.
OP: s ekvation är en algebraisk ekvation som kan omvandlas till en polynomekvation, som NSolve
sedan enkelt kan lösa. Eftersom NSolve
inte gör detta på egen hand, måste vi rationalisera “för hand.”Rationalisering tenderar att införa externa lösningar, så vi måste kontrollera svaren som returneras av NSolve
. De något otäcka koefficienterna på eq
gör det svårare,eftersom rationalisering av eq
leder till att några av kvadratrötterna inte avbryter som de borde. Vi kan fixa det genom att ställa in precisionen eq
till Infinity
.
vi kan hitta kvadratrötterna i ett uttryck expr
med
DeleteDuplicates @ Cases, Infinity]
sedan för varje kvadratrot kan vi multiplicera expr
med uttrycket med kvadratroten ersatt av dess negativa och förenkla.
nedan är resultatet av att metoden tillämpas på funktionen som erhålls genom att ersätta Equal
i eq
med Subtract
med Apply
(@@
). Vi kan lösa ekvationen exakt.
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}}*)
det finns några signifikanta skillnader mellan de två lösningarna:
(ζ /. rootsRat) - (ζ /. rootsRatExact) // Abs // Max(* 3.15226*10^-10*)
välja rötterna för den givna ekvationen
Uppdatering: här missade jag några nollor.
jag använde Pick
för att välja rötter där den ursprungliga ekvationen eq
har ett litet värde, mindre än något litet tröskelvärde, säg 10^-1
. Vi kan kontrollera senare att var och en är en rot, men jag kollade inte efter andra rötter. De ursprungliga rötterna var:
rootsEq = Pick < 10^-1 /. rootsRat](* {{ζ -> -3.78042 - 5.50655 I}, {ζ -> -3.20562 + 5.39914 I}, {ζ -> 6.98478 + 0.493405 I}}*)
dessa motsvarar rötterna 7, 9 och 20 i rootsRat
:
Position < 10^-1 &)](* {{7}, {9}, {20}}*)
om jag kontrollerar den exakta ekvationen för nollor med PossibleZeroQ
får jag fler rötter:
Position(* {{1}, {7}, {9}, {14}, {20}, {21}}*)
tar evigt på dessa potentiella rötter, eller åtminstone längre än jag var villig att vänta.]
intressant är att dessa ytterligare rötter är kopplade till skillnaderna mellan resultaten som returneras av NSolve
och Solve
:
Position, _?Positive](* {{1}, {2}, {3}, {14}, {21}, {22}, {23}, {24}}*)
Låt oss låta våra nya rötter vara följande och vi kan kontrollera dem nedan:
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}}*)
kontroller
först, at ser inte för bra ut för de nya tilläggen:
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} *)*)
problemet är att utvärdera ekvationen med maskinprecision.
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}*)
varningen verkar inte vara signifikant. Bumping up $MaxExtraPrecision
till 2000
producerar fortfarande varningen och nollorna med noggrannhet runt 2020
till 2040
. Om de är nollor kommer N
förmodligen alltid att producera varningen, eftersom inte kan ha en
Precision
(annat än MachinePrecision
).