Problemi di Fisica con l'uso del PC

XI Lezione: elementi di programmazione ed applicazioni. Moto in 2 dimensioni.

Iterazioni

L'istruzione Do permette di iterare un set di operazioni. Esempio:

sum = 0;
Do[sum = sum + i; Print[sum], {i, 1, 5}]

Esempio: costruiamo una sequenza equivalente a NDSolve. 
Consideriamo il caso del moto in una dimensione di un corpo soggetto ad una forza elastica:

                (( (d^2) x)/(d (t^2)))=(F/m)    ;       F = -k x

Dati e condizioni iniziali:

Clear[f]
f[x_, v_, t_] = -2*x;
m =  .5;
t0 = 0. ;
t1 = 10. ;
x0 = 0. ;
v0 =  .5;

Troviamo la soluzione con NDSolve:

Clear[x]
sol = NDSolve[{x''[t] - f[x[t], x'[t], t]/m , x'[0], x[0]}==
{ 0, v0, x0}, x, {t,  t0,  t1}]

Flatten[sol]

\(\(Clear[xes];\)	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)
xes[t_] := x[t] /.  Flatten[sol]

xes[t]

Plot[xes[t],  {t,  t0,  t1}, 
      PlotRange -> All]

Il programma equivalente si ottiene iterando la sequenza:

(( d x(t))/(d t))=v(t)    ;(( d v(t))/(d t))=((F(x, v, t))/m)
cioe'
x(t + dt) =  x(t)  + v(t) dt   ;  v(t + dt)= v(t) +((F(x, v, t))/m) dt   

con dt fissato e sufficientemente piccolo:

npoint = 200;
dt = ((t1 - t0)/npoint);
Clear[xx, vv];
xx[0] = x0;
vv[0] = v0;

Do[
xx[k + 1] = 
        xx[k] + vv[k]*dt; 
vv[k + 1] = 
        vv[k] + f[xx[k], vv[k], t]/m*dt, 
{k, 0, 
        npoint}]

xvec = Table[{t0 + k*dt, xx[k]}, {k, 0, npoint}];

ListPlot[xvec, PlotRange -> All]

Show[
      Plot[xes[t] ,  {t,  t0,  t1}, PlotRange -> All, 
        PlotStyle -> RGBColor[1, 0, 0]], 
      ListPlot[xvec, PlotRange -> All]]

Notare come al crescere di npoint, cioe' al diminuire di dt la soluzione approssimi sempre meglio quella esatta. 
Abbiamo ottenuto una soluzione xvec per punti. Interpoliamola:

Clear[xap];
xap[100] = Interpolation[xvec]

xap[100][3.67]

Plot[xap[100][t] ,  {t,  t0,  t1}, 
      PlotRange -> All]

npoint va scelto come compromesso tra velocita' e precisione. Perfezioniamo il programma: partiamo da un npoint piccolo e lo raddoppiamo ogni volta finche' l' "errore medio" definito dalla media del  valore assoluto della differenza di due soluzioni successive

(((∫_to^t1) | (Subscript[xap, 2 npoint])(t) - (Subscript[xap, npoint])(t)|d t)/(t1 - t0))

non e' uguale almeno a un valore prefissato ϵ .

xap[200] = Interpolation[xvec]

Plot[{xap[100][t], xap[200][t]} ,  {t,  t0,  t1}, 
      PlotRange -> All]

Plot[Abs[xap[100][t] - xap[200][t]],  {t,  t0,  t1}, 
      PlotRange -> All]

Inseriamo per questo il Do precedente in un nuovo Do nel quale si calcola l'errore e che viene interrotto dal controllo If appena l'errore prefissato viene raggiunto:

ϵ =  .01;  
nini = 10;   
npoint = 
        (nini/2);
Clear[
        xapold];
xapold[t_] = 
        0;
\(Do[
npoint = 
            2*npoint; 
dt = \(t1 - t0\)\/npoint; 
Clear[xx, vv]; 
xx[0] = 
            x0; 
vv[0] = 
            v0; 

Do[
xx[k + 1] = xx[k] + vv[k]*dt; 
vv[k + 1] = 
              vv[k] + f[xx[k], vv[k], t]\/m*dt, 
{k, 0, 
              npoint}]; 

xvec = 
            Table[{t0 + k*dt, xx[k]}, {k, 0, 
                npoint}]; 
Clear[
            xapnew]; 
xapnew = 
            Interpolation[xvec]; 
err = 
            1\/\(t1 - t0\)*
              NIntegrate[
                Abs[xapnew[tt] - xapold[tt]], {tt, t0, 
                  t1}]; 
Print[npoint, "    ", 
            err]; 
Clear[
            xapold]; 
xapold[t_] = 
            xapnew[t]; 
If[err ≤ ϵ, 
            Break[]], 
{iter, 1, 20}];\)

Clear[
        xap];
xap[t_] = xapold[t];

Show[
      Plot[xes[t] ,  {t,  t0,  t1}, PlotRange -> All, 
        PlotStyle -> RGBColor[1, 0, 0]], 
      Plot[xap[t] ,  {t,  t0,  t1}, PlotRange -> All]]

Esercizio: aumentare gradualmente il valore della costante elastica: valutare come cambia
l'errore e spiegare.

Esercizio: rifare l'esempio aggiungendo una forza di attrito .

Pendolo.

Studiamo il moto di un corpo appeso a un filo lungo l . La forza del filo e' nulla se il filo non e' teso; elastica, con costante elevata, se e' teso: se (Overscript[((Subscript[x, 1])),  ->]) = {x1 , y1 }sono le coordinate del corpo appeso all'estremo libero del filo ( l'altro estremo e' fisso nell'origine del sistema di coordinate) la forza del filo vale:

(Overscript[((Subscript[F, el])),  ->])  = - k ( |(Overscript[((Subscript[x, 1])),  ->]) | - l ) (((Overscript[((Subscript[x, 1])),  ->]) )/(|(Overscript[((Subscript[x, 1])),  ->]) |))       con    k > 0       se | (Overscript[((Subscript[x, 1])),  ->])  | > l

(Overscript[((Subscript[F, el])),  ->])  = 0      altrimenti.

Clear[fel];
fel[x1_, y1_, vx1_, vy1_, t_] = 
        If[(Sqrt[(((x1))^2) + (((y1))^2)]) > 
              l, -100000. *((Sqrt[(((x1))^2) + (((y1))^2)]) - l), 
            0]*{x1, y1}/(Sqrt[(((x1))^2) + (((y1))^2)]);

Sul corpo agisce la forza (Overscript[((Subscript[F, el])),  ->])   piu' la forza di gravita'  m (Overscript[g,  ->])  :

(Subscript[((Overscript[F,  ->]) ), 1]) = (Overscript[((Subscript[F, el])),  ->])  +  m (Overscript[g,  ->])

Clear[f1];
f1[x1_, y1_, vx1_, vy1_, t_] = {0, -m1*g} + 
          fel[x1, y1, vx1, vy1, t];

Dati e condizioni iniziali:

m1 = 1. ;
g = 9.8;
l = 10;

    
t0 = 0. ;
t1 = 30. ;

    
{x10, 
          y10} = { .365*l*Sin[(π/4)], -l*
            Cos[(π/4)]};
{vx10, vy10} = {0, 0};

Soluzione : qui ho chiamato

vec = { m1 x1''(t) - (Subscript[f, x]) ,  m1 y1''(t) - (Subscript[f, y]) ,  x1(0) , y1(0) ,  x1'(0) ,  y1'(0)  }

vec0 = {0 , 0 , x10 , y10 , vx10 , vy10 }

Clear[eq, x1, y1, t];
eq = 
      D[{m1*x1[t], m1*y1[t]}, {t, 2}] - 
        f1[x1[t], y1[t], x1'[t], y1'[t], 
          t]
ini = {x1[0], y1[0], x1'[0], y1'[0]}
vec = Join[eq, ini]
vec0 = {0, 0, x10, y10, vx10, vy10}

sol  =  
        NDSolve[vec == vec0,  {x1, y1},  {t,  t0,  t1}, 
          MaxStepSize ->  .001, StartingStepSize ->  .001, 
          MaxSteps -> 1000000];
\(\(Clear[xx1, 
        yy1];\)	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)
xx1[t_] := 
        x1[t] /.  Flatten[sol][[1]];
yy1[t_] := 
        y1[t] /.  Flatten[sol][[2]];

    
Plot[xx1[t], {t, t0, t1}, 
      PlotRange -> All]
Plot[yy1[t], {t, t0, t1}, PlotRange -> All]

Notare come sia differente la soluzione nei casi in cui si parte col filo teso e non . Per fare un'animazione della soluzione utilizzeremo alcuni oggetti grafici  ( vedere Graphics su Help ):

Graphics[{RGBColor[1, 0, 0], Disk[{2, 3},  .2]}]

Show[%, AspectRatio -> 1, Axes -> True, 
      PlotRange -> {{-10, 10}, {10, -10}}]

Graphics[Line[{{0, 0}, {4, 5}}]]

Show[%, AspectRatio -> Automatic, Axes -> True, 
      PlotRange -> {{-10, 10}, {10, -10}}]

Costruiamo la sequenza delle immagini utilizzando la funzione Do:

niter = 30;
dt = ((t1 - t0)/niter);
t = t0;
Do[
t = t + dt; 
punto1 = 
          Graphics[{RGBColor[1, 0, 0], 
              Disk[{xx1[t], yy1[t]},  .2]}]; 
curva = 
          ParametricPlot[{xx1[tt], yy1[tt]}, {tt, t0, t}, 
            DisplayFunction -> Identity, 
            PlotStyle -> RGBColor[1, 0, 0]]; 
line = 
          Graphics[Line[{{0, 0}, {xx1[t], yy1[t]}}]]; 
Show[
          punto1, curva, line, Axes -> True, 
          PlotRange -> {{-10, 10}, {0, -13}}, 
          AspectRatio -> Automatic], 
{i, 1, 
          niter}];

Moto di due corpi soggetti a forze generiche.

Le equazioni del moto in 2 dimensioni di 2 corpi soggetti ad una forza reciproca (Subscript[F, 12])
piu' una forza esterna su ciascuno , (Subscript[F, 1]) ed (Subscript[F, 2]) , sono:

                (( (d^2) (Subscript[x, 1]))/(d (t^2)))=(((Subscript[F, 1]) + (Subscript[F, 12]))/((Subscript[m, 1])))    ;      (( (d^2) (Subscript[x, 2]))/(d (t^2)))=(((Subscript[F, 2])  -  (Subscript[F, 12]))/((Subscript[m, 2])))

Risolviamole    ed  animiamo la soluzione . Qui c'e' il caso di due corpi soggetti ad una attrazione gravitazionale reciproca che si muovono nel campo gravitazionale di un terzo corpo ( p. es. un sistema pianeta-satellite nel campo di una stella fissa ).

Clear[f1, f2, f12, f21, 
        t]
m1 = 
        1;
m2 =  .2;
mm = 5;

r1 = (Sqrt[(((x1))^2) + (((y1))^2)]);
r2 = (Sqrt[(((x2))^2) + (((y2))^2)]);
f1[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, 
            t_] = (-((m1*mm)/(r1^2)))*{x1, y1}/r1;  
f2[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, 
              t_] = 
(-((m2*
                        mm)/(r2^2)))*{x2, y2}/r2;  


r = (Sqrt[(((x2 - x1))^2) + (((y2 - y1))^2)]);
f12[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, 
            t_] = (((m1*m2)/(r^2)))*{x2 - x1, y2 - y1}/r;  
f21[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, 
              t_] = -f12[x1, y1, x2, y2, vx1, vy1, vx2, vy2, t];  

t0 = 
        0. ;
t1 = 3000. ;

{x10, y10} = {0, 
          30};
{x20, y20} = {0, 
          39};
{vx10, vy10} = { .38, 
          0};
{vx20, vy20} = { .05, 
          0};

Clear[eq, x1, y1, x2, y2, t];
eq = 
        D[{m1*x1[t], m1*y1[t], m2*x2[t], m2*y2[t]}, {t, 
              2}] - 
Flatten[{f12[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t] + 
                f1[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t], 
              f21[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t] + 
                f2[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], 
                  t]}];
ini = {x1[0], y1[0], x2[0], 
        y2[0], x1'[0], y1'[0], x2'[0], y2'[
          0]}
vec = Join[eq, ini]
vec0 = {0, 0, 0, 0, x10, y10, x20, y20, vx10, vy10, vx20, 
        vy20}

sol  =  
        NDSolve[vec == vec0,  {x1, y1, x2, y2},  {t,  t0,  t1}, 
          MaxStepSize ->  .1, StartingStepSize ->  .1, 
          MaxSteps -> 100000];
\(\(Clear[xx1, yy1, xx2, 
        yy2];\)	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)
xx1[t_] := 
        x1[t] /.  Flatten[sol][[1]];
yy1[t_] := 
        y1[t] /.  Flatten[sol][[2]];
xx2[t_] := 
        x2[t] /.  Flatten[sol][[3]];
yy2[t_] := 
        y2[t] /.  Flatten[sol][[4]];

    
Plot[{xx1[t], xx2[t]}, {t, t0, t1}, 
      PlotRange -> All]
Plot[{yy1[t], yy2[t]}, {t, t0, t1}, PlotRange -> All]

Per animare la soluzione servono numerose immagini , troppe per essere tenute nella memoria del notebook.
Generiamo percio' la sequenza mediante le istruzioni che comandano il cursore ( vedere Help ) :

nb=SelectedNotebook[];

(...........
SelectionMove[nb,Previous,EvaluationCell];
SelectionEvaluate[nb]
)

Il risultato equivale a una sequenza ottenuta con Do, pero' tenendo in memoria una sola immagine ad ogni istante dell'animazione.

nb = 
          SelectedNotebook[];
niter = 
        1000;
ddd = 50;
{{xmin, xmax}, {ymin, ymax}} = {{-ddd, ddd}, {-ddd, 
              ddd}};

dt = 
        N[((t1 - t0)/niter)];
t = 
        N[t0];
Clear[run];

punto0 = 
          Graphics[{RGBColor[1, 1, 0], PointSize[ .09], Point[{0, 0}]}];

run := (
t = t + dt; 
punto1 = 
            Graphics[{RGBColor[1, 1, 1], PointSize[ .06], 
                Point[{xx1[t], yy1[t]}]}]; 
punto2 = 
            Graphics[{RGBColor[1, 0, 0], PointSize[ .03], 
                Point[{xx2[t], yy2[t]}]}]; 
Show[punto0, 
            punto1, punto2, Axes -> True, 
            PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
            AspectRatio -> Automatic, 
            Background -> RGBColor[0, 0, 0]]; 
If[
            t > t1 - 2*dt, Abort[]]; 
SelectionMove[nb, 
            Previous, EvaluationCell]; 
SelectionEvaluate[
            nb];
);

run

Esercizio: provare con diversi tipi di forze e diversi dati iniziali .

Urto anelastico.

Il programma e' analogo al precedente. Abbiamo ora solo due corpi sferici di raggio (Subscript[r, 1]) ed (Subscript[r, 2]): il primo si conficca nel secondo perdendo energia a causa di una forza di attrito per esempio proporzionale alla velocita' relativa  del primo corpo rispetto al secondo:

(Overscript[F,  ->])  = - k ( (Overscript[((Subscript[v, 1])),  ->])  - (Overscript[((Subscript[v, 2])),  ->])  )      con    k > 0       se | (Overscript[((Subscript[x, 1])),  ->]) -(Overscript[((Subscript[x, 2])),  ->])  | < (Subscript[r, 1])+ (Subscript[r, 2])

(Overscript[F,  ->])  = 0      altrimenti.

Tutte le altre forze sono nulle.

Clear[f1, f2, f12, f21, t, r, 
        vrel]
m1 = 
        1. ;
m2 = 
        7. ;
r1 = 1. ;
r2 = 5. ;

f1[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, t_] = {0, 0};  
f2[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, t_] = {0, 
              0};  

f12[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, t_] = 
          If[(Sqrt[(((x2 - x1))^2) + (((y2 - y1))^2)]) ≤ (r1 + 
                    r2), - .2, 0]*{vx1 - vx2, vy1 - vy2};  
f21[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, 
            t_] = -f12[x1, y1, x2, y2, vx1, vy1, vx2, vy2, t];

t0 = 
        0. ;
t1 = 300. ;

{x10, y10} = {0, 
          0};
{x20, y20} = {50. , 
          0};
{vx10, vy10} = {2. , 
          0};
{vx20, vy20} = {0, 0};

Clear[eq, x1, y1, x2, y2, t];
eq = 
        D[{m1*x1[t], m1*y1[t], m2*x2[t], m2*y2[t]}, {t, 
              2}] - 
Flatten[{f12[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t] + 
                f1[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t], 
              f21[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t] + 
                f2[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], 
                  t]}];
