6ª Aula Prática de TFCOMP-II

Comecem por criar e guardar em disco um Notebook com o nome TFC6_grupo.nb, onde 'grupo' é o nome do vosso grupo de trabalho. Neste Notebook criem uma célula de título identificando os elementos constituintes do grupo (nome, número e contactos).

Inicialização

<< Graphics`PlotField`

({{({{At[pos__, 
                      op_] := #1 /. #1[[
                            pos]] :> 
                          op[#1[[
                              pos]]] &}})}})

Problema 1

Considere uma roda de raio R e velocidade V que roda sem deslizar num movimento linear e uniforme . Deduza as expressões para a velocidade (OverVector[v]) e aceleração (OverVector[a]) de qualquer ponto da roda, quando visto do referencial do chão. Faça gráficos usando PlotVectorField mostrando uma distribuição destes vectores em todas as zonas da roda.

Solução 1

Inicialização

openlnk := 
        Install["/home/amaro/MATHEMATICA/MathLink/Arrow/arrowtip"];
openlnk;

Clear[Vect]Clear[Vect]

Vect[v_,x_,hs_:.1,s_:1] :=  {Line[{{0,0},s v}+Table[x,{i,1,2}]], 
				Polygon[Partition[Arrow[s v,hs],2]+Table[x,{i,1,6}]]}Vect[v_, x_, hs_:0.1, s_:1] := 
 
  {Line[{{0, 0}, s v} + Table[x, {i, 1, 2}]], 
 
   Polygon[Partition[Arrow[s v, hs], 2] + Table[x, {i, 1, 6}]]}

<< Utilities`Notation`

Symbolize[\(e\&⇀\__\)]

Symbolize[\(_\&⇀\)]

Symbolize[(Overscript[r, ⇀])]

Symbolize[\(_\_o\)]

