*/LIBRERIA DORMAND PRINCE Javascript - MAGGIO 2014.
Con estensione HTML posso ottenere persino scritte rosse! In rete:
http://www.elegio.it/calcolatrice/libreria-dormand-prince.js.html
http://www.elegio.it/calcolatrice/nuova-dormand-prince.html*/ // // Questa e' la funzione fondamentale. // Richiede in argomento fder(x,y,costa) ossia la funzione da usare // nel Runge Kutta e che deve produrre come risultato // un vettore. // Il primo argomento di fder e' il vettore delle incognite // ossia le posizioni e le velocita' // Il secondo argomento di fder e' il vettore delle variazioni // Il terzo argomento di fder e' il vettore delle costanti usate // che non vanno cambiate ma solo usate. // // Il significato degli altri argomenti e' questo: // // x0 == Valore iniziale della variabile indipendente. // y0 == Valore iniziale del vettore delle variabili dipendenti. // x1 == Valore finale desiderato della variabile indipendente. // costa == Costanti usate per calcolare le derivate prime delle // variabili dipendenti che all'inizio valgono y1. // peso == Vettori usato per equilibrare l'importanza degli errori // del calcolo. Spesso l'errore di una qualche variabile // ha importanza molto diversa dall'errore di una qualche // altra variabile. // tol == parametro usato per limitare l'errore. Piu' e' piccolo // e maggiore e' la precisione ossia minore la tolleranza // dell'errore totale ottenuto come somma degli errori // delle varie variabili pesati sul peso della loro importanza. // function odedopri(fder,x0,y0,x1,costa,peso,tol,maxiter){ var ha21,ha31,ha32,ha41,ha42,ha43,ha51,ha52,ha53,ha54, ha61,ha62,ha63,ha64,ha65,ha71,ha73,ha74,ha75,ha76; var e1,e3,e4,e5,e6,e7; var i,j,ifa,x,ny,h,hmax,hmin,verrore,erromax; var y=[0],z=[0],k1=[0],k2=[0],k3=[0],k4=[0],k5=[0],k6=[0],k7=[0]; x=x0; ny=y0.length; // alert("ny = "+ny); for(j=1;ny>j;j++)y[j]=y0[j]; h=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); // // alert("Stampo k1 \n"+k1); // for(i=1;maxiter>=i;i++){ ifa=i; ha21=34272564480*h; for(j=1;ny>j;j++)z[j]=y[j]+ha21*k1[j]; k2=fder(x+34272564480*h,z,costa); ha31=12852211680*h; ha32=38556635040*h; for(j=1;ny>=j;j++)z[j]=y[j]+ha31*k1[j]+ha32*k2[j]; k3=fder(x+51408846720*h,z,costa); ha41=167554759680*h; ha42=-639754536960*h; ha43=609290035200*h; for(j=1;ny>j;j++)z[j]=y[j]+ha41*k1[j]+ha42*k2[j]+ ha43*k3[j]; k4=fder(x+137090257920*h,z,costa); ha51=505965644800*h; ha52=-1987087872000*h; ha53=1683278643200*h; ha54=-49833907200*h; for(j=1;ny>j;j++)z[j]=y[j]+ha51*k1[j]+ha52*k2[j]+ ha53*k3[j]+ha54*k4[j]; k5=fder(x+152322508800*h,z,costa); ha61=487745760600*h; ha62=-1843448544000*h; ha63=1526229734400*h; ha64=47708967600*h; ha65=-46873096200*h; hc6=171362822400*h; for(j=1;ny>j;j++)z[j]=y[j]+ha61*k1[j]+ha62*k2[j]+ ha63*k3[j]+ha64*k4[j]+ha65*k5[j]; k6=fder(x+hc6,z,costa); ha71=15619007250*h; ha73=76982400000*h; ha74=111564337500*h; ha75=-55243291950*h; ha76=22440369600*h; for(j=1;ny>j;j++)z[j]=y[j]+ha71*k1[j]+ ha73*k3[j]+ha74*k4[j]+ha75*k5[j]+ha76*k6[j]; k7=fder(x+171362822400*h,z,costa); // erromax=0; for(j=1; ny>j; j++){verrore=peso[j]* Math.abs(e1*k1[j]+e3*k3[j]+e4*k4[j]+e5*k5[j]+ e6*k6[j]+e7*k7[j]); if(verrore>erromax)erromax=verrore;} // if(tol>erromax/171362822400){avanti=true;} else{ if(hc6>hmin){h=h/2;avanti=false;} else { avanti=true; } amen=false;} // if(avanti){ x=x+hc6; for(j=1; ny>j; j++){ y[j]=y[j]+ha71*k1[j]+ha73*k3[j]+ha74*k4[j]+ ha75*k5[j]+ha76*k6[j]; k1[j]=k7[j]}; if(amen){break}; if(hmax>h)h=1.05*h; } if(x+h*171362822400>x1){ h=(x1-x)/171362822400;amen=true}; // }; // // amen vale true se ha fatto tutto, altrimenti false. // ifa indica il numero di iterazioni fatte realmente. // return [0,amen,x,y,ifa]; } // // __________________________________________________________ // // CASI PROVA LIBRERIA DORMAND PRINCE // // Ogni caso e' basato su due funzioni. La funzione da integrare // e la funzione che fa l'inizializzazione dei dati del calcolo. // // Le seguenti sono le variabili d'ambiente utilizzate: // // var tol=1.0e-13,maxiter=250; var x_0,x_1,y_0,costa,peso; // // // Primo caso prova, elementare... // // function primafder(x,y,c){ var rr; rr=Math.exp(2*x)/2; return [0,1,rr+1.5*y[2]]; } // function ini_fder1(){ x_0=0; x_1=1; y_0=[0,0,1]; costa=[0,1]; peso=[0,1,1]; tol=1.0e-13; maxiter=250; return [0,x_0,y_0,x_1,costa,peso,tol,maxiter] } // // // Secondo caso prova leggermente piu' complicato. // // function secondafder(x,v,c){ var y,z; y=v[1]; z=v[2]; return [0,z/2,-2*y]; } // function ini_fder2(){ x_0=0; x_1=1.5*Math.PI; y_0=[0,0,2]; costa=[0,1]; peso=[0,1,1]; tol=1.0e-13; maxiter=2000; return [0,x_0,y_0,x_1,costa,peso,tol,maxiter] } // // // Terzo caso prova: moto piano ( bidimensionale ) di asteroide // attorno al Sole. // // function terzafder(t,v,c){ var x,y,vx,vy,mu,rq,rr,den,x_,y_,vx_,vy_; x=v[1]; y=v[2]; vx=v[3]; vy=v[4]; mu=c[1]; x_=vx; y_=vy; rq=x*x+y*y; rr=Math.sqrt(rq); den=1/(rr*rq); vx_=-mu*x*den; vy_=-mu*y*den; // // Il risultato e' il vettore delle quattro variabili ossia: // la derivata della posizione in ascissa, // la derivata della posizione in ordinata, // la derivata della velocita' in ascissa, // la derivata della velocita' in ordinata. // return [0,x_,y_,vx_,vy_]; } // function ini_fder3(){ // velocita' media, al perielio, all'afelio, eccentricita'... var vm,vp,va,ec; costa=[0,1.327581e20,149.61e9,30500.0]; vm=Math.sqrt(costa[1]/costa[2]); if(costa[3]>vm){vp=costa[3];va=vm*vm/vp;} else {va=costa[3];vp=vm*vm/va}; ec=(vp-va)/(vp+va); x_0=0; x_1=2*Math.PI*costa[2]/vm; y_0=[0,costa[2]*(1-ec),0,0,vp]; peso=[0,1.0e-9,1.0e-9,1.0e-4,1.0e-4]; tol=5.0e-14; maxiter=500; return [0,x_0,y_0,x_1,costa,peso,tol,maxiter]; } // //_________________________________________________________________ // // Problema di tre corpi nel piano ( Terra e Luna // attorno al Sole immobile ) // xt e yt sono l'ascissa e l'ordinata della Terra, // xl e yl sono l'ascissa e l'ordinata della Luna. // x sarebbe il tempo che pero' non interviene esplicitamente // nel calcolo delle derivate prime. // vxt e vyt sono la velocita' della Terra in ascissa e in ordinata. // vxl e vyl sono la velocita' della Luna in ascissa e in ordinata. // dunque le variabili sono in totale 8 e di queste la funzione // deve produrre le loro derivate prime. // Nel caso delle posizioni ossia le prime 4 variabili, le derivate // prime sono ovviamente le velocita'. // Nel caso delle velocita' ossia le variabili dalla quinta alla // ottava, le derivate prime sono ovviamente le forze. // // function quartafder(x,v,c){ var mu,ts,xt,yt,vxt,vyt,rqt,rrt,dent,xt_,yt_,vxt_,vyt_; var dentl,denlt,ls,xl,yl,vxl,vyl,rql,rrl,denl,xl_yl_,vxl_vyl_; // // Ricopia il vettore delle coordinate solo per rendere // piu' leggibile questa funzione. // xt=v[1]; yt=v[2]; xl=v[3]; yl=v[4]; vxt=v[5]; vyt=v[6]; vxl=v[7]; vyl=v[8]; // // mu e' la azione gravitazionale del Sole ossia la sua massa per G // ossia la costante di gravitazione universale. // ts e' il rapporto tra la massa della Terra e quella del Sole // ls e' il rapporto tra la massa della Luna e quella del Sole... // mu=c[1]; ts=c[2]; ls=c[3]; // // Ovviamente le derivate prime delle velocita' della Terra e della // Luna sono le velocita' ( bidimensionali, in ascissa e ordinata ) // xt_=vxt; yt_=vyt; xl_=vxl; yl_=vyl; // // rqt e' il quadrato della distanza della Terra dal Sole. // rqt=xt*xt+yt*yt; // // rrt e' la distanza della Terra dal Sole. // rrt=Math.sqrt(rqt); // // l'inverso del cubo della distanza tra la Terra e in Sole. // dent=1/(rrt*rqt); // // rql e' il quadrato della distanza della Luna dal Sole. // rql=xl*xl+yl*yl; // // rrl e' la distanza della Luna dal Sole. // rrl=Math.sqrt(rql); // // l'inverso del cubo della distanza tra la Luna e il Sole. // denl=1/(rrl*rql); // // rqtl e' il quadrato della distanza della Luna dalla Terra. // rqtl=(xt-xl)*(xt-xl)+(yt-yl)*(yt-yl); // // rrtl e' la distanza della Luna dalla Terra. // rrtl=Math.sqrt(rqtl); // // per velocizzare i calcoli precalcolo queste due quantita' // dentl=ls/(rqtl*rrtl); denlt=ts/(rqtl*rrtl); // // le forze sulla Terra in ascissa ed in ordinata. Per // ottenere una forza attrattiva la forza deve avere il // segno meno. // vxt_=-mu*(xt*dent+(xt-xl)*dentl); vyt_=-mu*(yt*dent+(yt-yl)*dentl); // // le forze sulla Luna in ascissa ed in ordinata // vxl_=-mu*(xl*denl+(xl-xt)*denlt); vyl_=-mu*(yl*denl+(yl-yt)*denlt); return [0,xt_,yt_,xl_,yl_,vxt_,vyt_,vxl_,vyl_]; } // _____________________________________________________________ // // Produce un vettore di dati di questo quarto test ossia: // 1== L'istante di tempo assunto come iniziale // 2== Il vettore delle condizioni iniziali ossia le posizioni // delle particelle e le velocita' delle particelle. // 3== L'istante di tempo che bisogna raggiungere. // 4== Le costanti da usare per calcolare le forze ossia // le accelerazioni. // 5== Il vettore di equilibramento delle tolleranze degli errori. // 6== Il parametro di controllo dell'errore, piccolissimo // se molto preciso. // 7== Il numero massimo delle iterazioni ammesse per fare // evolvere il sistema fino alla fine. // function ini_fder4(){ var vm,vp,va,ec; // // Le costanti usate sono: // 1== la costante gravitazionale solare ossia // la massa del sole per la G ossia costante di // gravitazione universale. // 2== Il rapporto tra la massa della Terra e la massa del Sole. // 3== Il rapporto tra la massa della Luna e la massa del Sole // e la massa della Luna e' circa un ottantesimo della massa // della terra per cui uso questo dato... // 4== La distanza iniziale della Terra dal Sole. // 5== La velocita' in metri al secondo della Terra. // costa=[0,1.327581e20, 1/333000.1, 1/(333000.1*80), 149.61e9, 30500.0]; vm=Math.sqrt(costa[1]/costa[4]); if(costa[5]>vm){vp=costa[5];va=vm*vm/vp;} else {va=costa[5];vp=vm*vm/va; }; ec=(vp-va)/(vp+va); // // x_0 e' l'istante iniziale del transitorio. // x_0=0; // // x_1 e' l'istante finale a cui vorrei arrivare nella simulazione. // x_1=2*Math.PI*costa[4]/vm; // // y_0 e' il vettore delle incognite: 4 posizioni e 4 velocita'. // y_0=[0,costa[4]*(1-ec),0,costa[4]*(1-ec),-360e6,0,vp,1100.0,vp]; // // il vettore peso controlla la tolleranza del calcolo // Il peso delle posizioni deve essere molto piccolo // rispetto ai pesi delle velocita' perche' le // posizioni si misurano in miliardi di metri e le // velocita' si misurano in decine di migliaia di metri al secondo. // peso=[0,1e-9,1e-9,1e-9,1e-9,1e-4,1e-4,1e-4,1e-4]; // // Il criterio di tolleranza piuttosto severo... // tol=0.5e-14; // // il massimo numero delle iterazioni consentite: // maxiter=2000; // // Ecco il vettore risultato usabile per inizializzare i calcoli. // Notare che tutti questi dati sono variabili globali // e non locali mentre il vettore risultato e' solo una raccolta // di queste variabili globali inizializzate da questa funzione. // return [0,x_0,y_0,x_1,costa,peso,tol,maxiter]; } // // ____________________________________________________________________ // //PER PILOTARE I CASI PROVA DELLA LIB. DORMAND PRINCE // // Devono esistere le seguenti MARCHE HTML: // // Per il primo caso prova : "x100", "x101", "x102". // Per il secondo caso prova : "x200", "x201", "x202". // Per il terzo caso prova : "x300", "x301", "x302". // Per il quarto caso prova : "x400", "x401", "x402", "x403". // function test_casi_prova_libreria(){ var obs,ob,p1,ris3,ris4,ris5,inforis1,inforis2; // // ====================================================== // Prima prova. // Il risultato finale dovrebbe essere [1,7.38905609893065] // p1=ini_fder1(); inforis1=[0,"var.indip.x", "x considerato dipendente","y","x finale","costa", "xpeso","ypeso","tolleranza","iterazioni massime"]; inforis2=[0,"Bene se true, male se false","x","xdip","y", "iterazioni fatte" ]; ob=document.getElementById("x100"); ob.innerHTML="Dati di inizializzazione del primo test "+ vedoris(p1,inforis1); ris3=odedopri(primafder,x_0,y_0,x_1,costa,peso,tol,maxiter); ob=document.getElementById("x101"); ob.innerHTML="Avendo fatto al massimo "+maxiter+" iterazioni "+ " non dovrebbe riuscire a trovare la soluzione "+ vedoris(ris3,inforis2); if(!ris3[1])ris4=odedopri(primafder,ris3[2],ris3[3],x_1, costa,peso,tol,4*maxiter); ob=document.getElementById("x102"); ob.innerHTML="Risultato finale avendo "+ " quadruplicato il massimo delle iterazioni ammesse "+ vedoris(ris4,inforis2); // // ====================================================== // // Seconda prova. // Il risultato finale dovrebbe essere [-1,0] // p1=ini_fder2(); inforis1=[0,"varindip.iniz.x","y","z","varindip.finale", "cost","ypeso","zpeso","tolleranza","maxpassi"]; inforis2=[0,"Bene se true, male se false","xfinale", "yfinale","zfinale","iterazionifatte"]; ob=document.getElementById("x200"); ob.innerHTML="Dati di inizializzazione del secondo test "+ vedoris(p1,inforis1); ris3=odedopri(secondafder,x_0,y_0,x_1,costa,peso,tol,maxiter); ob=document.getElementById("x201"); ob.innerHTML="Avendo fatto al massimo "+maxiter+" iterazioni "+ "non trova la soluzione"+ vedoris(ris3,inforis2); if(!ris3[1])ris4=odedopri(secondafder,ris3[2],ris3[3],x_1, costa,peso,tol,4*maxiter); ob=document.getElementById("x202"); ob.innerHTML="Risultato finale con il quadruplo di iterazioni"+ " che non usa se non quando servono al criterio "+ " di cambio ampiezza passo "+vedoris(ris4,inforis2); // // ============================================================== // // Terza prova. // Il risultato finale dovrebbe essere : // [146079760576.14456,0,0,30500] // p1=ini_fder3(); inforis1=[0,"t-iniziale","xaster","yaster","vxaster","vyaster", "t_finale","costaG*M","costa_amaggiore","costavmedia", "xpeso","ypeso","vxpeso","vypeso","tolleranza","Maxiterazioni"]; inforis2=[0,"Bene se true, male se false", "t_finale","xfinale","yfinale","vxfinale","vyfinale", "iterazionifatte"]; ob=document.getElementById("x300"); ob.innerHTML="Dati di inizializzazione del terzo test "+ vedoris(p1,inforis1); ris3=odedopri(terzafder,x_0,y_0,x_1,costa,peso,tol,maxiter); ob=document.getElementById("x301"); ob.innerHTML="Avendo fatto al massimo "+maxiter+" iterazioni "+ " in questo caso riesce al primo colpo "+ vedoris(ris3,inforis2); if(!ris3[1])ris4=odedopri(terzafder,ris3[2],ris3[3],x_1, costa,peso,tol,maxiter) else ris4=ris3; ob=document.getElementById("x302"); ob.innerHTML="Risultato finale "+vedoris(ris4,inforis2); // // ============================================================ // // Quarta prova che commento in dettaglio... // Il risultato finale dovrebbe essere ??? // // Innanzi tutto calcola i dati che servono per impostare i calcoli // p1=ini_fder4(); inforis1=[0,"t-iniziale","xterra","yterra","xluna","yluna","vxterra", "vyterra","vxluna","vyluna","t-finale","costa G*Msole", "Mterrasusole","Mlunasusole","Distadalsole","Vtangen", "pesoxt","pesoyt","pesoxl","pesoyl","pesovxt", "pesovyt","pesovxl","pesovyl","tollera","Maxpassi"]; inforis2=[0,"Bene se true ma male se false","t-finale", "xtfin","ytfin","xlfin","ylfin","vxtfin","vytfin","vxlfin", "vylfin","Iterfatte"]; ob=document.getElementById("x400"); ob.innerHTML="Dati di inizializzazione del quarto test "+ vedoris(p1,inforis1); // // Poi fa un numero massimo di 20 iterazioni. // ris3=odedopri(quartafder,x_0,y_0,x_1,costa,peso,tol,20); ob=document.getElementById("x401"); ob.innerHTML="Avendo fatto preliminarmente "+20+" iterazioni "+ vedoris(ris3,inforis2); // // Se le sole 20 iterazioni non sono bastate ad ottenere true // fa realmente un altro gruppo di maxiter iterazioni, ma se sono // bastate le prime 20 ad arrivare al tempo x_1, non ripete il calcolo... // Notare che ris3[2] e' il tempo a cui e' riuscito ad arrivare con 20 // iterazioni e ris3[3] e' il vettore delle condizioni iniziali. // if(!ris3[1])ris4=odedopri(quartafder,ris3[2],ris3[3],x_1, costa,peso,tol,maxiter) else ris4=ris3; ob=document.getElementById("x402"); ob.innerHTML="Avendo fatto preliminarmente "+maxiter+ " iterazioni non ce la fa "+vedoris(ris4,inforis2); // // Se non e' riuscito a raggiungere il tempo x_1 rifa' un altro // blocco di maxiter iterazioni sperando che almeno questo terzo tentativo // gli consenta di arrivare alla conclusione ossia al tempo x_1 desiderato. // if(!ris4[1])ris5=odedopri(quartafder,ris4[2],ris4[3],x_1, costa,peso,tol,maxiter) else ris5=ris4; // ob=document.getElementById("x403"); ob.innerHTML="Risultato finale ottenuto ridando "+ "i due dati, tempo e valori ottenuti "+ "e facendo un altro blocco di iterazioni con successo "+ vedoris(ris5,inforis2); // } // //_________________________________________________________________________ // // Una funzione di servizio per trasformare vettori // e vettori di vettori in stringhe. // function vedoris(rr,inforr){ // Visualizza il vettore rr var j,pn=[],npn,paragrafo,zz,ta=typeof([0]), ss=" {\u2192 ",nn=rr.length; if(arguments.length>2)paragrafo=arguments[2]; else {paragrafo=[],inforr[0]=0}; npn=paragrafo.length; for(j=0;npn>j;j++)pn[j]=paragrafo[j]; pn[npn]=0; for(j=1;nn>j;j++){ pn[npn]=j; if(typeof(rr[j])==ta) { zz=vedoris(rr[j],inforr,pn); ss+="[\u00ab"+pn+"\u00bb: "+zz+"] "; } else { zz=rr[j]; inforr[0]=Math.min(inforr[0]+1,inforr.length-1); ss+="[\u00ab"+pn+":"+inforr[inforr[0]]+ "\u00bb: "+zz+"] ";} }; return ss+" \u2190} "; } // // // Funzione generica per salvare nella libreria voluta ossia lib, // nell'elemento qui ( se e' stato creato oppure in coda ) // i dati contenuti in rr che puo' essere un dato o una lista o // una lista contenente altre liste... // function salvinlib(lib,qui,rr){ var nn=Math.min(qui,lib.length); lib[nn]=[]; salvorr(lib[nn],rr); } // function salvorr(lista,rr){ var j,mm,nn=rr.length,ta=typeof([0]); // Salva tutto quello che sta in rr for(j=0;nn>j;j++){ mm=lista.length; if(typeof(rr[j])==ta){ lista[mm]=[]; salvorr(lista[mm],rr[j]);} else lista[mm]=rr[j]; } } // /*AMEN LIBRERIA DORMAND PRINCE - MAGGIO 2014