MathStat2
Relazioni funzionali

G. D'Agostini

Primo sguardo ai dati

I seguenti dati sperimentali sono relativi al semplice esperimento della caduta di un grave legato ad un volano tramite un filo inestensibile. I valori della seguente lista rappresentano i tempi in secondi registrati in corrispondenza di traguardi elettronici (fotocellule) equidistanziati `esattamente' 5 centimetri (valore nominale). Essi sono stati effettivamente raccolti in laboratorio il 16 ottobre 2000 e rappresentano la prima di 10 serie di misure, ovvero essi non sono stati assolutamente selezionati. Le altre serie verranno utilizzati nel seguito. Per adesso vediamo cosa si può imparare da questa prima serie.

Tempi = {0. , 1.016434, 1.556013, 1.976810, 2.329293, 2.644171, 
          2.927985, 3.194652, 3.437351, 3.665107, 3.883512, 4.091030, 
          4.289556, 4.482196, 4.663909, 4.842332, 5.011102, 5.181138, 
          5.344091, 5.507233};

Grafichiamo i valori in successione.

ListPlot[Tempi, Prolog -> PointSize[0.018]];

Che legge fisica seguono i punti? Per metterlo in evidenza, trasformiamo il numero d'ordine in coordinata spaziale (che indicheremo genericamente con s) e creiamoci una lista bidimensionale tempo-spazio.

DeltaX = 5. ;
ST = 
        Table[{Tempi[[i]], (i - 1)*DeltaX}, {i, 20}];

Ecco come si presenta l'andamento:

GraficoST = 
        ListPlot[ST, Prolog -> PointSize[0.02], 
          PlotRange -> {{0, 6}, {0, 100}}, 
          AxesLabel -> {"t (s)", "s (cm)"}];

L'andamento sembra parabolico. Diciamo `sembra' perché onestamente tale affermazione è più legata ai nostri pregiudizi (leggi `conoscenze di fisica generale') che alle nostre capacita` di giudicare ad occhio che i punti giacciono su una parabola che su una curva qualsiasi altra famiglia. Per convincerci, facciamo le seguenti operazioni sui dati.

Linearizzazioni e fit lineari

Come è ben boto, è abbastanza facile convincersi se dei punti giacciono approssimativamente su una  retta, in quanto l'occhio umano (sarebbe meglio dire il cervello) è sensibile agli oggetti allineati o non allineati. Ben più difficile è giudicare se i punti giacciono su una curva di un'altra famiglia. Si usa quindi la tecnica di linearizzare degli andamenti trasformando opportunamente le variabili basandosi sui pregiudizi teorici legati alla fenomenologia che si sta studiando (o andando per tentativi in caso di indagini esplorative).

Ritenendo l'andamento del tipo s = (1/2)a (t^2), si vede che, graficando lo spazio percorso in funzione di (t^2), si otterrà un andamento lineare, di coefficiente angolare pari ad a/2.

ST2 = 
        Table[{Tempi[[i]]^2, (i - 1)*DeltaX}, {i, 20}];

GraficoST2 = 
        ListPlot[ST2, Prolog -> PointSize[0.02], 
          PlotRange -> {{0, 35}, {0, 100}}, 
          AxesLabel -> {"t^2 (s^2)", "s (cm)"}, 
          DisplayFunction -> Identity];

RettaST2 = Fit[ST2, {1, t2}, t2]

PlotRettaST2 = 
        Plot[RettaST2, {t2, 0, 36}, 
          DisplayFunction -> Identity];

Show[GraficoST2, PlotRettaST2, 
        DisplayFunction -> $DisplayFunction];

Come si vede, l'accordo fra punti e retta non è completamente soddisfacente. In questi casi è buona norma cercare di raffinare il modello teorico per ottenere un migliore accordo, anche se, come vedremo, bisogna diffidare dei `fit perfetti'. Raffinare il modello non significa sostituire una retta ad una parabola o altro polinomio di ordine superiore, ma cercare l'effetto fisico che può giustificare tale discrepanza. Bisogna anche cominciare ad imparare che discrepanze residue rimmarranno sempre e che grandi discrepanze possono verificarsi per caso anche se non ci sono altri effetti fisici, a parte `rumore' nelle misure, da tener conto.

Nel costro caso, l'effetto è banale.  Il modello usato non tiene conto di una possibile, seppur piccola velocità iniziale. Nel tale caso, l'equazione del moto diventa s = (1/2)a (t^2) + v0 t  ("+ s0" nel caso generale, ma non nel nostro esperimento, nel quale abbiamo fatto partire il cronometro elettronico in corrispondenza di s=0). Un altro modo di linearizzare l'andamendo, pur in presenza di una velocità iniziale diversa da zero è di  graficare s/t in funzione di t, chiaramente per t≠0. Il coefficiente angolare darà a/2 e l'intercetta sarà pari a v0. La grandezza s/t  ha chiaramente il significato di velocità media fino al tempo t. Indichiamo quindi la lista con il nome VmT.

VmT = 
        Table[{Tempi[[i]], (i - 1)*DeltaX/Tempi[[i]]}, {i, 2, 
            20}];

GraficoVmT = 
        ListPlot[VmT, Prolog -> PointSize[0.015], 
          PlotRange -> {{0, 6}, {0, 18}}, 
          AxesLabel -> {"t(s)", "s/t (cm/s)"}, 
          DisplayFunction -> Identity];

RettaVmT = Fit[VmT, {1, t}, t]

Da cui l'accelerazione in cm/s:

2  Part[RettaVmT, 2]/t

Plottiamo la retta del fit e savrapponiamola ai punti sperimentali.

PlotRettaVmT = 
        Plot[RettaVmT, {t, 0, 6}, DisplayFunction -> Identity];

Show[GraficoVmT, PlotRettaVmT, 
        DisplayFunction -> $DisplayFunction];

L'accordo fra punti sperimentali e andamento retta sembra ottimo e quindi i valori di accelerazione e velocità ricavati in questo modo quelli più attendibili. Per curiosità vediamo cosa si sarebbe ottenuto utilizzando soltanto il primo e l'ultimo punto.

RettaVmT2Punti = 
      Fit[Drop[VmT, {2, Length[VmT] - 1}], {1, t}, t]

RettaVmT

L'accordo è più che soddisfacente. Esso insegna che i punti più importanti ai fini della valutazione della retta sono il primo e l'ultimo, in quanto fra di essi c'è il maggior braccio di leva. Se invece si prendono coppie di punti consecutivi si osservano maggiori fluttuazioni rispetto ai valori determinati dalla globalità dei dati:

For[i = 1, i < Length[VmT],  i++,  
      r = Fit[Take[VmT, {i, i + 1}], {1, t}, t]; 
      Print["Retta fra i punti ", i, " e ", i + 1, "  -> ", 
        r]]

Fit Quadratico

Abbiamo visto come i parametri di un andamento lineare possono essere presi come misure delle grandezze fisiche ad essi associati, a condizione che tale andamento modellizzi realisticamente il fenomeno fisico di interesse. Ripetiemo che l'importanza delle linearizzazioni è legata alla nostra facilità di rendersi conto di allineamienti e, qualora si facessero conti senza il sussidio di computer, alla facilità di tracciare linee dritte con in righello. È comunque possibile adattare ai punti sperimentali curve più complicate, mediante opportuni algoritmi che, per dirla alla buona, "minimizzino la distanza" fra curva e punti, in analogia a quando si traccia  con il righello una retta tra due punti. Vediamo come adattare un andamento del tipo c1 t + c2 (t^(2 )) ai dati sperimentali, dove c1 e c2 sono simboli generici per i coefficienti.

QuadST = Fit[ST, {t, t^2}, t]

Ne segue un'accelerazione di

2  Part[QuadST, 2]/t^2

Facciamo il plot e sovrapponiamolo ai pubti sperimentali.

PlotQuadST = 
        Plot[QuadST, {t, 0, 6}, DisplayFunction -> Identity];

Show[GraficoST, PlotQuadST, 
        DisplayFunction -> $DisplayFunction];

Come si vede, anche questo andamendo è in ottimo accordo con i dati sperimentali. Confrontiamo i valori di accelerazione e di velocità iniziale con quanto ottenuto precedentemente dall'andamendo che ci sembrava ugualmente corretto.

Print["Retta   Vm(t) :  v0 = ", 
      Part[RettaVmT, 1], " cm/s,  a = ", 
      2  Part[RettaVmT, 2]/t, " cm/s^2"]
Print["Parabola s(t) :  v0 = ", 
      Part[QuadST, 1]/t, " cm/s,  a = ", 
      2  Part[QuadST, 2]/t^2, " cm/s^2"]

Test dei due metodi su dati simulati con andamento ideale

Sebbene i valori siano molto simili, essi non sono esattamente uguali. Per capirne il motivo simuliamo dei dati che seguano perfettamente l'andamento s = (1/2) a (t^2)+ v0 t, con a = 5.5 cm/(s^2) e v0 = 2.2 cm/s.

f[t_] := 5.5/2*t^2 + 2.2*t
STsim = Table[{0.29*i, f[0.29*i]}, {i, 0, 19}];
VmTsim = 
        Table[{STsim[[i, 1]], STsim[[i, 2]]/STsim[[i, 1]]}, {i, 
            2, 20}];

RettaVmTsim = Fit[VmTsim, {1, t}, t]

QuadSTsim = Fit[STsim, {t, t^2}, t]

Come si vede, in questo caso ideale si ottengono risultati esattamente identici. Quindi, la differenza fra i due risultati nel caso reale è da imputare alla dispersione dei dati sperimentali intorno all'andamento ideale.

Cosa fa Fit? `Distanze' minimizzabili.