A posição dos pontos da roda é muito convenientemente indicada por ρ (Subscript[((Overscript[e,  ⇀]) ), ρ]), com ρ ∈ [0, R] e  \(TraditionalForm\`\(e\& ⇀ \)\_ρ, \ \(e\& ⇀ \)\_θ\) vectores polares ligados ao centro da roda. No caso de movimento uniforme, sem deslizamento, com velocidade (Overscript[V,  ⇀])  = 
        V (Subscript[((Overscript[e,  ⇀]) ), x]) relativa ao chão, a posição do centro é dada por (Overscript[R,  ⇀])  = 
        V t (Subscript[((Overscript[e,  ⇀]) ), x]) + R (Subscript[((Overscript[e,  ⇀]) ), y]). Assim, o vector posição em qualquer instante dum ponto da roda é, no referencial do chão,  \(TraditionalForm\`\(r\& ⇀ \) = \(R\& ⇀ \) + ρ\ \(e\& ⇀ \)\_ρ\). A velocidade e aceleração desse ponto determinam-se como sempre calculando \(d\ \(r\& ⇀ \)\/dt\) e (((d^2)(Overscript[r,  ⇀]) )/(d(t^2))), tendo em conta que para um ponto fixo da roda \(\(ρ\& . \)\  = \ \(\(ρ\&¨\)\(\ \)\(=\)\)\) 0, e se esta rola com velocidade angular constante \(θ\& . \) = -(V/R) então (Overscript[θ, ¨])  =  0.

ClearAll[(Overscript[r,  ⇀]) , (Overscript[v,  ⇀]) , (Overscript[a,  ⇀]) ]
\(\(e\& ⇀ \)\_ρ = {Cos[θ], 
        Sin[θ]}; 		\(e\& ⇀ \)\_θ = {\(-Sin[θ]\), Cos[θ]};\)
\(\(R\& ⇀ \) = {V\ t, 
        R}; 				\ \(r\& ⇀ \) = \(R\& ⇀ \) + ρ\ \(e\& ⇀ \)\_ρ;\)
\(\(\(\(v\& ⇀ \) = 
      Dt[\(r\& ⇀ \), 
        t]\)\(;\)\(\ 			\ \ \ \ \)\(\(a\& ⇀ \) = 
      Dt[\(v\& ⇀ \), 
        t]\)\(;\)\(
\)\)\)

rva=Thread[{HoldForm[(Overscript[r, ⇀])], 
              HoldForm[(Overscript[v, ⇀])], HoldForm[(Overscript[a, ⇀])]}==({(Overscript[r, ⇀]), (Overscript[v, ⇀]), 
                  (Overscript[a, ⇀])}	/.	({{\(Dt[ρ, t] -> \(ρ\& . \)\)}, {\(Dt[θ, t] -> \(θ\& . \)\)}, {Dt[ρ, {t, 2}] -> (Overscript[ρ, ¨])}, {Dt[θ, {t, 
                              2}] -> (Overscript[θ, ¨])}, {Dt[V, t] -> (Subscript[A, x])}, {Dt[R, _] -> 0}, {Dt[V, {t, 2}] -> 0}}) 	)]//TableForm

rva	/.	 ({{\(θ\& . \)->-(V/R)}, {\(\(_\& . \) -> 0\)}, {\(_\&¨ -> 0\)}, {(Subscript[A, x]) -> 0}, {Sin[θ]->((-R + y)/ρ)}, {Cos[θ]->((-t V + x)/ρ)}})//Simplify//TableForm
ClearAll[(Overscript[r, ⇀]), 
      (Overscript[v, ⇀]), 
      (Overscript[a, ⇀])]
(%% // ReleaseHold) /. 
        Equal -> Set;

Podemos agora integrar as equações de movimento para um ponto genérico {x0, y0} da roda, usando DSolve. Depois de obter a função temporal POSITION, as suas derivadas explícitas em ordem a t determinam VELOCITY e ACCELERATION em cada POSITION.

\(POSITION = \({x[t], y[t]} /. 
            DSolve[{Dt[\(r\& ⇀ \), {t, 
                            2}] == \(a\& ⇀ \)\ \  /. {x -> x[t], y -> y[t]} // 
                    Thread, 
\ \ \ 					\(Dt[\(r\& ⇀ \), 
                            t] == \(v\& ⇀ \)\ \ \ \ \ \  /. {x -> x[t], y -> y[t]}\) /. t -> 0 // 
                    Thread, 
\ \ \ 					\(\(r\& ⇀ \) == {x\_o, y\_o}\ \ \ \ \ \  /. {x -> 
                            x[t], y -> y[t]}\) /. t -> 0 // 
                    Thread} // Flatten, {x[t], y[t]}, t] // Flatten\) // 
        Simplify\)

\(VELOCITY = ∂\_t POSITION			 // Simplify\)

\(ACCELERATION = ∂\_\(t, t\)POSITION	\ \  // 
        Simplify\)

No seguinte gráfico, estão representadas as trajectórias (ciclóides) dos pontos da roda inicialmente sobre a vertical do ponto de contacto no chão. As respectivas velocidades (a cores) e acelerações (a preto) mostram que a aceleração é sempre radial no caso do movimento uniforme.

\(\(\((Show[#1, DisplayFunction -> \((Display[$Display, #1] &)\), 
              Axes -> {True, False}] &)\)[\(\(({Append[
                  Table[
ppl[
                        POSITION, {t, 0, 
                          6}, 
PlotPoints -> 30, 
                        DisplayFunction -> Identity, 
                        AspectRatio -> 2.5/8, 
                        PlotRange -> {{\(-1\), 7}, {0, 2.5}}, 
                        PlotStyle -> TH[ .4]] //. {x\_o -> 
                          r\ Cos[θ], 
                        y\_o -> R - r\ Sin[θ], V -> 1, 
                        R -> 1}, {θ, \(-\(π\/2\)\), 
                      2  π, π}, {r, 0, 1,  .2}], 
                  ppl[{R\ Cos[θ], R - R\ Sin[θ]} /. 
                      R -> 1, {θ, 0, 2  π}, 
                    DisplayFunction -> Identity]], 
                Table[Graphics[{Hue[\(r\ θ\)\/\(2\ π\)], 
                        vect[ .5\ VELOCITY, POSITION, 
                          0.05]}]			 //. {x\_o -> r\ Cos[θ], 
                      y\_o -> R - r\ Sin[θ], V -> 1, 
                      R -> 
                        1}, {θ, \(-\(π\/2\)\), π\/2, π}, {r, 0, 1,  .2}, {t, 0, 6,  .5}], 
                Table[Graphics[{vect[\ 1\ ACCELERATION, POSITION, 
                          0.05]}]	 //. {x\_o -> r\ Cos[θ], 
                      y\_o -> R - r\ Sin[θ], V -> 1, 
                      R -> 
                        1}, {θ, \(-\(π\/2\)\), π\/2, π}, {r, 0, 1,  .2}, {t, 0, 6,  .5}]
})\) /. 
            vect -> Vect\) /. ppl -> ParametricPlot];\)\)

Uma 'fotografia' da roda daria uma imagem instantânea das velocidades e acelerações da seguinte forma:

wheel = 
        ParametricPlot[{R Cos[θ], R - R Sin[θ]} /. 
              R -> 1 // Release, {θ, 0, 2  π}, 
          PlotRange -> {{-1, 1}, {0, 2}}, AspectRatio -> 1, 
          DisplayFunction -> Identity];

Show[GraphicsArray[{
{Graphics[Table[If[(# . # &@(POSITION - (Overscript[R,  ⇀]) )) ≤ (R^2), 
                                    vect[0.2` VELOCITY, POSITION, 
                                        0.03] // 
                                      Release, {}]//.({{(Subscript[x, o]) -> 
                                        t V + r Cos[θ]}, {(Subscript[y, o]) -> r Sin[θ]}, {V -> 1}, {R -> 1}, {t -> 0}}),
{θ, 0, 
                                  2  π, (π/30.)`},{r, 0, 2, 0.1`}],