ini = {x1[0], y1[0], x2[0], 
          y2[0], x1'[0], y1'[0], x2'[0], y2'[
            0]};
vec = Join[eq, ini];
vec0 = {0, 0, 0, 0, x10, y10, x20, y20, vx10, vy10, vx20, 
          vy20};
sol  =  
        NDSolve[vec == vec0,  {x1, y1, x2, y2},  {t,  t0,  t1}, 
          MaxStepSize ->  .1, StartingStepSize ->  .1, 
          MaxSteps -> 100000];
\(\(Clear[xx1, yy1, xx2, 
        yy2];\)	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)
xx1[t_] := 
        x1[t] /.  Flatten[sol][[1]];
yy1[t_] := 
        y1[t] /.  Flatten[sol][[2]];
xx2[t_] := 
        x2[t] /.  Flatten[sol][[3]];
yy2[t_] := 
        y2[t] /.  Flatten[sol][[4]];

    
Plot[{xx1[t], xx2[t]}, {t, t0, t1}, 
      PlotRange -> All]
Plot[{yy1[t], yy2[t]}, {t, t0, t1}, PlotRange -> All]

nb = 
          SelectedNotebook[];
niter = 200;

dt = 
        N[((t1 - t0)/niter)];
t = 
        N[t0];
Clear[run];
  