È ora il momento di capire meglio cosa fa esattamente Fit:

?Fit

Tale funzione usa il metodo dei minimi quadrati, ottenuto (vedi Help di Mathematica) minimizzando la somma dei quadrati fra f(x_i) e y_i, dove x_i è l'i-mo valore  grandezza graficata sulle ascisse, y_i  l'analogo per le ordinate e f() è la funzione teorica che lega x e y ideali. Quindi, quella che noi chiamavamo genericamente "distanza', è in realtà la somme di quadrati di distanze calcolati lungo la y. Per ora non tentiamo di dare alcuna giustificazione teorica di questo algoritmo. Anzi, lo possiamo considerare come uno dei tanti possibili, avente il vantaggio formale che è più minimizzare tale oggetto matematico piuttosto che altre "distanza" più complicate (somma dei moduli degli scarti lungo la y, somma dei quadrati o dei moduli degli scarti lungo la x, somma dei quadrati o dei moduli delle distanze sul piano x-y, etc.). Come si può capire intuitivamente, a parte il caso ideale di punti perfettamente allineati,  ad ogni distanza corrispondono diverse stime dei parametri  (basta fare esempi estremi, ad esempio si pensi ad un fit lineare di quattro punti posti ai vertici di un quadrato).

Vediamo di riprodurre i risultati di Fit per l'andamento lineare della velocità media in funzione del tempo (list VmT). Per ora, utilizzando la potenza di Mathematica risolviamo il problema a `forza bruta' (seppur analiticamente!). Vedremo a tempo debito come arrivare a formule più maneggevoli per l'applicazione del metodo dei minimi quadrati.

f =. 
f[x_] := m x  +  c
DistMQ = 
        Sum[(f[VmT[[i, 1]]] - VmT[[i, 2]])^2, {i, 1, 
            Length[VmT]}];
DistMQ = Simplify[DistMQ]

Solve[{D[DistMQ, m] == 0, D[DistMQ, c] == 0}, {m, c}]   
MinDistMQ = FindMinimum[DistMQ, {m, 2, 3}, {c, 2, 3}]

Come si vede, si ottengono esattamente gli estessi risultati di Fit, a prova che i due algoritmi sono equivalenti. La ricerca del minimo è stata eseguita sia analiticamente, che numericamente tramite Findminimum. Infatti, volendo controllare, per curiosità, cosa si ottiene minimizzando altre 'distanze', non sarà più possibile usare calcoli analitici.

Vediamo cosa succede scegliendo la somma delle distanze lungo la y.

DistModY = 
        Sum[Abs[f[VmT[[i, 1]]] - VmT[[i, 2]]], {i, 1, 
            Length[VmT]}];
MinDistModY = FindMinimum[DistModY, {m, 2, 3}, {c, 2, 3}]