PlotRange -> {{-1, 1}, {0, 2}},AspectRatio -> 1,DisplayFunction -> Identity,PlotLabel -> "Velocidade"],wheel}/.vect -> Vect//Show,
{Graphics[Table[If[(# . # &@(POSITION - (Overscript[R,  ⇀]) )) ≤ (R^2), 
                                    vect[0.2` ACCELERATION, POSITION, 
                                        0.03`] // Release, {}]//.({{(Subscript[x, o]) -> 
                                        t V + r Cos[θ]}, {(Subscript[y, o]) -> 
                                        R + r Sin[θ]}, {V -> 1}, {R -> 1}, {t -> 0}}),
{θ, 0, 
                                  2  π, (π/30.)`},{r, 0, 2, 0.1`}],
PlotRange -> {{-1, 1}, {0, 2}},AspectRatio -> 1,DisplayFunction -> Identity,PlotLabel -> "Aceleraç\
ão"],wheel}/.vect -> Vect//Show}],DisplayFunction -> (Display[$Display, #1] &)];

Problema 2

Defina uma função EulerSolve para obter listas de valores de soluções aproximadas de equações diferenciais ordinárias.

Note que uma equação diferencial de 1ª ordem ((dx)/(dt))=f (x, t)de facto pode ser entendida como um campo de direcções, representado pelo versor (Subscript[(OverVector[e]), f])=\({1\/\@\(\(f \((x, t)\)\)\^2 + 1\), \(f \((x, t)\)\)\/\@\(\(f \((x, t)\)\)\^2 + 1\)}\)

Construa uma representação do campo em que as linhas de campo sejam usadas para origem dos vectores, i.e. a lista de pontos obtida por uma sequência de condições iniciais (Subscript[x, o]), (Subscript[t, o]) suficientemente grande para cobrir a área a representar.

Considere a seguinte equação diferencial: ((dx)/(dt))=t - (x^2). Utilize a função EulerSolve para estudar as suas soluções para vários (Subscript[x, o]) na vizinhança de t = 0.  Faça um gráfico usando ListPlot e PlotVectorField em que seja visível o comportamento assimptótico das várias soluções. Verifique que para t suficientemente grande o método de Euler falha, obtendo-se um comportamento aparentemente caótico dos pontos calculados. Mostre que existe um passo δt suficientemente pequeno para que EulerSolve volte a funcionar, dando soluções perfeitamente regulares.

Solução

