5ª Aula Prática de TFCOMP-II

Comecem por criar e guardar em disco um Notebook com o nome TFC5_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).

<< Utilities`Notation`

Symbolize[(Subscript[θ, o])]

Problema 1

Faça um estudo do pêndulo simples de comprimento L no campo gravítico g SEM a aproximação sin (θ) ≈ θ, i.e. para movimentos angulares arbitrários.			
Resolva numéricamente a equação diferencial deste pêndulo usando NDSolve com as apropriadas condições iniciais. Faça gráficos representando o tipo de movimentos obtidos com diferentes condições iniciais. 
Calcule o período T deste pêndulo verdadeiro em função da posição inicial parada  (Subscript[θ, o]), usando o facto de que uma oscilação completa leva quatro vezes o tempo de ir de  θ=(Subscript[θ, o]) a θ = 0.
Compare o seu valor com a expressão :  \(TraditionalForm\`T = 
        2\ π\ \@\(L\/g\)\ 1\/\(AG \((1, cos \((θ\_o\/2)\))\)\)\)		
(AG (a, b) a função nativa ArithmeticGeometricMean de a e b) e represente gráficamente o desvio relativamente à fórmula T = 2  π (Sqrt[ (L/g)]) usada para pequenas oscilações.

Solução

Eqn= θ''[t]==-((g Sin[θ[t]])/L);

SolReal=θ[t]/.NDSolve[{θ''[t]==-((g Sin[θ[t]])/L),θ[0] == 1/2,θ'[0] == 0}/.g -> 9.8/.L -> 1,θ[t],{t, 0, 4}];

SolAprox = (Subscript[θ, o ]) 
                Sin[(Sqrt[(g/L)]) t + (π/2)] /. g -> 9.8 /. 
            L -> 1 /. (Subscript[θ, o]) ->  .5;

Plot[{SolAprox, SolReal}, {t, 0, 4}, 
        PlotStyle -> {{Dashing[ .01 {1, 1, 1}], rd}, bl}, 
        Ticks -> {{#  π/4, #  π/4} & /@ Rg[0, 4], 
            Automatic}];

O pêndulo é um sistema conservativo, por isso a energia total Ε = (m/2)L ((θ')^2)- m g L Cos[θ] é constante Ε  =  - 
            m g L Cos[(Subscript[θ, o])].  Daí que a equação para a velocidade angular se deduza de

θEqn=(Solve[{Ε==(1/2) m (((L θ'))^2)-m g L Cos[θ],Ε == -m g L Cos[(Subscript[θ, o])]},θ'][[1, 1]]/.θ -> θ[t]/.θ[t]' -> θ'[t])//At[2, 4, Sp[#, g > 0] &]//PowerExpand//(# /. Rule -> Equal /. 
              Cos[(Subscript[θ, o])] -> a &)

A partir desta equação diferencial linear homogénea obtém-se uma solução para θ[t] na forma de uma função inversa: θ = InverseFunction[f][t] é o valor que verifica f[θ] = t.

θSol = 
      DSolve[{θEqn, θ[0] == (Subscript[θ, o])}, θ[t], 
        t]

θSol = (θSol[[1, 
              1]] // 
            At[-1, -1, Sp[#, a < 1] &]) /. 
        a -> Cos[(Subscript[θ, o])]

O tempo correspondente à primeira passagem por θ = 0 é um quarto de período (Subscript[T, 1/4])(para um movimento oscilatório), por isso

θPeriod = 
      Solve[θ[t] == 0 /. θSol, t] /. 
          t -> (Subscript[T, 1/4]) // Fl

Para os valores indicados L = 1,  g = 9.8 tem-se (Subscript[ω, o]) = (Sqrt[g/L])  ->  
            (Subscript[T, (Subscript[ω, o])]) = ((2  π)/(Subscript[ω, o]))

<< Graphics`Legend`

Plot[{2  π (Sqrt[ (L/g)]) 
                    (1/ArithmeticGeometricMean)[1, Cos[((Subscript[θ, o])/2)]], 4 
                           (Subscript[T, 1/4]) /. θPeriod  /. 
                      g -> 9.8 /. L -> 1, 
                  2  π (Sqrt[ (L/g)]) } /. g -> 9.8 /. L -> 1 // 
            Re // Release, {(Subscript[θ, o]),  0, π}, 
        PlotPoints -> 30, PlotRange -> {0, 6}, 
        PS -> {{Thickness[0.017], gr}, {Thickness[0.007], rd}, GL[0]}, 
        Ticks -> {{#  π/4, #  π/4} & /@ Range[0, 4], 
            Automatic}, 
        GridLines -> {#  π/4 & /@ Range[0, 4], 
            None}, 
PlotLegend -> {(Subscript[T, AG]), 
            4  (Subscript[T, 1/4]), (Subscript[T, (Subscript[ω, o])])}, LegendSize ->  .5, 
        LegendPosition -> {1, 0}];

Problema 2

Analise o sistema do problema anterior no espaço de fases com coordenadas {θ,(OverDot[ θ])}, integrando um sistema de duas equações diferenciais de 1ª ordem. Represente gráficamente o movimento neste espaço. Visualise os três tipos de curva que pode obter e interprete-as.		
Determine a velocidade angular (Subscript[(OverDot[θ]), o]) na posição inicial  (Subscript[θ, o]) necessária  para que o pêndulo atinja a posição de equilíbrio instável  θ = π. Despreze atritos no movimento.

Solução

Symbolize[\(θ\& . \)]

T[θ_] = 2  π (Sqrt[ (L/g)]) 
              (1/ArithmeticGeometricMean)[1, Cos[(θ/2)]] /. 
            L -> 1 /. g -> 9.8;

phase[θinit_, vinit_]:=NDSolve[{θ'[t]==\(θ\& . [t]\),(\(θ\& . \))'[t]==-9.8` Sin[θ[t]],θ[0] == θinit,\(θ\& . [0] == vinit\)},\({θ[t], θ\& . [t]}\),{t, -(T[θinit] // Abs), 
                T[θinit] // Abs}];

\(\(\(Table[\((pp[{θ[t], θ\& . [t]} /. 
                    phase[θ\_o, 
                      v\_o], {t, \(-\((T[θ\_o] // Abs)\)\), 
                    T[θ\_o] // Abs}, 
                  DisplayFunction -> 
                    Identity, 
PlotRange -> {{\(-2\) π, 2  π}, {\(-9.5\), 9.5}}, 
AspectRatio -> 1, 
                  Frame -> True, 
                  AxesLabel -> {θ, θ\& . }])\), {θ\_o, \(-6\), 6, 2}, {v\_o, \(-7\), 7, 2. }] /. 
            pp -> ParametricPlot // Fl\) // SH[#, Dspl] &;\)\)

\(Table[\((ppp[{Cos[θ[t]], 
                        Sin[θ[t]],  .5\ θ\& . [t]} /. 
                      phase[θ\_o, v\_o] // 
                    Fl, {t, \(-\((T[θ\_o] // Abs)\)\), 
                    T[θ\_o] // Abs}, DisplayFunction -> Identity, 
                  BoxRatios -> {1, 1, 1}, Boxed -> False, 
                  PlotRange -> {{\(-π\), π}, {\(-π\), π}, {\(-2\) π, 2  π}}, 
                  PlotPoints -> 220])\), {θ\_o, \(-6\), 6, 
                2}, {v\_o, \(-7\), 7, 2. }] /. 
            ppp -> PointParametricPlot3D // Fl\)//SH[#/.Graphics3D[z__] :> 
                    Graphics3D[{PointSize[ .003], {z}[[1]]}, Drop[{z}, 1]],DisplayFunction -> $DisplayFunction]&;

Pode-se experimentar uma visualizaçáo dinâmica com

RealTimeGraphicsDisplayFunction[%,  Mesh -> False, 
      Boxed -> True, BoxRatios -> {10,  .1, 1. }]

A velocidade angular (Subscript[(OverDot[θ]), o]) na posição inicial  (Subscript[θ, o]) necessária  para que o pêndulo atinja a posição de equilíbrio instável  θ = π é

\(Solve[\(1\/2\) m\ \(L\^2\) \(θ\& . \)\_o\^2 - 
          m\ g\ L\ Cos[θ\_o] == \(-m\)\ g\ L\ Cos[π], \(θ\& . \)\_o]\)

Problema 3

Faça um desenvolvimento em série de Taylor da função nativa Hypergeometric2F1[a,  b,  c,  x] na vizinhança de (Subscript[x, o]) = 0 e deduza a expressão para o termo geral (Subscript[A, k ]) (x^k).  

Use a função Limit para calcular o raio de convergência desta série, i.e. o valor R para o qual se verifica que a série converge quando |x| < R e diverge quando |x| > R. 

Determine a condição sobre a,  b e c para que a série seja convergente quando x = R.

Para examinar o comportamento da série quando  x = R use o teste de Gauss: uma série de termos positivos (∑_k) (Subscript[A, k]) é convergente se, a partir de um valor k ≥ N, 

                                                        ((Subscript[A, k + 1])/(Subscript[A, k])) ≈ 
            1 - (α/k) + β[k]/(k^(1 + δ))
onde α > 1 e δ > 0 são constantes e ({{ }, {lim}, {k -> ∞}}) β[k] < ∞.

Use as funções OpenWrite e Write para armazenar num ficheiro "hyper.dat" os valores reais com 5 dígitos da função Hypergeometric2F1[2, 0, 4, x] avaliada em 256 pontos no domínio de convergência |x| < R.  Use Read para fazer um ListPlot desses dados.

Problema 4

Construa um operador genérico HamEq[ham_, vars_]  que, dada uma função ℋ[(Overscript[q, ⇀]),  
        (Overscript[p, ⇀])] de coordenadas generalizadas (Overscript[q, ⇀]) = {(Subscript[q, 1]), …,  (Subscript[q, n])} e respectivos momentos (Overscript[p, ⇀]) = {(Subscript[p, 1]), …,  (Subscript[p, n])}, deduza as componentes do respectivo campo vectorial hamiltoneano (Subscript[(Overscript[X, ⇀]), ℋ])={\(q\& . \_1\),…, \(q\& . \_n\),(Subscript[(OverDot[ p]), 1]),…,(Subscript[(OverDot[ p]), n])}, cujas componentes obedecem às Equações de Hamilton  
                                                    ({{((d(Subscript[q, 1]))/(dt))=((∂ℋ)/(∂(Subscript[p, 1])))}, {⋮}, {((d(Subscript[q, n]))/(dt))=((∂ℋ)/(∂(Subscript[p, n])))}})    e       ({{((d(Subscript[p, 1]))/(dt))=-((∂ℋ)/(∂(Subscript[q, 1])))}, {⋮}, {((d(Subscript[p, n]))/(dt))=-((∂ℋ)/(∂(Subscript[q, n])))}})
                                                    
Use as funcionalidades de  Graphics`PlotField`  para representar o campo hamiltoneano para funções de duas variáveis, e em particular aplique-o ao Hamiltoneano (i.e. função Energia Total) do pêndulo simples do primeiro problema.

Solução

HamEq[ham_, vars_] := 
      Block[{n = Length[vars]/2, 
          rls = Thread[vars -> Through[vars[t]]]}, 
        Table[{(∂_t)(vars[[
                      i]][
                    t]) == (((∂_(vars[[i + n]]))ham) /. 
                  rls), (∂_t)(vars[[
                      i + n]][
                    t]) == -(((∂_(vars[[i]]))ham) /. rls)}, {i, 1, 
            n}]]

\(ℒ = \(1\/2\) 
          m\ \((L\ \(θ\& . \))\)\^2 - \((\(-m\)\ g\ L\ Cos[θ])\)\)

Symbolize[(Subscript[p, θ])]

\(p\_θ == ∂\_\(θ\& . \)ℒ\)

\(momentum = 
      Solve[p\_θ == ∂\_\(θ\& . \)ℒ, \(θ\& . \)] // Fl\)

\(ℋ = \(1\/2\) m\ \((L\ \(θ\& . \))\)\^2 - 
          m\ g\ L\ Cos[θ] /. momentum\)

HamEq[-g L m Cos[θ] + 
          ((Subscript[p, θ]^2)/(2 (L^2) m)), {θ, (Subscript[p, θ])}] // 
      TF

<< Graphics`PlotField`

PlotVectorField[{((Subscript[p, θ])/((L^2) m)), -g L m Sin[θ]} /. L -> 1 /. m -> 1 /. 
        g -> 9.8, {θ, -1.5 π, 
        1.5  π}, {(Subscript[p, θ]), -12, 12}, AR -> 1, 
      Frame -> True, {ScaleFactor -> None, 
        MaxArrowLength -> 4.4`, ScaleFunction -> ( .1 #1 &), 
        HeadScaling -> Absolute, HeadWidth -> 0.41`, 
        HeadLength -> 5.51`, PlotPoints -> 30, 
        ColorFunction -> (RGBColor[#1,  .1  #1, 
                1 - #1] &) }]