run := (
t = 
            t + dt; 
punto1 = 
            Graphics[{RGBColor[1, 0, 0], 
                Disk[{xx1[t], yy1[t]}, r1]}]; 
punto2 = 
            Graphics[{RGBColor[0, 0, 1], 
                Disk[{xx2[t], yy2[t]}, r2]}]; 
Show[punto2, 
            punto1, Axes -> True, 
            PlotRange -> {{-5, 140}, {-20, 20}}, 
            AspectRatio -> Automatic]; 
If[
            t > t1 - 2*dt, Abort[]]; 
SelectionMove[nb, 
            Previous, EvaluationCell]; 
SelectionEvaluate[
            nb];
);

run

Esercizio: dimostrare graficamente che la quantita' di moto del sistema si conserva.
Dimostrare graficamente che l'energia cinetica del sistema non si conserva.

Urto elastico.

Abbiamo due corpi di raggio r1 ed r2 : la forza tra essi e' repulsiva di modulo molto grande e diretta lungo la congiungente i loro centri (Overscript[((Subscript[x, 1])),  ->])  e (Overscript[((Subscript[x, 2])),  ->])  non appena si toccano ; nulla altrimenti. Detta r la distanza tra i centri possiamo rappresentarla

(Overscript[F,  ->])  = - k ( (Overscript[((Subscript[x, 2])),  ->])  - (Overscript[((Subscript[x, 1])),  ->])  )      con    k > 0       se | (Overscript[((Subscript[x, 1])),  ->]) -(Overscript[((Subscript[x, 2])),  ->])  | < (Subscript[r, 1])+ (Subscript[r, 2])

