v_20080412
Attenzione: questo file
non funziona con IE8...
Funziona con Firefox,
con Opera, con Safari, con Chrome....
In rete

Un bel modo di risolvere le ODE....

Asta penzolante con smorzamento viscoso interno:

r′ = v
s = λ² − r • ( r + ω · v )
f = h + σ · s · r
v′ = μ · f

Il sistema di equazioni differenziali è quello riportato nel titolo ed è facile intuire il significato di scalari e vettori in gioco.

La lunghezza dell'asta è λ mentre la sua rigidità è data da σ e la sua viscosità interna è ω.
Maggiore è il valore di σ e più l'asta tende a comportarsi in modo rigido ossia come un corpo rigido mentre l'asta è tanto più dissipativa quanto più elevata è la sua viscosità interna.
Lo scalare μ rappresenta la mobilità dell'estremo pesante ossia l'inverso della sua massa.
Tutti gli scalari qui indicati con lettere greche sono positivi.
Si è usato il primo ( ′ = ′ ) per indicare la derivazione rispetto al tempo, l'operatore pallino ( • = • ) rappresenta il prodotto scalare tra vettori che vengono indicati in grassetto. L'operatore moltiplicazione non è mai sottinteso per consentire l'uso di simboli anche di due o più caratteri; uso il carattere puntinmezzo ( · = · )
Il vettore r rappresenta le coordinate dell'estremo pesante e non vincolato dell'asta mentre l'estremo vincolato è incernierato nell'origine. Il vettore v rappresenta la velocità dell'estremo pesante e il vettore f rappresenta la forza agente su tale punto.
Il vettore h, costante nel tempo, rappresenta la forza di gravità mentre 1/μ è la massa dell'estremo pesante non vincolato e dunque pendolante. Pertanto, se la forza di gravità è diretta in direzione opposta all'asse dell'ordinate e se assumiamo che il moto sia bidimensionale, otteniamo h = [ 0, − g / μ ] essendo g l'accelerazione di gravità ovvero 9.80665 m·s⁻² nel Sistema Internazionale.

Il modello dell'asta con estremi pesanti e dotata di rigidità non infinita e viscosità interna non nulla è molto importante perchè rappresenta l'elemento base per rappresentare in modo realistico strutture reticolari a parametri concentrati.

Nel caso bidimensionale il sistema di ODE consta di quattro equazioni differenziali ed in Javascript potrebbe essere codificato in questo modo:

function fa_derivate(p){
     var ss,fx,fy; // var.esterne: ql,om,sg,hy,mu.
     d1 = p[3]; 
     d2 = p[4];
     ss = ql - p[1]*(p[1] + om*p[3]) - p[2]*(p[2] + om*p[4]);
     fx = sg*ss*p[1];
     fy = hy + sg*ss*p[2];
     d3 = mu*fx;
     d4 = mu*fy;
     return [ 0, d1, d2, d3, d4, 1 ];
     }

ossia la funzione fa_derivate(...) riceve il vettore della posizione e velocità dell'estremo pesante e produce come risultato il vettore delle derivate prime.
In questa funzione si è rispettata la convenzione di non utilizzare il primo elemento del vettore risultato e di generare come ultima derivata, quella della variabile tempo che, ovviamente, vale sempre 1.

== Massimo passi
== Viscosità interna
== Rigidezza
== Massa
== Tempo finale



== x ,y 

== vx,vy

== Tempo

Note:

Questa pagina serve da esempio su come fare per inserire un disegno SVG in una pagina HTML.
La soluzione però non funziona con Internet Explorer 7 e neppure con IE8, ma con Firefox 2.0.0.13 ( e successive ), con Safari 3.1 e con Opera 9.27 ( anche se con una segnalazione di errore).
Questa pagina, sperimentata ad inizio 2011 funziona anche con Chrome,

Può risultare interessante vedere come è stato implementato in JavaScript l'algoritmo di J.R.Cash e A.H.Karp ossia un algoritmo di tipo Runge Kutta, di alta precisione numerica e universalmente consigliato per l'integrazione dei sistemi di equazioni differenziali ordinari ( ODE == Ordinary differential equations ).Nel caso si ritenesse opportuno si consulti: en.wikipedia.org/wiki/Numerical_ordinary_differential_equations

---

Runge Kutta di Cash e Karp

Si tratta di un algoritmo di elevata precisione ( quint'ordine) che fornisce anche una soluzione di minore precisione (quart'ordine) da usarsi per il controllo del passo di avanzamento temporale.
Riporto solo la versione a parametri interi che ho dedotto dalla versione tradizionale. Specificando solo parametri interi viene eliminata qualsiasi approssimazione nella definizione dell'algoritmo che non dipende più dal numero di cifre significative disponibili a livello software. L'algoritmo è dunque così definito:

k1 = h · f( y ) u2 = 16632 · k1 k2 = h · f( y + 11776 · u2 ) u3 = 6237 · ( k1 + 3 · k2 ) k3 = h · f( y + 11776 · u3 ) u4 = 24948 · ( k1 − 3 · k2 + 4 · k3 ) k4 = h · f( y + 11776 · u4 ) u5 = 1540 · ( −11 · k1 + 135 · k2 − 140 · k3 + 70 · k4 ) k5 = h · f( y + 11776 · u5 ) u6 = 3262 · k1 + 37800 · k2 + 4600 · k3 + 44275 · k4 + 6831 · k5 k6 = h · f( y + 8855 · u6 ) u7 = 9361 · k1 + 38500 · k3 + 20125 · k4 + 27648 · k6 yn = y + 10240 · u7 u8 = 39550 · k1 + 148600 · k3 + 94675 · k4 + 7479 · k5 + 96768 · k6 ynt = y + 2530 · u8

Ad ogni passo il tempo incrementa di   979292160 · h = 4096·27·5·7·11·23 · h   per cui il valore di h deve essere circa un miliardesimo dell'incremento che si desidera ottenere.
Il valore di y_n rappresenta l'effettivo risultato mentre il vettore y_nt deve essere usato solo a scopo di test; se si discosta eccessivamente ( ritenuto eccessivo in base ad opportuni criteri più o meno euristici ) da quello di y_n bisogna rieffettuare il calcolo riducendo opportunamente h.
Con f(y) si intende la funzione che, ricevuto in argomento il vettore delle variabili considerate, fornisce come risultato il vettore delle loro derivate prime. Si noti che il vettore delle variabili considerate deve includere anche la variabile considerata come indipendente, ossia ,generalmente, il tempo. Per la variabile indipendente la derivata prima è, per definizione, 1. Nella implementazione in Javascript ho imposto la convenzione che la variabile indipendente sia sempre l'ultima e quindi il vettore delle derivate prime contiene sempre 1 come valore del suo ultimo elemento.

Bibliografia:
J.R.Cash, A.H.Karp - ACM Transaction on Mathematical Software, Vol. 16, No. 3, September 1990, Pages 201-222.
Vedere il cap. 16.2 di Numerical Recipes

Contatti (G.Bottoni)