Ci ricaviamo prima la formula della distanza fra punto e retta nel caso bidimensionale (l'uso intermedio di tre dimensioni è per poter usare Cross)

\(P1r = {x, m\ x + c, 0}\ ; \ \ \ \ \ \ \ \ \ \  (*\ 
      generico\ punto\ sulla\ retta\ *) 

    vr = {c/m, c, 0}\ ; \ \ \ \ \ \ \ \ \  (*\ 
      vettore\ parallelo\ alla\ retta\ *) \ \ 

    vrmod = Sqrt[
        vr . vr]; \ \  (*\ \(...\)\(\ \)\(suo\)\(\ \)\(modulo\)\ *) 

    P0 = {x0, y0, 0}\ ; \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      punto\ del\ piano\ x\ y\ distante\ `dist'\ *) 

    dv = P1r - P0\ \ ; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      vettore\ da\ P0\ a\ P1r\ \ *) 

    vz = Cross[vr, dv]\ ; \ \ \ \ \ \ \  (*\ 
      prodotto\ vettoriale\ fra\ vr\ e\ dv\  -> \ lungo\ z*) 

    crossmod = 
      Sqrt[vz . 
          vz]\ ; \ \ \  (*\ \(...\)\(\ \)\(suo\)\(\ \)\(modulo\)\ \ \ *) \)
dist = crossmod/vrmod

Riscriviamo e semplifichiamo a mano (Mathematica se la cava molto bene con problemi complicati, ma a volte si perde in banalità ...).

fundist[x0_, y0_, m_, c_] := 
      Abs[c +  m x0 -  y0]/(Sqrt[m^2  + 1])

Minimizziamo ora distanza e distanza al quadrato (si provi a togliere il punto e virgola dalla prossima istruzione per avere un'idea dell'orrenda espressione...). Mettiemo i risultati in MinDistDist e MiDistDist2, che mostreremo successivamente insieme a quelli dei due metodi già incontrati.

DistDist = 
        Sum[fundist[VmT[[i, 1]], VmT[[i, 2]], m, c], {i, 1, 
            Length[VmT]}];
MinDistDist = 
        FindMinimum[DistDist, {m, 2, 3}, {c, 2, 3}];

DistDist2 = 
        Sum[fundist[VmT[[i, 1]], VmT[[i, 2]], m, c]^2, {i, 1, 
            Length[VmT]}];
MinDistDist2 = 
        FindMinimum[DistDist2, {m, 2, 3}, {c, 2, 3}];

Ecco, insieme i parametri delle quattro rette.

Print["Minimi Quadrati standard:     ", 
      MinDistMQ[[2]]]
Print["Somma distanze lungo y:       ", 
      MinDistModY[[2]]]
Print["Somma distanze nel piano:     ", 
      MinDistDist[[2]]]
Print["Somma quadr. dist. nel piano: ", 
      MinDistDist2[[2]]]

I valori sembrano praticamente simili, anche se diciamo per il momento `sembrano' in quanto dobbiamo ancora valutare l'incertezza da attribuire alla misura e quindi capire se le differenze fra i metodi sono `significative' o no. Altro punto da chiarire è se sia possibile ottenere una delle procedure incontrate in base ad argomentazioni generali sull'analisi dei dati, oppure se dobbiamo ritenerle più o meno ugualmente valide, insieme alle altre infinite minimizazioni possibili. Per ora seguitiamo a lavoare con i minimi quadrati, che sarà con opportune modifiche il metodo meglio fondato teoricamente, sotto ipotesi ben definite che andranno vagliate di volta in volta.

Concludiamo facendo osservare come, invece, la minimizzazione della 'distanza' (nel senso specificato) punti-retta, e che sembrerebbe il migliore dal punto di vista puramente geometrico, presenta un problema di interpretazione. Qualcuno avrà già notato che nella formula della distanza compare a denominatore 1+(m^2). Dal punto di vista fisico m ha le dimensioni [Y]/[X], nel nostro caso cm/s^2. Quindi il suo quadrato non può essere sommato ad 1. Il problema nasce dal fatto che è la distanza nel piano ad essere mal definita nel caso ascissa e ordinata abbiano diverse dimensioni  (se non si possono sommare mele e patate, tanto meno si possono sommare i loro quadrati). Quindi tutte le operazioni fatte precedentemente si intendono riferite ai valori numerici delle grandezze e non alle grandezze.

Analisi dei residui

Diamo ora uno sgardo ravvicinato all'accordo `soddisfacente' fra fit e dati sperimentali. Per fare questo, grafichiamo i residui, ovvero le differenze fra la componente `y' dei dati e la funzione fittata. Se le deviazioni dei dati sperimentali e l'andamento ideale è dovuto a piccoli effetti aleatori, ci aspettiamo residui positivi e negativi distribuiti `a caso' intorno allo zero.

ResVmT = 
        Table[{VmT[[i, 1]], 
            VmT[[i, 2]] - RettaVmT /. t -> VmT[[i, 1]]}, {i, 1, 
            19}];

GraficoResVmT = 
        ListPlot[ResVmT, PlotRange -> {{0, 6}, {-0.1, 0.1}}, 
          Prolog -> PointSize[0.02], 
          AxesLabel -> {"t(s)", "Res Vm (cm/s)"}];

ResST = 
        Table[{ST[[i, 1]], 
            ST[[i, 2]] - QuadST /. t -> ST[[i, 1]]}, {i, 1, 
            20}];

GraficoResST = 
        ListPlot[ResST, PlotRange -> {{0, 6}, {-0.1, 0.1}}, 
          Prolog -> PointSize[0.02], 
          AxesLabel -> {"t(s)", "Res s (cm)"}];

Entrame le distribuzioni dei residui mostrano qualcosa di sospetto, anche se le leggi del caso non proibiscono che delle cose a caso mostrino apparenti strutture (si pensi alle costellazioni del firmamento). A questo punto, se vogliamo capire qualcosa di più dobbiamo passare ad analizzare le altre serie di misure, in quanto "Einmal ist Keinmal".

Analisi di tutte le serie di misure.

Riportiamo tutte le misure in una sola grande lista bidimensionale, che richiamiamo Tempi.

Tempi = {Tempi, {0, 1.085254, 1.630131, 2.053967, 2.408044, 2.724131, 
            3.008754, 3.276158, 3.519370, 3.747497, 3.966447, 4.174128, 
            4.373025, 4.565777, 4.747730, 4.926470, 5.095287, 5.265103, 
            5.427422, 5.588757}, {0, 1.068144, 1.61106, 2.035277, 2.389096, 
            2.704712, 2.989123, 3.256609, 3.499438, 3.727759, 3.946807, 
            4.154152, 4.353041, 4.545679, 4.727502, 4.906120, 5.074691, 
            5.244506, 5.406430, 5.567767}, {0. , 1.078358, 1.621914, 
            2.046783, 2.400443, 2.716117, 3.000911, 3.268559, 3.511345, 
            3.739757, 3.958853, 4.166161, 4.365218, 4.557905, 4.739650, 
            4.918255, 5.086945, 5.256917, 5.418815, 5.580020}, {0, 1.051640, 
            1.593325, 2.016730, 2.369708, 2.684988, 2.969174, 3.236645, 
            3.479116, 3.707301, 3.92648, 4.133710, 4.332379, 4.525146, 
            4.706941, 4.885343, 5.054073, 5.223947, 5.385940, 5.547225}, {0, 
            1.087837, 1.632937, 2.057339, 2.411299, 2.726711, 3.011155, 
            3.279085, 3.521723, 3.750176, 3.969745, 4.177082, 4.375825, 
            4.568780, 4.750554, 4.928937, 5.097858, 5.268059, 5.430183, 
            5.591583}, {0, 1.113437, 1.658983, 2.084863, 2.43970, 2.754896, 
            3.039868, 3.308008, 3.550610, 3.779192, 3.998620, 4.20609, 
            4.404865, 4.598049, 4.779779, 4.958168, 5.127261, 5.297282, 
            5.459234, 5.620765}, {0, 1.085539, 1.628906, 2.054344, 2.408466, 
            2.723703, 3.008533, 3.276344, 3.518958, 3.747268, 3.966672, 
            4.173993, 4.372523, 4.565659, 4.747412, 4.925650, 5.094479, 
            5.264495, 5.426269, 5.587634}, {0, 1.085826, 1.628915, 2.05598, 
            2.409042, 2.724125, 3.008811, 3.276969, 3.519487, 3.747961, 
            3.967412, 4.174745, 4.373302, 4.566422, 4.748140, 4.926387, 
            5.094913, 5.265138, 5.426985, 5.588394}, {0, 1.044745, 1.585056, 
            2.008122, 2.361843, 2.675765, 2.960571, 3.228319, 3.470171, 
            3.698223, 3.917823, 4.124650, 4.323320, 4.516258, 4.697902, 
            4.876048, 5.044906, 5.214840, 5.376391, 5.537727}};
Length[Tempi]

Avendo una discreta massa di dati, cerchiamo di operare in modo compatto.

Carichiamo nelle liste bidimensionali i dati di interesse (spazio e velocità media in funzione del tempo).

ST = 
        Table[Table[{Tempi[[j, i]], (i - 1)*DeltaX}, {i, 20}], {j, 
            10}];
VmT = 
        Table[Table[{Tempi[[j, i]], (i - 1)*
                DeltaX/Tempi[[j, i]]}, {i, 2, 20}], {j, 
            10}];

Facciamo i rispettivi fit.

\(For[j = 1, j <= Length[ST], \(j++\), 
      RetteVmT = Table[Fit[VmT[\([i]\)], {1, t}, t], {i, Length[VmT]}]; 

      QuadST = Table[
          Fit[ST[\([i]\)], {t, t^2}, t], {i, Length[ST]}]]\)

Nellla lista v0 poniamo le 10 coppie di velocità iniziali delle 10 serie, ottenute con i due metodi (in ordine velocità media in funzione del tempo e spazio in funzione del tempo).

v0 = 
        Table[{Part[RetteVmT[[j]], 1], Part[QuadST[[j]], 1]/t}, {j, 
            Length[ST]}];
Print["Risultati fit per v0 (lineare e quadratici) \n serie     VmT   \
     ST "]
For[j = 1, j <= 10, j++, 
      Print["   ", j, "     ", v0[[j, 1]], "    ", 
        v0[[j, 2]]]]

Idem per le accelerazioni (lista a).

a = 
        Table[{Part[RetteVmT[[j]], 2]/t*2, 
            Part[QuadST[[j]], 2]/t^2*2}, {j, Length[ST]}];
Print["Risultati fit per a (lineare e quadratici) \n serie     VmT    \
    ST "]
For[j = 1, j <= 10, j++, 
      Print["   ", j, "     ", a[[j, 1]], "    ", 
        a[[j, 2]]]]

Facciamo un po` di statistica (descrittiva) sui risultati, caricando l'opportuno add-on.

<< Statistics`DescriptiveStatistics`

Differenza fra v0 calcolate con i due metodi (parabola-retta)

Deltav0 = 
      Table[v0[[j, 2]] - v0[[j, 1]], {j, Length[ST]}]

Il secondo metodo dà sempre valori maggiori. Facciamo medie e deviazione standard delle differenza:

Mean[Deltav0]
StandardDeviation[Deltav0]

Idem per le accelerazioni.

Deltaa = 
      Table[a[[j, 2]] - a[[j, 1]], {j, Length[ST]}]

Mean[Deltaa]
StandardDeviation[Deltaa]

Come si vede, i due metodi danno valori sistematicamente diversi. Più esattamente, il fit quadratico tende a dare velocità iniziali maggiori del 3% circa, mentre dà accelezioni minori dello 0.5%. Le piccole deviazioni standard indicano che la dispersione fra le differenze è molto piccola.

Per quanto riguarda invece la dispersione dei valori di v0 e a sulle 10 serie di misure abbiamo, a parità di metodo.

Print["v0:  Fit lineare VmT"]
\(Print["Media: ", 
      media = Mean[\(Transpose[v0]\)[\([1]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              v0]\)[\([1]\)]], "; Coeff. var.: ", 

      devst/media]\)

Print["v0: Fit quadratico ST"]
\(Print["Media: ", 
      media = Mean[\(Transpose[v0]\)[\([2]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              v0]\)[\([2]\)]], "; Coeff. var.: ", 

      devst/media]\)

Print["a: Fit lineare VmT"]
\(Print["Media: ", 
      media = Mean[\(Transpose[a]\)[\([1]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              a]\)[\([1]\)]], "; Coeff. var.: ", 

      devst/media]\)

Print["a: Fit Quadratico ST"]
\(Print["Media: ", 
      media = Mean[\(Transpose[a]\)[\([2]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              a]\)[\([2]\)]], "; Coeff. var.: ", 

      devst/media]\)

Dal confronto dei risultati si evince quanto segue. 
   -  La dispersione dei valori di v0 è molto grande. Si noti come il parametro significativo di dispersione non è la dispersione assoluta dei dati (quantificata dalla deviazione standard), ma la dispersione relativa al valore centrale (quantificata dalla deviazione standard diviso il valore medio, rapporto che prende il nome di coefficiente di variazione). I risultati su v0 mostrano una dispersione relativa dell 11% circa, giustificato che v0  vari da una serie di misure all'altra. Esso dipende, infatti, dall'esatta posizione da cui il peso viene rilasciato. La differenza di risultato fra i due metodi applicati agli stessi dati è molro minore della dispersione intrinseca di v0. 
   - La dispersione dei valori di accelerazione è, invece, molto piccola, consistentemente con il fatto che l'accelerazione con cui il peso cade è una caratteristica del sistema (peso,volano,forza di gravità) e può considerrsi costante, a parte piccoli effeti aleatori. Si noti anche come la differenza fra i due metodi è dello stesso ordine di grandezza della dispersione dei valori. Diventa quindi cruciale decidewre quale dei due metodi è il più accurato.

Analisi globale dei residui

Riprendiamo il dubbio lasciato in sospeso a proposito della distribuzione dei residui della prima serie di misure. Occupiamoci soltanto dei residui del plot delle velocità medie in funzione del tempo.

Per mettere tutti i residui nello stesso plot (con MultipleListPlot), applichiamo delle traslazioni multpliple di 0.1 ai valori dei residui.

ResVmT = 
        Table[Table[{VmT[[j, i, 
                  1]], (VmT[[j, i, 2]] - RetteVmT[[j]] /. 
                    t -> VmT[[j, i, 1]]) + (j - 1)*0.1}, {i, 1, 
              19}], {j, 10}];

<< Graphics`MultipleListPlot`

MultipleListPlot[ResVmT[[1]], ResVmT[[2]], ResVmT[[3]], 
        ResVmT[[4]], ResVmT[[5]], ResVmT[[6]], ResVmT[[7]], 
        ResVmT[[8]], ResVmT[[9]], ResVmT[[10]], 
        PlotRange -> {{0, 6}, {-0.1, 1. }}, Prolog -> PointSize[0.02], 
        AxesLabel -> {"t(s)", "Res Vm (cm/s)"}];

ResST = 
        Table[Table[{ST[[j, i, 
                  1]], (ST[[j, i, 2]] - QuadST[[j]] /. 
                    t -> ST[[j, i, 1]]) + (j - 1)*0.1}, {i, 1, 
              20}], {j, Length[ST]}];

MultipleListPlot[ResST[[1]], ResST[[2]], ResST[[3]], 
        ResST[[4]], ResST[[5]], ResST[[6]], ResST[[7]], 
        ResST[[8]], ResST[[9]], ResST[[10]], 
        PlotRange -> {{0, 6}, {-0.1, 1. }}, Prolog -> PointSize[0.02], 
        AxesLabel -> {"t(s)", "Res s (cm)"}];

Effettivamente i residui sembrano avere tutti la stessa forma, ad indicare che, anche se di poco, il modello usato non descrive bene i dati. Si può allora pensare di introdurre un effetto di una forza dipendente dalla velocità, causato dalla resistenza dell'aria. Innanzitutto si noti come quest'effetto va nella direzione indicata dai dati sperimentali. Infatti, l'attrito smorza diminuisce l'accelerazione, fino ad annullarla quando la forza di attrito è pari alla somma delle altre forze esterne (forza di fravità e tensione del filo). Di conseguenza, al crescere del tempo, sia lo spazio percorso che la velocità tenderanno ad essere inferiori a quanto calcolato in assenza di attrito. Quindi le rispettive curve tenderanno a "curvarsi verso il basso". Prendiamo il caso del fit della velocità media in funzione del tempo, più perché (ma analogo ragionamento vale anche per lo spazio percorso in funzione del tempo). Siccome ogni piccola deviazione alla linarità può vista, in prima approssimazione, come un contributo quadratico alla curva, l'effetto di attrito trasformerà la retta in una parabola avente derivata seconda negativa. Cercando di adattare una retta fra punti aventi andamento parabolico, i punti esterni presenteranno residui negativi, quelli centrali positivi, come si osserva nei grafici. Una volta che convinti che l'effetto sia questo e, ammesso di essere interessati a calcolare l'accelerazione iniziale (per esempio perché essa è legata in modo abbastanza diretto al momento di inerzia del volano) si possono analizzare i dati per i quali l'effetto di resistenza dell'aria è trascurabile, ad esempio i primi cinque o dieci punti (servono un po` di punti per avere un discreto braccio di leva, come abbiamo visto precedentemente).

Fit limitati ai primi dieci traguardi

Fit (in VmT si prendono i primi 9 valori, in quanto  si ricorda che il primo traguardo non era usato per quasta analisi).

For[j = 1, j <= Length[ST], j++, 
      RetteVmTInizio = 
        Table[Fit[Take[VmT[[i]], 9], {1, t}, t], {i, Length[VmT]}]; 
      QuadSTInizio = 
        Table[Fit[Take[ST[[i]], 10], {t, t^2}, t], {i, 
            Length[ST]}]]

Coppie di risultati per v0

v0Inizio = 
      Table[{Part[RetteVmTInizio[[j]], 1], 
          Part[QuadSTInizio[[j]], 1]/t}, {j, Length[ST]}]

... e per a.

aInizio = 
      Table[{Part[RetteVmTInizio[[j]], 2]/t*2, 
          Part[QuadSTInizio[[j]], 2]/t^2*2}, {j, Length[ST]}]

Print["v0: Fit lineare VmT (primi punti)"]
\(Print["Media: ", 
      media = Mean[\(Transpose[v0Inizio]\)[\([1]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              v0Inizio]\)[\([1]\)]], "; Coeff. var.: ", 

      devst/media]\)

Print["v0: Fit quadratico ST (primi punti)"]
\(Print["Media: ", 
      media = Mean[\(Transpose[v0Inizio]\)[\([2]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              v0Inizio]\)[\([2]\)]], "; Coeff. var.: ", 

      devst/media]\)

Print["a: Fit lineare VmT (primi punti)"]
\(Print["Media: ", 
      media = Mean[\(Transpose[aInizio]\)[\([1]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              aInizio]\)[\([1]\)]], "; Coeff. var.: ", 

      devst/media]\)

Print["a: Fit Quadratico ST (primi punti)"]
\(Print["Media: ", 
      media = Mean[\(Transpose[aInizio]\)[\([2]\)]], 
"; Dev. St.:  ", 
      devst = StandardDeviation[\(Transpose[
              aInizio]\)[\([2]\)]], "; Coeff. var.: ", 

      devst/media]\)

Osserviamo ora i residui (mettiamo anche i punti non inclusi nel fit per avere un quadro generale).

ResVmT = 
        Table[Table[{VmT[[j, i, 
                  1]], (VmT[[j, i, 2]] - RetteVmTInizio[[j]] /. 
                    t -> VmT[[j, i, 1]]) + (j - 1)*0.1}, {i, 1, 
              19}], {j, 10}];

MultipleListPlot[ResVmT[[1]], ResVmT[[2]], ResVmT[[3]], 
        ResVmT[[4]], ResVmT[[5]], ResVmT[[6]], ResVmT[[7]], 
        ResVmT[[8]], ResVmT[[9]], ResVmT[[10]], 
        PlotRange -> {{0, 6}, {-0.1, 1. }}, Prolog -> PointSize[0.02], 
        AxesLabel -> {"t(s)", "Res. Vm (cm/s)"}, 
        PlotLabel -> "    [Fit sui primi 9 punti]"];

ResST = 
        Table[Table[{ST[[j, i, 
                  1]], (ST[[j, i, 2]] - QuadSTInizio[[j]] /. 
                    t -> ST[[j, i, 1]]) + (j - 1)*0.1}, {i, 1, 
              20}], {j, Length[ST]}];

MultipleListPlot[ResST[[1]], ResST[[2]], ResST[[3]], 
        ResST[[4]], ResST[[5]], ResST[[6]], ResST[[7]], 
        ResST[[8]], ResST[[9]], ResST[[10]], 
        PlotRange -> {{0, 6}, {-0.1, 1. }}, Prolog -> PointSize[0.02], 
        AxesLabel -> {"t(s)", "Res s (cm)"}, 
        PlotLabel -> "    [Fit sui primi 10 punti]"];

Come si vede, in entrambi i casi la distribuzione dei residui è molto `buona' dove le funzioni sono ipotizzate descrivere bene i dati, mentre per tempi grand, i grafici mettono in maggiore evidenza un ulteriore effetto.

Fit "perfetti" e ragioni per diffidarne

Un modo a bruta forza per ottenere un fit "perfetto" è quello di adattare fra i punti una curva abbastanza complicata. Come è noto, fra due punti nel piano passa una e una sola retta, per tre una e una sola parabola, eccetera. Quindi potremmo essere tentati di fittare n punti non un polinomio di grado n-1. I residui verranno esattamente zero e il fit sarà perfetto. Proviamo a trasformare il nostro fit lineare sui dati della prima serie di misure in un fit con polinomio di 18 grado (con Mathematica ce lo possiamo permettere).

Fit18 = Fit[VmT[[1]], Table[t^i, {i, 0, 18}], t]

Questo è il risultato del fit "perfetto" nella regione nella regione compresa fra il tempo minimo e il tempo massimo usati nel fit.

PlotFit18 = 
        Plot[Fit18, {t, VmT[[1, 1, 1]], VmT[[1, 19, 1]]}, 
          DisplayFunction -> Identity];
Show[GraficoVmT, PlotFit18, 
      DisplayFunction -> $DisplayFunction]

Sebbene la curva tocchi esattamente i punti, essa ha un comportamento alquanto strano fra i punti misurati.  Se poi ci si allontana dai punti misurati il problema `esplode', nel senso che, in base alla nostra esperienza, non crediamo assolutamente che l'estrapolazione corrisponderebbe alla previsione della grandezza di interesse.

PlotFit18 = 
        Plot[Fit18, {t, 0, 6}, DisplayFunction -> Identity];
Show[GraficoVmT, PlotFit18, PlotRange -> {{0, 6}, {0, 30}}, 
        DisplayFunction -> $DisplayFunction];

Quindi, non basta modificare la formula per avere un accordo perfetto con quanto si è osservato. Bisogna trovare una funzione motivata fisicamente che descriva ragionevolmente quanto osservato e preveda con un certo grado di fiducia quanto non è ancora stato osservato.

Raffinamento del modello

Equazioni del moto

s

Caduta nel vuoto

Ricaviamoci ora la legge più generale che tenga conto della resistenza del mezzo modellizzata come forza in direzione del moto e verso opposto e di modulo proporzionle alla velocità (approssimazione valida a basse velocità). Siccome dobbiamo risolvre delle equazioni differenziali con Mathematica, cominciamo dal semplice di caduta in vuoto. Per utilizare dei simboli intuitivi, resettiamo la variabile a usata precedentemente.

\(\(\(\(ClearAll["Global`a"]\)\(\ \)\)\(;\)\(\ \)\(Remove[
      a]\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)\)

Ricordiamo che l'equazione del moto usata è soluzione della seguente equazione differenziale, ottenuta assumendo che, durante la caduta, il corpo sia soggetto ad una forza costante, data da due contributi conrapposti: la forza di gravità mg e la tensione del filo rivolta verso l'alto T (ignoriamo i dettagli di cosa trattiene il filo verso l'alto, ovvero volano e carrucola).

sol1 = DSolve[m s''[t] == m g - T, s[t], t];
s = s[t] /. sol1[[1, 1]]

ove le costanti C[1] e C[2] rappresentano posizione e velocità iniziali. Velocità e accelerazione sono data da

v = D[s, t]

a = D[v, t]

Attrito dell'aria

Passiamo ora al modello più realistico (anche se, per non entrare nei dettagli del moto del volano, tratteremo T come una costante). Sia b il coefficiente di proporzionalità fra forza e velocità.

ClearAll["Global`a"]; ClearAll["Global`s"];  
    ClearAll["Global`v"]; ClearAll["Global`v0"];
Remove[a, s, v, 
      v0]

sol2 = 
        DSolve[m v'[t] == m g - T - b v[t], v[t], t];
v = v[t] /. sol2[[1, 1]]

Eliminiamo la costante C[1], esprimendola in termini della velocità iniziale v0 e delle altre grandezza in gioco.

sol = Solve[(v /. t -> 0) == v0, C[1]]

v = v /. sol[[1, 1]];
v = Simplify[v]

Imponiamo ora che l'accelerazione iniziale sia pari a a0 e sostituiamo T:

solT = Solve[Limit[D[v, t], t -> 0] == a0, T]

v = v /. solT[[1, 1]]

Controlliamo il risultato  generale facendo dei limiti a situazioni più  semplici.

Limit[v, b -> 0]

v /. a0 -> g /. v0 -> 0

Questo è il solo caso (caduta di corpo senza reazioni vincolari e velocità iniziale nulla)  del quale si trova soluzione in normali libri di testo. Ricontrolliamone il limite per b->0.

Limit[%, b -> 0]

Calcoliamo anche l'espressione dell'accelerazione nel caso generale:

a = D[v, t]

La velocità limite è data da

vlim = Simplify[v /. (b/m t) ->  Infinity]

Si noti come in questo caso - banale! - ci sono problemi a fare il limite con Mathematica, in quanto non sembra ci sia un modo per convincere il programma a considerare positivo il rapport b/m che compare all'esponente.

Dall'espressione della velocità in funzione del tempo otteniamo l'equazione del moto s(t).

sol3 = DSolve[s'[t] == v, s[t], t];
s = s[t] /. sol3[[1, 1]]

Imponiamo la condizione s[t=0]=0 per eliminare C[1].

sol = Solve[Limit[s, t -> 0] == 0, C[1]]

s = s /. sol[[1, 1]]

Simplify[%]

Sostituiamo il valore numerico di m (in Kg).

s = % /. m  -> 0.600

Utilizziamo nel seguito una funzione ad assegnazione ritardata (l'alternativa ad usare la funzione Evaluate sarebbe stata di ricopiare l'espressione, ottenendo lo stesso risultato, ma precludendo la possibilità di eseguire automaticamente l'intero notebook in caso di modifiche della notazione o del valore numerico di m).

sAria[t_] := Evaluate[s]

Idem per la velocità velocità media e istantanea:

vmAria[t_] := sAria[t]/t
vAria[t_] := Evaluate[v /. m -> 0.600]

Rianalisi dei dati della posizione in funzione del tempo

Questa volta non possiamo usare la semplice Fit. Potremmo minimizzare i quadrati degli scarti usando funzioni elementari di Mathematica come fatto precedentemente, ma preferiamo usare una funzione degli add-on, che ci fa tutto automaticamente anche nei casi più complicati quale questo. Si tratta di NonlinearRegress di Statistics`NonLinearFit`.

<< Statistics`NonLinearFit`

FitSTAria = 
      Table[NonlinearRegress[ST[[j]], sAria[t], 
          t, {{a0, 5.5, 0, 10}, {v0, 2. , 1, 3}, {b, 0.001, 0, 0.01}}, 
          RegressionReport -> BestFitParameters], {j, 10}]

Come si vede, Si ottengono valori di a leggermente più alti di quelli ottenuti precedentemente e valori di v0 un po' più bassi (consistente con quanto deve succedere, si pensi al perché).

Check: T e (Subscript[F, a]).

Facciamo ora alcune verifiche: calcoliamo T e il valore  massimo della forza dovuta all'attrito dell'aria. Approfittiamo per vedere come Mathematica gestisce le unità di misura.

mgValues={mg  -> {g -> 980, 
            m -> 0.600}}

T /. 
          solT[[1, 1]] /. (mg /. 
            mgValues) /. (BestFitParameters /. 
          FitSTAria[[1]])

Aggiungiamo le unità di misura, facendo uso del package Miscellaneous`Units`.

<< Miscellaneous`Units`

%% Kilogram Centimeter/ Second^2

Sistema Internazionale.

SI[%]

Trasformiamo in Newton.

TFitST = Convert[%, Newton]

Questa forza fa confrontata con quella di gravità:

Fg = 
      Convert[(m g /. (mg /. 
                  mgValues)) Kilogram  Centimeter Second^-2, 
        Newton]

Per la forza di attrito abbiamo:

(Subscript[F, a]) = 
      SI[(b vAria[t] /. (BestFitParameters /. FitSTAria[[1]])) 
          Kilogram  Centimeter Second^-2]

Per t=0 e t=tmax abbiamo, rispettivamente

FAriaFitSTmin = Convert[(Subscript[F, a]) /. t -> 0, Newton]
FAriaFitSTmax = 
      Convert[(Subscript[F, a]) /. t -> ST[[1, 20, 1]], Newton]

Come si vede, con una analisi accurata si possono metter in risalto effetti abbastanza piccoli (si tratta, al massimo,  del 10% circa della forza netta iniziale).

Residui

shift = 0.1;
ResSTAria = 
        Table[Table[{ST[[j, i, 1]], 
              ST[[j, i, 
                    2]] - (sAria[
                      ST[[j, i, 1]]] /. (BestFitParameters /. 
                        FitSTAria[[j]])) + (j - 1)*shift}, {i, 1, 
              20}], {j, 10}];

MultipleListPlot[ResSTAria[[1]], ResSTAria[[2]], 
        ResSTAria[[3]], ResSTAria[[4]], ResSTAria[[5]], 
        ResSTAria[[6]], ResSTAria[[7]], ResSTAria[[8]], 
        ResSTAria[[9]], ResSTAria[[10]], 
        PlotRange -> {{0, 6}, {-0.1, 1. }}, Prolog -> PointSize[0.02], 
        AxesLabel -> {"t(s)", "Res s (cm)"}];

In media la situazione è molto migliore della precedente, anche se ci sono un paio di punti che risultano sistematicamente fuori. Potrebbe essere a fluttuazioni nel posizionamento delle fotocellule o delle letture dei tempi (si noti come la struttura dei residui si ripete: le fluttuazioni sono quindi sistematiche del sistema e non casuali e variabili da una lettura all'altra). Quello che è importante è che non si vedono andamenti che facciano pensare che sia la legge fisica ad essere inadeguata, come succedeva precedentemente.

sST = 
      sAria[t] /. (BestFitParameters /. FitSTAria[[1]])

PlotSTAria = 
        Plot[sST, {t, 0, 6}, DisplayFunction -> Identity];

Show[GraficoST, PlotSTAria, 
        DisplayFunction -> $DisplayFunction];

Terminiamo questa parte facendo una statistica generale sui residui, togliendo quelli relativi al tempo zero, i quali devono essere identicamente nulli avendo imposto s[t=0]=0 nell'equazione del moto. Ovviamente togliamo anche gli shift [(j-1)*0.1] aggiunti precedentemente per ragioni grafiche.

NewRes = 
        Table[Drop[ResSTAria[[j]], 1], 
{j, 
            Length[ResSTAria]}];
AllRes = 
        Flatten[Join[
            Table[(Transpose[
                      NewRes[[j]]][[2]] - (j - 1)*0.1), {j, 
                Length[ResSTAria]}]]];
Length[AllRes]
Min[AllRes]
Max[AllRes]
Mean[AllRes]
StandardDeviation[AllRes]

Come si vede, su 190 valori, gli scarti minimi e massimi sono dell'ordine del millimetro, con una deviazione standard di 4 decimi di millimetro.

Rianalisi dati velocità medie in funzione del tempo

Ripetiamo l'analisi anche per la velocità media in funzione del tempo.

FitVmTAria = 
      Table[NonlinearRegress[VmT[[j]], vmAria[t], 
          t, {{a0, 5.5, 0, 10}, {v0, 2. , 1, 3}, {b, 0.001, 0, 0.01}}, 
          RegressionReport -> BestFitParameters], {j, 10}]

Residui

Mostriamo ora l'insieme dei residui, per convircerci che effettivamente stiamo descrivendo bene il fenomeno.

ResVmTAria = 
        Table[Table[{VmT[[j, i, 
                  1]], (VmT[[j, i, 2]] - 
                      vmAria[VmT[[j, i, 1]]] /. (BestFitParameters /. 
                        FitVmTAria[[j]])) + (j - 1)*0.1}, {i, 1, 
              19}], {j, 10}];

MultipleListPlot[ResVmTAria[[1]], ResVmTAria[[2]], 
        ResVmTAria[[3]], ResVmTAria[[4]], ResVmTAria[[5]], 
        ResVmTAria[[6]], ResVmTAria[[7]], ResVmTAria[[8]], 
        ResVmTAria[[9]], ResVmTAria[[10]], 
        PlotRange -> {{0, 6}, {-0.1, 1. }}, Prolog -> PointSize[0.02], 
        AxesLabel -> {"t(s)", "Res Vm (cm/s)"}];

Valgono gli stessi commenti fatti a proposito dei residui della distribuzione dello spazio in funzione del tempo. Grafichiamo punti sperimentali e risultato del fit:

vMA = 
      vmAria[t] /. (BestFitParameters /. FitVmTAria[[1]])

PlotVmTAria = 
        Plot[vMA, {t, 0, 6}, DisplayFunction -> Identity];

Show[GraficoVmT, PlotVmTAria, 
        DisplayFunction -> $DisplayFunction];

Incertezze di misura - primissimo approccio

A questo punto avremmo desiderio  di giungere a delle conclusioni sul valore delle grandezze fisiche che volevamo misurare, velocità e accelerazione iniziali. Abbiamo visto come, ripetendo l'esperimento non riproduciamo mai gli stessi valori numerici. Inoltre, i valori cambiano sistematicamente, anche se di poco, con il metodo di analisi. A quale dobbiamo credere?

Riportiamo qui valore medio e deviazione standard di v0 e a ottenuti dalle 10 serie di misure con i due metodi, tenendo conto della resistenza dell'aria.

v0VmTAria = 
      Table[v0 /. (BestFitParameters /. FitVmTAria[[j]]), {j, 
          10}]

Print["Media ", Mean[v0VmTAria], "; StDev ", 
      StandardDeviation[v0VmTAria]]

v0STAria = 
      Table[v0 /. (BestFitParameters /. FitSTAria[[j]]), {j, 
          10}]

Print["Media ", Mean[v0STAria], "; StDev ", 
      StandardDeviation[v0STAria]]

aVmTAria = 
      Table[a0 /. (BestFitParameters /. FitVmTAria[[j]]), {j, 
          10}]

Print["Media ", Mean[aVmTAria], "; StDev ", 
      StandardDeviation[aVmTAria]]

aSTAria = 
      Table[a0 /. (BestFitParameters /. FitSTAria[[j]]), {j, 
          10}]

Print["Media ", Mean[aSTAria], "; StDev ", 
      StandardDeviation[aSTAria]]

Volendo arrotondare, possiamo dire di sentirci sicuri di dire che l'accelerazione iniziale vale circa 5.7 cm/s^2 e la velocità iniziale intorno a 1.5--2 cm/s. Per questa seconda diciamo "intorno" in quanto, come abbiamo detto, essa non è costante da misura a misura. Concentriamoci soltanto sull'accelerazione, la grandezza che sembra essere misurata più precisamente (nel senso che il linguaggio comune dà a questa parola). Possiamo dare il risultato in modo più preciso?

Considerando valori medi e dispersione dei risultati ottenuti con i due metodi, possiamo dire che il valore vero di a0 potrebbe trovarsi con molta probabilità nell'intervallo fra 5.67 e 5.74 cm/s^2 ottenuto considerando la minore delle due medie meno la relativa deviazione standard e la maggiore delle due medie più la relativa deiazione standard. Oppure si può considerare il range entro cui fluttuano i valori di accelerazione (ma chi ci dice che se rifacessiomo altre serie di misure il range non si allargherebbe?):

Min[Join[aVmTAria, aSTAria]]
Max[Join[aVmTAria, aSTAria]]

Queste affermazioni sono abbastanza vaghe, anche se esprimono una certa confidenza che il valore vero si trovi in quei range. Ma per trasformare questa affermazione in una affermazione quantitativa di uso scientifico bisogna anche dare il livello di probabilità. Molto probabile significa al 50%, al 90% o al 99.99%?  Altrimenti si potrebbe dare sempre un intervallo di valori molto ampio (da 0 a 980 cm/s^2, per esempio) dichiarandoci praticamente sicuri che il valore vero sia compreso all'interno di esso. Non che l'affermazione sia sbagliata. È semplicemente non informativa, nel senso che non abbiamo imparato niente dall'esperimento. Il giusto mezzo fra fornire intervalli eccessivamente conservativi, ma non informativi, e l'ideale della misura esatta è dare un intervallo associando ad esso un valore di probabilità che indichi quanto noi crediamo in tale affermazione.

Purtroppo per arrivare ad affermzioni probabilistiche sensate bisogna capire meglio cosa si intende per probabilità, modellizzare il processo di misura in modo probabilistico, capire il processo di inferenza  dai dati sperimenati alle grandezze fisiche e imparare, finalmente, a dare affermazioni probabilistiche sui valori delle grandezze fisiche. Quindi, a questo punto ci dobbiamo fermare e studiare la teoria della probabilità.

Volendo comunque riassumere le cause che, nell'esempio che stiamo studiando,  ci inducono ad essere incerti sul risultato della misura abbiamo:
- ripetendo più la stessa misura in condizioni apparentemente identiche non otteniamo due risultati uguali;
- i risultati dipendono dal metodo di analisi e dal modello teorico alla luce del quale i dati sono interpreatati.

Appendice: ripesare i punti sperimenatli

Anticipando concetti che vedremo nel seguito, ritorniamo al problema dei valori diversi ottenuti dal fit lineare e dal fit quadratico, e che abbiamo imputato alla disersione dei punti sperimentali rispetto all'andamento ideale. 
Il metodo dei minimi quadrati implementato in Mathematica tratta in modo democratico tutti i punti sperimenatali, nel senso che, quando prova ad adattare la curva fra di essi, li considera tutti con lo  stesso peso.  Proviamo,
per esercizio, a dare dei pesi diversi ai punti sperimentali della velocità media in funzione del tempo. Senza tentare di giustificare la ragione, diamo ad essi dei pesi proporzionali al quadrato del tempo, e vediamo cosa otteniamo (NonLinearRegress permette tale flessibilità, il ché fa capire che pesare i dati in qualche modo rappresenta un modo consueto di anlzzare i dati).

pesi = (Transpose[VmT[[1]]][[1]])^2

\(FitVmTAriaPesi = 
      Table[NonlinearRegress[VmT[\([j]\)], vmAria[t], 
          t, {{a0, 5.5, 0, 10}, {v0, 2. , 1, 3}, {b, 0.001, 0, 0.01}}, 
          Weights -> pesi, 
			RegressionReport -> 
            BestFitParameters], {j, 10}]\)

aVmTAriaPesi = 
        Table[a0 /. (BestFitParameters /. FitVmTAriaPesi[[j]]), {j, 
            10}];
Print["Media ", Mean[aVmTAriaPesi], "; StDev ", 
      StandardDeviation[aVmTAriaPesi]]

v0VmTAriaPesi = 
      Table[v0 /. (BestFitParameters /. FitVmTAriaPesi[[j]]), {j, 
          10}]
Print["Media ", Mean[v0VmTAriaPesi], "; StDev ", 
      StandardDeviation[v0VmTAriaPesi]]

Si vede come si ottengono valori di accelerazione e velocità iniziale molto prossimi a quelli ottenuti dai dati di spazio in funzione del tempo, anche se non proprio perfettemante identici. Ad esempio, le medie delle accelerazioni iniziali sono uguali entro 2 parti su 10000 e le medie delle velocità iniziali uguali entro 8 parti su 10000. Anche le deviazioni standard sono praticamente identiche. Sembra quindi che il ripesare i punti sperimentali in questo modo senz'altro poco intuitivo, per ora, ha funzionato. Per giustificarlo si richiedono ulteriori conoscenze di probabilità e statistica inferenziale.

Altre linearizzazioni notevoli

Ci sono due linearizzazioni particolarmente importanti e che si incontrano frequentemente nella sperimentazione. Come al solito, cominciamo con degli esempi.

fun1 = 10.  E^(-t/4);
\(		\(fun2 = 15. \ E^\((\(-t\)/3)\);\)\)
fun3 = 0.1 E^(t/3);
\(		\(fun4 = 20. *10^\((\(-t\)/4)\);\)\)
\(		\(fun5 = 30. *10^\((\(-t\)/3)\);\)\)
fun6 = 0.5*10^(t/5);

Plot[{fun1, fun2, fun3, fun4, fun5, fun6}, {t, 0, 25}, 
        PlotRange -> {0, 32}, AspectRatio -> 0.6, 
        PlotLabel -> "        Andamenti esponenziali in scala normale", 
        AxesLabel -> {"t", "f(t)"}, 
        PlotStyle -> {{}, {}, {}, Dashing[{0.02}], Dashing[{0.02}], 
            Dashing[{0.02}]}];

Mettiamo ora in evidenza la loro natura esponenzialemediante LogPlot di Graphics`Graphics.

<< Graphics`Graphics`

LogPlot[{fun1, fun2, fun3, fun4, fun5, fun6}, {t, 0, 25}, 
        PlotRange -> {{0, 25}, {0.001, 32}}, 
        PlotLabel -> "   ... con ordinata in scala logaritmica", 
        AxesLabel -> {"t", "f(t)"}, 
        PlotStyle -> {{}, {}, {}, Dashing[{0.02}], Dashing[{0.02}], 
            Dashing[{0.02}]}];

Ogni andamento del tipo y = a (b^(c x)) apparirà lineare se l'ordinata viene graficata su scala logaritmica,ovvero la distanza tra due punti non è proporzionale allla loro reale differenza (principio sul quale si basano le mappe), ma alla differenza fra i loro logaritmi. Tale proprietà si ricava prendendo i logaritmi di entrambi i membri di y = a (b^(c x)).

\(\(fun7 = 2\ Sqrt[x];\)

\)
\(\(fun8 = \ 0.5\ \ x^1.5;\)
\)
\(\(fun9 = 20/x^\((3/2)\);\)

\)
fun10 = 10/Sqrt[x];

Plot[{fun7, fun8, fun9, fun10}, {x, 0, 25}, PlotRange -> {0, 32}, 
        PlotLabel -> "   Leggi di potenza su scale normali", 
        AxesLabel -> {"x", "f(x)"}, 
        PlotStyle -> {{}, Dashing[{0.015}], Dashing[{0.03}], 
            Dashing[{0.01, 0.03}]}];

LogLogPlot[{fun7, fun8, fun9, fun10}, {x, 0.001, 25}, 
        PlotRange -> {{0.05, 25}, {0.001, 10000}}, 
        PlotLabel -> "    Leggi di potenza su scale logaritmiche", 
        AxesLabel -> {"x", "f(x)"}, 
        PlotStyle -> {{}, Dashing[{0.015}], Dashing[{0.03}], 
            Dashing[{0.01, 0.03}]}];

Come si vede, sempre per le proprietà dei logaritmi, ogni andamento del tipo y = a (x^b) appare lineare se entrambe le coordinate sono plottate su scala logaritmica.

Opportune carte millimetrate logaritmiche permettono di fare agevolmente questi tipi di grafici, senza dovere calcolare i logaritmi dei valori da riportare. Come abbiamo visto, le funzioni di Mathematica LogPlot e LogLogPlot permettono di plottare in tal modo delle funzioni. Per graficare punti sperimentali esistono LogListPlot e LogLogListPlot, equivalenti di ListPlot.