(Overscript[F,  ->])  = 0      altrimenti.

Clear[f1, f2, f12, f21, t, r, xx, rr, ll, vvx, 
        vvy, attr]
m1 = 
        1. ;
m2 = 
        1. ;
r1 = 2. ;
r2 = 2. ;

f1[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, 
            t_] = 
{0, 0};  
f2[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, t_] = {0, 0};  
r = (Sqrt[(((x2 - x1))^2) + (((y2 - y1))^2)]);
f12[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, t_] = 
          If[r ≤ (r1 + r2), -(10^4), 0]*{x2 - x1, 
              y2 - y1};  
f21[x1_, y1_, x2_, y2_, vx1_, vy1_, vx2_, vy2_, 
            t_] = -f12[x1, y1, x2, y2, vx1, vy1, vx2, vy2, t];

t0 = 
        0. ;
t1 = 2. ;

{x10, y10} = {6, 
          9};
{x20, y20} = {18, 
          7};
{vx10, vy10} = {20.3, 
          0.9};
{vx20, 
          vy20} = {0. ,  .0};
\(If[\@\(\((x20 - x10)\)\^2 + \((y20 - y10)\)\^2\) ≤ \((r1 + r2)\), Print["ERRORE"]];\)


Clear[eq, x1, y1, x2, 
        y2, t];
