Problemi di Fisica con l'uso del PC

VI Lezione

Forza elastica

F = - k x

Energia cinetica 	Ec = (1/2)m (v^2)

Energia potenziale 	Ev = (1/2)k (x^2)

Energia totale 		Et = Ec + Ev = cost

Oscillazioni armoniche

- k x[t] = -m x''[t]

oppure il sistema
x'[t]  =  v[t]
v'[t]  =  -(ω^2) x[t]
(ω^2)  =  k/m

Clear["Global`*"];
DSolve[{x'[t]  ==  
            v[t], v'[t] == -ω ω x[t]}, {x[t], 
          v[t]}, t] // FullSimplify

C[1] = A cosϕ  C[2] = -ω A sinϕ
x[t] = A cosϕ cosωt - A sinϕ sinωt = A cos(ωt + ϕ)
ϕ = ϕ' - π/2
x[t] = A cos(ωt + ϕ' - π/2) = A sin(ωt + ϕ)
ω : pulsazione
T = 2π/ω : periodo
ν = 1/T = ω/2π : frequenza

Clear["Global`*"];
sol = 
        DSolve[{x'[t]  ==  
                v[t], v'[t] == -(ω^2) x[t], 
              x[0] == x0, v[0]  ==  v0}, {x[t], v[t]}, t] // 
          FullSimplify;
x[t_] = x[t] /. sol[[1]]
v[t_] = v[t] /. sol[[1]]

{ω, x0, v0} = {0.314, 10, 0};
Plot[x[t], {t, 0, 40}, PlotRange -> {-10, 10}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico spazio-tempo"]
Plot[v[t], {t, 0, 40}, PlotRange -> {-8, 8}, 
      AxesLabel -> {"t (s)", "v (m/s)"}, 
      PlotLabel -> "Grafico velocità-tempo"]

m = 2.5
k = m (ω^2)
Ec[v_] = (1/2) m (v^2)
Ev[x_] = (1/2) k (x^2)
Et[v_, x_] = Ec[v] + Ev[x]

    
Plot[Ev[x[t]], {t, 0, 40}, PlotRange -> {0, 15}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico Energia potenziale-tempo"]
Plot[Ec[v[t]], {t, 0, 40}, PlotRange -> {0, 15}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico Energia cinetica-tempo"]
Plot[{Ec[v[t]], Ev[x[t]], Ec[v[t]] + Ev[x[t]]}, {t, 0, 40}, 
      PlotRange -> {0, 15}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico Energia totale-tempo"]

Oscillazioni smorzate

- k x[t] = -m x''[t]  -  c x'[t]

oppure il sistema
x'[t]  =  v[t]
v'[t]  =  -(ω^2) x[t]  -  2  b x'[t]
(ω^2)  =  k/m
2  b  =  c/m

Clear["Global`*"];
sol = 
      DSolve[{x'[t]  ==  
              v[t], v'[t] == -(ω^2) x[t]  -  
                2 b v[t], x[0] == x0, v[0]  ==  v0}, {x[t], 
            v[t]}, t] // FullSimplify
x[t_] = x[t] /. sol[[1]]
v[t_] = v[t] /. sol[[1]]

{ω, b, x0, v0} = {0.628, 0.1, 10, 
          0};
Plot[{x[t], x[0] Exp[-b t], -x[0] Exp[-b t]}, {t, 0, 
        40}, PlotRange -> {-10, 10}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico spazio-tempo"]
Plot[v[t], {t, 0, 40}, PlotRange -> {-10, 10}, 
      AxesLabel -> {"t (s)", "v (m/s)"}, 
      PlotLabel -> "Grafico velocità-tempo"]

m = 2.5
k = m (ω^2)
Ec[v_] = (1/2) m (v^2)
Ev[x_] = (1/2) k (x^2)
Et[v_, x_] = Ec[v] + Ev[x]

    
Plot[Ev[x[t]], {t, 0, 40}, PlotRange -> {0, 50}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico Energia potenziale-tempo"]
Plot[Ec[v[t]], {t, 0, 40}, PlotRange -> {0, 50}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico Energia cinetica-tempo"]
Plot[{Ec[v[t]], Ev[x[t]], Ec[v[t]] + Ev[x[t]]}, {t, 0, 40}, 
      PlotRange -> {0, 50}, 
      AxesLabel -> {"t (s)", "x (m)"}, 
      PlotLabel -> "Grafico Energia totale-tempo"]

Oscillazioni forzate

- k x[t] = -m x''[t] - c x'[t] + h Cos[Ω t]

oppure il sistema
x'[t]  =  v[t]
v'[t] = -(ω^2) x[t] - 2  b x'[t]+ f Cos[Ω t]
(ω^2)  =  k/m
2  b  =  c/m
f  =  h/m

Soluzione numerica

Clear["Global`*"];
{ω, 
          b, f, Ω, x0, v0} = {0.628, 0.1, 1, 0.65, 10, 
          0};
tmax = 100
{x[t_], v[t_]}={x[t], v[t]}/.NDSolve[{x'[t]  ==  v[t],v'[t]==-(ω^2) x[t] - 2 b v[t] + f Cos[Ω t],x[0] == x0,v[0]  ==  v0},{x[t], v[t]},{t, 0, tmax}][[1]]
Plot[
        x[t], {t, 0, tmax}, PlotRange -> {-10, 10}, 
        AxesLabel -> {"t (s)", "x (m)"}, 
        PlotLabel -> "Grafico spazio-tempo"];
Plot[v[t], {t, 0, tmax}, 
        PlotRange -> {-6, 6}, 
        AxesLabel -> {"t (s)", "v (m/s)"}, 
        PlotLabel -> "Grafico velocità-tempo"];

Soluzione analitica

Clear[sol, ω, Ω, b, f, x, x0, x01, x02, x03, xas, 
      x1, x2, x3, A1, A2, A, B, F, φ, θ, xo, vo, xx, rul1, 
      rul2]
sol = 
      DSolve[x''[t] + 2  b x'[t] + (ω^2) x[t] == 
          f Cos[Ω t], x[t], t]

x01[t_] = x[t] /. sol[[1]] /. 
        f -> 0
xas[t_] = x[t] /. sol[[1]] /. {C[1] -> 0, C[2] -> 0} // 
        Simplify

Se ω < Ω la soluzione dell'omogenea associata diventa oscillante:

x02[t_] = 
      ExpToTrig[
          x01[t] /. {(Sqrt[(b^2) - (ω^2)]) -> I (Sqrt[ (ω^2) - (b^2)]), 
              C[1] -> A (E^(I θ)), 
              C[2] -> 
                A (E^(-I θ))}] // 
        FullSimplify

x1[t_] = x01[t] + xas[t]
x2[t_] = x02[t] + xas[t]

Calcoliamo le regole che permettono di ricavare  A e θ  o C[1] e C[2] dalla condizioni iniziali  x(0)=x0 v(0)=v0

Clear[xo, vo, δ, ω, Ω, f, xx]
rul1 = 
      MapAt[Together, Solve[{x1[0] == xo, x1'[0] == vo}, {C[
                  1], C[2]}][[1]], {{1, 2}, {2, 2}}] // 
        Simplify
sol = Solve[{x2[0] == xo, x2'[0] == vo} /. { 
                  Sin[θ] -> (u/A), 
                  Cos[θ] -> (v/A)}, {u, v}][[1]] // 
          Simplify;
rul2 = {A -> 
          Sqrt[((u^2) + (v^2)) /. sol // Simplify], θ -> 
              ArcTan[Numerator[(v/u)], Denominator[(v/u)]] /. sol // 
          Simplify}

Calcoliamo le soluzioni per b=0 e  b = ω

xx[t_] = 
        xx [t] /. DSolve[{xx''[t] + (ω^2) xx[t] == 
                    Cos[Ω t], xx[0] == xo, xx'[0] == vo}, 
                xx[t], t][[1]] // Simplify;
x3[t_] = 
      Limit[x1[t] /. rul1, ω -> b] // FullSimplify
x03[t_] = 
      Limit[x01[t] /. rul1, ω -> b] // FullSimplify

La soluzione asintotica può essere scrtitta nella forma:

(A1 Sin[Ω t] + 
            A2 Cos[Ω t]) /. {A1 -> 
            F Cos[φ], A2 -> F Sin[φ]} // 
      FullSimplify

Con a e b dati da:

A1 = 
      Simplify[Coefficient[xas[t], 
          Sin[Ω t]]]
A2 = Simplify[Coefficient[xas[t], Cos[Ω t]]]

Ricaviamo infine F e φ

F[f_, b_, ω_, Ω_] = 
      PowerExpand[Simplify[(Sqrt[(A1^2) + (A2^2)])]]
φ[δ_, ω_, Ω_] = 
      ArcTan[Numerator[(A2/A1)], Denominator[(A2/A1)]]

\(Se\ poniamo\ ω = \(1\ e\ f = 1\)\ 

      				Ω\ \ \ \ viene\ misurato\ in\ unità\ di\ ω\ 

      				b\ \ \ \ viene\ misurato\ in\ unità\ di\ ω\ 

      				F\ \ \ \ in\ unità\ f\/ω\^2\)

ω = 1; f = 1;
\(\(Plot[
        Evaluate[
          Table[F[f, \ b\ , ω, Ω\ ], {b, 0.1, 0.6, 
              0.1}]], {Ω, 0, 2.5}, PlotRange -> All, 
        AxesOrigin -> {0, 0}, 
        AxesLabel -> {\*""(Ω/ω)"", \*""((ω^2)/f)F""}];\)\)
\(\(Plot[
        Evaluate[
          Table[φ[b, ω, Ω], {b, 0.1, 0.6, 
              0.1}]], {Ω, 0, 2.5}, 
        PlotRange -> {0, π}, 
        Ticks -> {Automatic, {0, π\/2, π}}, 
        AxesLabel -> {\*""(Ω/ω)"", "φ"}];\)\)

b = ω; Ω = 0.5;
xo = 0;
vo = 1;
x0[t_] := 
        Which[b == ω, x03[t], b > ω, x01[t] /. rul1, 
          b < ω, x02[t] /. rul2];
x[t_] :=  
        Which[b == ω, x3[t], b > ω, x1[t] /. rul1, 
          b < ω, x2[t] /. rul2];
\(\(Plot[{x0[t], xas[t], x[t]}, {t, 0, 4\/b}, 
        PlotStyle -> {RGBColor[0, 0, 1], RGBColor[0, 1, 0], 
            RGBColor[1, 0, 0]}, 
        Ticks -> {Table[α, {α, 0, 4\/b, π}], 
            Automatic}, 
        AxesLabel -> {"ω t", \*""((ω^2)/f)x""}];\)\)

b = 0.1;
Plot[{xx[t]}, {t, 0, (4/b)}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        PlotPoints -> 100];