Il modo migliore per risolvere equazioni parametriche non lineari stupidamente complicate con riduci / risolvi ecc

Aug 16 2020

Ho una funzione parametrica non lineare che è letteralmente un incubo. So che le radici esistono, sono reali e i due parametri p,esono entrambi positivi. Quello che mi aspettavo da matematica era di ottenere una soluzione (sotto forma di radice, nessuna forma chiusa) ma anche lasciando che il programma funzionasse tutta la notte mi sono arreso. Non riesco a capire se sono io a inquadrare il problema in modo non efficiente, è il mio PC che necessita di aggiornamenti seri o il problema che è semplicemente troppo difficile per metodi come Riduci o Risolvi. Se il caso è il secondo, immagino di essere condannato ... Qualche suggerimento verso gli altri due? Grazie per l'aiuto.

I miei tentativi e l'equazione:

f[x_]:=(1/(8 (p^2+x^2)^3))p^2 (-2 p^6+p^5 (4-8 x)+2 p^3 (3-8 x) x^2-6 p x^4+p^4 (80000+2 x-9 x^2)+2 p^2 x (40000+60000 x+x^2-5 x^3)-3 x^3 (-80000+40000 x+x^3)+(4 Sqrt[10] e x (p^2+x^2)^2 (2 p x^3+p^4 (-2+3 x)+2 p^3 x (-3+4 x)+x^2 (-80000+40000 x+x^3)+2 p^2 (40000-60000 x-x^2+2 x^3)))/Sqrt[-e p^2 (-1+x) x^2 (p^2+x^2)^2 (-40000+p^2+2 p x+x^2)])
Reduce[f[x]==0 && x>=0 &&p>=0 && e>=0,x,Reals] (*stuck running*)
Solve[f[x]==0 && x>=0 &&p>=0 && e>=0,x,Reals] (*stuck running*)

Risposte

2 CarlWoll Aug 18 2020 at 04:31

Se sei soddisfatto di una risposta approssimativa, puoi provare a utilizzare NDSolveValue. La tua funzione:

f[x_] := (1/(8 (p^2+x^2)^3))p^2 (-2 p^6+p^5 (4-8 x)+2 p^3 (3-8 x) x^2-6 p x^4+p^4 (80000+2 x-9 x^2)+2 p^2 x (40000+60000 x+x^2-5 x^3)-3 x^3 (-80000+40000 x+x^3)+(4 Sqrt[10] e x (p^2+x^2)^2 (2 p x^3+p^4 (-2+3 x)+2 p^3 x (-3+4 x)+x^2 (-80000+40000 x+x^3)+2 p^2 (40000-60000 x-x^2+2 x^3)))/Sqrt[-e p^2 (-1+x) x^2 (p^2+x^2)^2 (-40000+p^2+2 p x+x^2)])

Per utilizzare NDSolveValue, abbiamo bisogno di conoscere una condizione al contorno. Ad esempio, ecco il valore di xquando pè 1:

x1 = x /. Block[{p=1}, First @ Solve[f[x] == 0, x]]

Radice [256006399839996 + 1023948800640 e + (255942399200020-3071999998080 e) # 1 + (511955203840004 + 2304217598880 e) # 1 ^ 2 + (1279846402079976-2048025605760 e) # 4 ^ 3 + (-960 +2799842) (192129617159615 + 4095897593600 e) # 1 ^ 5 + (-384151995680463 - 512486397600 e) # 1 ^ 6 + (-3455678391520375 + 2047846414080 e) # 1 ^ 7 + (2880427177039798 - 332788480 + e -59) # 1 ^ 624 - 576 1024102385280 e) # 1 ^ 9 + (43199520578 + 256064008800 e) # 1 ^ 10 + (-14402879738 - 25598080 e) # 1 ^ 11 + (-359829 + 12802240 e) # 1 ^ 12 + 360087 # 1 ^ 13 + ( 9 + 160 e) # 1 ^ 14 + 9 # 1 ^ 15 e, 1]

Ora possiamo usare NDSolveValue:

sol = NDSolveValue[
    {
    D[f[x[p,e]]==0, p], x[1, e] == x1},
    x,
    {p,.1,100},
    {e,.1,10000},
    MaxStepFraction->.0005,
    PrecisionGoal->10
]; //AbsoluteTiming

{19.2292, Null}

Controlla alcuni campioni casuali:

Block[{p = 50, e = 200}, f[sol[p, e]]]
Block[{p = 10, e = 2000}, f[sol[p, e]]]

6.42413 * 10 ^ -9

8.0893 * 10 ^ -9

Visualizzazione:

Plot3D[sol[p,e], {p,.1,100}, {e,.1,10000}]

2 MichaelE2 Aug 17 2020 at 04:15

FWIW, ecco una mezza soluzione: prendi il numeratore, razionalizzalo in modo che sia un polinomio e trova le radici. Quello che resta da fare è selezionare quelli che sono positivi quando pe esono positivi. Questo passaggio richiede molto tempo (ammesso che sia possibile), tranne quando vengono forniti valori numerici specifici a pe e.

num = Simplify[ff, x >= 0 && p >= 0 && e >= 0] // Together // 
      Numerator // Simplify[#, x >= 0 && p >= 0 && e >= 0] & // 
    FactorList // #[[-1, 1]] & // Simplify
(*
p Sqrt[-e (-1 + x) (-40000 + p^2 + 2 p x + x^2)] (2 p^6 + 6 p x^4 + 
    p^5 (-4 + 8 x) + 2 p^3 x^2 (-3 + 8 x) + 
    p^4 (-80000 - 2 x + 9 x^2) + 3 x^3 (-80000 + 40000 x + x^3) + 
    2 p^2 x (-40000 - 60000 x - x^2 + 5 x^3)) - 
 4 Sqrt[10]
   e (2 p x^5 + 4 p^3 x^3 (-1 + 2 x) + p^6 (-2 + 3 x) + 
    2 p^5 x (-3 + 4 x) + p^2 x^3 (-80000 - 2 x + 5 x^2) + 
    x^4 (-80000 + 40000 x + x^3) + 
    p^4 (80000 - 120000 x - 4 x^2 + 7 x^3))
*)

Verifica che ci siano due termini (il primo ovviamente contenente il radicale):

Length@num
(*  2  *)
rat = num*MapAt[-# &, num, 1] // Expand // Simplify
(*
e (p^2 (-1 + x) (-40000 + p^2 + 2 p x + x^2) (2 p^6 + 6 p x^4 + 
      p^5 (-4 + 8 x) + 2 p^3 x^2 (-3 + 8 x) + 
      p^4 (-80000 - 2 x + 9 x^2) + 3 x^3 (-80000 + 40000 x + x^3) + 
      2 p^2 x (-40000 - 60000 x - x^2 + 5 x^3))^2 + 
   160 e (2 p x^5 + 4 p^3 x^3 (-1 + 2 x) + p^6 (-2 + 3 x) + 
      2 p^5 x (-3 + 4 x) + p^2 x^3 (-80000 - 2 x + 5 x^2) + 
      x^4 (-80000 + 40000 x + x^3) + 
      p^4 (80000 - 120000 x - 4 x^2 + 7 x^3))^2)
*)
Solve[rat == 0, x]
(* <15 Root objects> *)

Queste radici contengono soluzioni estranee e sembra che sia necessario sostituire valori numerici affinché i parametri funzionino con esse. Se un tale approccio è utile, forse sarebbe meglio sostituire i parametri f[x]e trattare l'equazione risultante f[x] == 0.