eq = 
        D[{m1*x1[t], m1*y1[t], m2*x2[t], m2*y2[t]}, {t, 
              2}] - 
Flatten[{f12[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t] + 
                f1[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t], 
              f21[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], t] + 
                f2[x1[t], y1[t], x2[t], 
                  y2[t], x1'[t], y1'[t], x2'[t], y2'[t], 
                  t]}];
ini = {x1[0], y1[0], 
          x2[0], y2[0], x1'[0], y1'[0], x2'[0], y2'[
            0]};
vec = 
        Join[eq, ini];
vec0 = {0, 0, 0, 0, x10, 
          y10, x20, y20, vx10, vy10, vx20, 
          vy20};
sol  =  
        NDSolve[vec == vec0,  {x1, y1, x2, y2},  {t,  t0,  t1}, 
          MaxStepSize ->  .1, StartingStepSize ->  .1, 
          MaxSteps -> 100000];
Clear[xx1, yy1, xx2, yy2];	                      
xx1[t_] := 
        x1[t] /.  Flatten[sol][[1]];
yy1[
          t_] := y1[
            t] /.  Flatten[
              sol][[2]];
xx2[t_] := 
        x2[t] /.  Flatten[sol][[3]];
yy2[t_] := y2[t] /.  Flatten[sol][[4]];

Plot[{xx1[t], 
        xx2[t]}, {t, t0, t1}]
Plot[{yy1[t], 
        yy2[t]}, {t, t0, t1}]

nb = 
          SelectedNotebook[];
niter = 
        80;
dt = 
        N[((t1 - t0)/niter)];
t = 
        N[t0];
Clear[
        run];
run := (
t = 
            t + dt; 
punto1 = 
            Graphics[{RGBColor[0, 0, 1], 
                Disk[{xx1[t], yy1[t]}, r1]}]; 
punto2 = 
            Graphics[{RGBColor[1, 0, 0], 
                Disk[{xx2[t], yy2[t]}, r2]}]; 
Show[punto1, 
            punto2, PlotRange -> {{0, 40}, {-10, 30}}, 
            AspectRatio -> Automatic, 
            Background -> RGBColor[0, 1, 0]]; 
If[
            t > t1 - 2*dt, Abort[]]; 
SelectionMove[nb, 
            Previous, EvaluationCell]; 
SelectionEvaluate[
            nb];
);

run

Esercizio: dimostrare graficamente che la quantita' di moto e l'energia cinetica del sistema si conservano.