m = 10;
\(\(clr := {Off[Unset::"norep"]; 	
\ \ ClearAll[t, 
            x]; 	x\_n_ =. ; 	t\_n_ =. ; 
	
          ReleaseHold[\((HoldForm[x\_#1 =. ] &)\) /@ Range[0, m]]; 
	
          ReleaseHold[\((HoldForm[t\_#1 =. ] &)\) /@ 
              Range[0, m]];};\)\)

\(EulerSolve[p0_List, δt_ :  .1, 
        k_ : 10] := 
		Block[{}, 	m = k; 	
        clr; 
			x\_0 = 
          p0[[1]]; 	
        t\_0 = p0[[2]]; 
			t\_n_ := \(t\_n = 
            t\_\(n - 1\) + δt\); 	
			x\_n_ := \(x\_n = x\_\(n - 1\) + 
              f[x\_\(n - 1\), 
                  t\_\(n - 1\)]\ δt\); 
			elist = \({t\_#1, x\_#1} &\) /@ \ Range[0, m]
			] (*\ 
      kcolB\ *) 
\)
\(\(lp\_n_ := 
        ListPlot[\(({t\_#1, x\_#1} &)\) /@ 
            Range[0, n], 
			DisplayFunction -> 
            Identity, 
          Axes -> 
            True, 
			PlotStyle -> {Thickness[
                0.005], RGBColor[1, 0, 0]}];\)\)

O método de Euler melhorado seria:

\(\(\(EulerSolve2[p0_List, δt_ :  .1, k_ : 10] := 
      Block[{}, 	
			m = 
          k; 		\ clr; 
			x\_0 = 
          p0[[1]]; 	
        t\_0 = p0[[2]]; 

			t\_n_ := \(t\_n = 
            t\_\(n - 1\) + δt\); 
			x\_n_ := \(x\_n = x\_\(n - 1\) + 
              1\/2\ \((f[
                      x\_\(n - 1\) + 
                        f[x\_\(n - 1\), t\_\(n - 1\)]\ δt, t\_n] + 
                    f[x\_\(n - 1\), 
                      t\_\(n - 1\)])\)\ δt\); 	
			elist = \({t\_#1, x\_#1} &\) /@ 
            Range[0, m]
			]\)\( (*\ 
      kcolB\ *) \)\)\)

Exemplo:

f[x_, t_] := -5 t (x^2)

EulerSolve[{1, 0},  .05, 30];

\(\(Show[{lp\_20, 
          PlotVectorField[\(#\/\@\(# . #\) &\)@{1, f[x, t]}, {t, 0, 2}, {x, 
              0, 1}, DisplayFunction -> Identity, 
            ScaleFunction -> None, ScaleFactor ->  .05, 
            HeadScaling -> Relative, HeadWidth ->  .5, 
            HeadLength ->  .3]}, 
        DisplayFunction -> \((Display[$Display, #1] &)\), 
        AspectRatio ->  .5];\)\)

f[x_, t_] := t - (x^2)

\(\((\(ListPlot[EulerSolve2[{#, 0},  .01, 550], NoDspl, 
              PlotRange -> {\(-2\), 5}] &\) /@ Rg[\(- .9\), 2,  .3])\) // 
      SH[{#, PlotVectorField[\(#\/\@\(# . #\) &\)@{1, f[x, t]}, {t, 0, 
                5}, {x, \(-1.9\), 5}, DisplayFunction -> Identity, 
              ScaleFunction -> None, ScaleFactor ->  .2, 
              HeadScaling -> Relative, HeadWidth ->  .5, 
              HeadLength ->  .3, ColorFunction -> \((Hue[#] &)\), 
              PlotPoints -> 25]}, 
          DisplayFunction -> \((Display[$Display, #1] &)\)] &\)

Instabilidade numérica cria um aparente mecanismo de duplicação de período conducente a um movimento caótico, mas esse comportamento é um artifício resultante da falta de precisão do método de integração numérica: a diminuição do passo δt resolve essa situação dentro do intervalo de tempo considerado.

EulerSolve[{0, 0},  .04, 150000. ];

ListPlot[elist, Dspl, 
      PlotJoined -> 
        False, 
PlotStyle -> {PointSize[0.000000001], 
          RGBColor[1, 0, 0]}, 
      PR -> {{0, 1420}, {-40, 70}}]

EulerSolve[{0, 0},  .01, 150000. ];

ListPlot[elist, Dspl, 
      PlotJoined -> 
        False, 
PlotStyle -> {PointSize[0.000000001], 
          RGBColor[1, 0, 0]}, 
      PR -> {{0, 1420}, {-40, 70}}]

Problema 3

Mostre que a expansão da solução exacta x[t] duma equação diferencial ((dx)/(dt))==f (x, t) até à ordem 2 em δt na vizinhança de (Subscript[x, o]) coincide com um passo do método de Euler melhorado a menos de termos de ordem 𝒪[δ(t^3)]. Mostre que as mesmas expansões levadas  até à ordem 4 validam o método de Runge-Kutta de ordem 4.

Problema 4

Construa um operador genérico LagrEq[lagr_, vars_]  que, dada uma função lagrangeana ℒ[(Overscript[q, ⇀]), 
        (Overscript[v, ⇀])] de coordenadas generalizadas (Overscript[q, ⇀]) = {(Subscript[q, 1]), …,  (Subscript[q, n])} e respectivas velocidades (Overscript[v, ⇀]) = {(Subscript[v, 1]), …,  (Subscript[v, n])}, deduza explícitamente as equações de Lagrange

                                                     \(TraditionalForm\`\(d\/dt\) \((∂ℒ\/∂v\_i)\) - ∂ℒ\/∂q\_i = 𝒬\_i\)      

onde (Subscript[𝒬, i])representam forças generalizadas.   Aplique-o para deduzir equações de movimento para o pêndulo de Atwood, consistindo em duas massas (Subscript[m, 1]) e (Subscript[m, 2]), ligadas por um fio inextensível passando livremente por um curto tubo horizontal fixo. A massa (Subscript[m, 2]) está inicialmente em equílibrio, enquanto (Subscript[m, 1]) tem uma condição inicial diferente do equilíbrio. Integre numéricamente as equações e represente gráficamente o movimento do pêndulo.

Solução

<< Utilities`Notation`

Symbolize[\(_\& . \)]
Symbolize[\(_\__\)]

O lagrangeano deste pêndulo, com a restrição de movimentos planos da massa (Subscript[m, 1]) e movimentos verticais na massa (Subscript[m, 2]), é função de coordenadas generalizadas  ℓ, θ e suas velocidades \(TraditionalForm\`\(ℓ\& . \), \(θ\& . \)\).

\(\(ℒ = \(1\/2\) \((\((m\_1 + 
                      m\_2)\) ℓ\& . [t]\^2 + \(m\_\(\(1\)\(\ \)\)\) ℓ[t]\^2\ θ\& . [t]\^2)\) + \(m\_1\) 
            g\ ℓ[t]\ Cos[θ[t]] + \(m\_2\) 
            g \((L - ℓ[t])\);\)\)

As respectivas equações de Euler-Lagrange  ((d(\(∂\_\(ℓ\& . \)ℒ\)))/(dt))-((∂ℒ)/(\(∂ℓ\)))==0,((d(\(∂\_\(θ\& . \)ℒ\)))/(dt))-((∂ℒ)/(\(∂θ\)))==0 fornecem as equações de movimento

(Eqns=(\({∂\_t\((∂\_\(ℓ\& . [t]\)ℒ)\) - ∂\_\(ℓ[t]\)ℒ == 0, 
∂\_t\((∂\_\(θ\& . [t]\)ℒ)\) - ∂\_\(θ[t]\)ℒ == 0}\)//
Solve[#, {(\(ℓ\& . \))'[t],(\(θ\& . \))'[t]}]&//
Simplify//
ReplaceAll[#, 
                  Rule -> Equal] &))//TableForm

A função atw gera soluções numéricas destas equações a partir de (Subscript[θ, o]),  (Subscript[m, 1]),  (Subscript[m, 2]) para um pêndulo cujo comprimento inicial do lado de (Subscript[m, 1]) é ℓ = 1.

\(atw[θi_, m1_, m2_] := 
      Block[{m\_1 = m1, 
          m\_2 = m2}, ℓ[
                t] {Sin[θ[t]], \(-Cos[θ[t]]\)} /. 
            NDSolve[{Eqns /. g -> 9.8 // 
                    Release, \(ℓ'\)[t] == ℓ\& . [
                      t], \(θ'\)[t] == θ\& . [
                      t], θ[0] == θi, ℓ[
                      0] == 1, \ ℓ\& . [0] == 
                    0, \ θ\& . [0] == 2} // 
                Fl, {θ[t], ℓ[t], θ\& . [
                  t], ℓ\& . [t]}, {t, 0, 
                24}, 
MaxSteps -> 5000] // 
          Fl]\)

ParametricPlot[atw[π  .2, #, 1] // Release, {t, 0, 24}, 
            PlotRange -> {{-1, 1}, {-1.5,  .5}}, 
            PlotPoints -> 1500] & /@ Rg[ .1,  .9,  .05];

mtn = atw[π  .1,  .1, 1] // Release;

mvie = 
        Table[ParametricPlot[mtn, {t, k, k +  .1}, 
            PlotRange -> {{-1, 1}, {-1.5,  .5}}, 
            AspectRatio -> 1, DisplayFunction -> Identity, 
            Ticks -> None, PlotPoints -> 5], {k, 0, 
            23,  .1}];

MakeGifMovie[
      mvie[[#]] &, "movie3.gif", 1, 
      231]

<< "/home/amaro/TFCOMP/GMT/GifMovieTool.m"

Considere um cenário epidémico em que uma fracção x duma população tem uma doença infecciosa. Se não se tomarem medidas de quarentena, de forma que os membros da população contactam livremente, deduza a taxa de aumento da infecção d(x/d)t, supondo que existe uma probabilidade p de que cada encontro entre um doente e uma pessoa saudável resulte num contágio. Estude as soluções dessa equação e tire conclusões a respeito do futuro da população. Faça gráficos suportando a sua análise para várias condições iniciais.