Per la messa a punto del Runge Kutta di Dormand Prince:
• http://www.elegio.it/calcolatrice/nuova-dormand-prince.html
• http://www.elegio.it/calcolatrice/asteroide-solare.html
• http://www.elegio.it/calcolatrice/libreria-dormand-prince.js.html
• http://www.elegio.it/calcolatrice/libreria-dormand-prince-provo.html
• http://www.elegio.it/calcolatrice/libreria-dormand-prince_libhtml.html con spiegazioni dei casi test...
• http://www.elegio.it/calcolatrice/dormand-prince-mathematichese.pdf
• http://www.elegio.it/mc2/rk/ == altri tipi di Runge Kutta

Dormand Prince in Javascript

Trascrivo qui i risultati di quattro casi prova di crescente difficoltà

// Il risultato finale dovrebbe essere [1,7.38905609893065]

...100

...101

...102

// Il risultato finale dovrebbe essere [-1,0]

...200

...201

...202

// Il risultato finale dovrebbe essere [146079760576,0,0,30500] ossia la distanza dal Sole dovrebbe essere circa 146 miliardi di metri e la velocità tangenziale, lungo y dovrebbe essere quella data all'inizio ossia 30500. L'orbita dovrebbe essere completata e la durata dell'orbita ossia il tempo finale è stato determinato usando la soluzione esatta del moto ellittico. Usando le formule si è fissata la lunghezza del transitorio ossia un anno ovvero un intero periodo di rivoluzione che è risultato essere di 31556606 secondi...31 milioni circa...

Notare che questo è un test più impegnativo dei precedenti dato che richiede la soluzione di 4 equazioni differenziali.

...300

...301

...302

// Un sistema a tre corpi non ammette una soluzione analitica per cui il solo test possibile è il confronto tra calcoli fatti a diversa tolleranza degli errori...

...400

...401

...402

...403

//

• Uso questa funzione:

  ..

• Questi sono i casi prova:

  ..

Dormand Prince in maximese

Qui sono le funzioni originarie che ho realizzato in Maxima programmandole in questo discutibile modo...

nucleode():=block(
   ha21:34272564480*h,
   k2:fder(x+34272564480*h,y+ha21*k1,costa),
   ha31:12852211680*h,
   ha32:38556635040*h,
   k3:fder(x+51408846720*h,y+ha31*k1+ha32*k2,costa),
   ha41:167554759680*h,
   ha42:-639754536960*h,
   ha43:609290035200*h,
   k4:fder(x+137090257920*h,y+ha41*k1+ha42*k2+ha43*k3,costa),
   ha51:505965644800*h,
   ha52:-1987087872000*h,
   ha53:1683278643200*h,
   ha54:-49833907200*h,
   k5:fder(x+152322508800*h,y+ha51*k1+ha52*k2+
      ha53*k3+ha54*k4,costa),
   ha61:487745760600*h,
   ha62:-1843448544000*h,
   ha63:1526229734400*h,
   ha64:47708967600*h,
   ha65:-46873096200*h,
   hc6:171362822400*h,
   k6:fder(x+hc6,y+ha61*k1+ha62*k2+ha63*k3+
      ha64*k4+ha65*k5,costa),
   ha71:15619007250*h,
   ha73:76982400000*h,
   ha74:111564337500*h,
   ha75:-55243291950*h,
   ha76:22440369600*h,
   k7:fder(x+171362822400*h,
      y+ha71*k1+ha73*k3+ha74*k4+ha75*k5+ha76*k6,costa)
   )$

Utilizzata in questa funzione:

odedopri(x0,y0,x1,costa,peso,tol,maxiter):=block(
  [ha21,ha31,ha32,ha41,ha42,ha43,ha51,ha52,ha53,ha54,
  ha61,ha62,ha63,ha64,ha65,ha71,ha73,ha74,ha75,ha76,
   e1,e3,e4,e5,e6,e7,
   i,j,ifa,x,y,ny,h,hmax,hmin,k1,k2,k3,k4,k5,k6,k7,
  verrore,erromax],
  x:float(x0),
  y:float(y0),
  ny:length(y),
  h:float(10*(x1-x0)/(maxiter*171362822400)),
  hmax:h,
  hmin:h/100,
  erromax:0,
  amen:false,
  e1:211228479,
  e3:-728766720,
  e4:6336854370,
  e5:-8716163841,
  e6:7180918272,
  e7:-4284070560,
  k1:fder(x,y,costa),
  for i:1 thru maxiter do(
     ifa:i,
     nucleode(),
     verrore:peso*abs(e1*k1+e3*k3+e4*k4+e5*k5+e6*k6+e7*k7),
     erromax:0,
     for j:1 thru ny do(
         if verrore[j]>erromax then erromax:verrore[j]
         ),
     if tol>erromax/171362822400 then avanti:true
     else ( if hc6>hmin then ( h:h/2,
              if debug then 
                print(["riduce",h,"i=",i,"erromax/171362822400=",
                    erromax/171362822400,"tol=",tol]),
               avanti:false ) 
             else avanti:true,
           amen:false
           ),
     if avanti then (
         x:x+hc6,
         y:y+ha71*k1+ha73*k3+ha74*k4+ha75*k5+ha76*k6,
         k1:float(k7),
         if amen then return(x),
         if hmax>h then h:1.05*h
         ),
     if x+h*171362822400>x1 
         then (h:(x1-x)/171362822400,amen:true)
     ),
  [amen,x,y,ifa] 
  )$

Bibliografia

• http://it.wikipedia.org/wiki/Problema_di_Cauchy
• http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method
• http://en.wikipedia.org/wiki/Butcher_tableau#Explicit_Runge.E2.80.93Kutta_methods

Panoramica dell'algoritmo

Ovviamente non voglio mettere in competizione questo solutore di equazioni differenziali alle derivate ordinarie con le funzioni realizzate da competenti in Maxima o in Mathematica 9 ma penso che questa versione del Runge Kutta di Dormand e Prince serva a ... fini didattici per capire come si può affrontare questo tipo di problemi

Insomma il mio è solo il desiderio di essere capace di camminare con le mie gambe ed evitare di annegare anche se, ovviamente, pur sapendo camminare e stare a galla, è molto meglio usare una bicicletta o una motocicletta o una automobile e, in mare, stare su una robusta e potente nave...

Giampaolo Bottoni gpbottoni@gmail.com http://www.alumni.polimi.it/it/Wall ( ing. nucleare 1972 )