8ª Aula Prática de TFCOMP-II

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

Problema 1

Usando os nomes de dois elementos do seu grupo, crie duas imagens com letras do tipo {"Times-Bold", 64} e transforme-as em duas matrizes quadradas de números reais (use as funcionalidades de Display e/ou Export). Usando Fourier e InverseFourier,   represente com ListDensityPlot ou Raster o efeito nas imagens de alterar a fase ou amplitude dos respectivos coeficientes de Fourier. Observe também a convolução e a correlação destas imagens uma com a outra. Que conclusões tira sobre a fase e a amplitude dos coeficientes de Fourier?

Problema 2

Determine a solução da equação de difusão (((∂^2) u)/(∂(x^2)))=((∂u)/(∂t)) com a condição inicial u[x, 0] = 1 se \(TraditionalForm\`|x| ≤ 1\) e zero fora deste intervalo. Faça uma animação mostrando a evolução da solução para t > 0.

Solução 2

Solução usando transformadas de Fourier

Para a solução deste problema veja o documento incluído na página da cadeira sob o título "Equações de Difusão e Fourier"

Condição inicial e sua transformada.

f[x_] = If[Abs[x] ≤ 1, 1, 0];

ℱ[
          k] == (F[k] = FourierTransform[f[x], x, k]) // 
      FullSimplify

Transformada da Equação de Fourier (atenção ao HoldForm: sem ele o resultado é errado)

HEqn = (FourierTransform[(∂_(x, x))u[x, t], x, k] == 
            HoldForm[(∂_t) FourierTransform[u[x, t], x, k]]) /. 
        FourierTransform[u[x, t], x, k] -> 𝒰[k, 
            t]

FSol = 
      DSolve[HEqn // ReleaseHold, 𝒰[k, t], {k, t}, 
          DSolveConstants -> (ℱ &)] // 
        Flatten

A solução obtém-se pela transformada inversa de 𝒰[k, t] quando se substitui ℱ[k] pela transformada da condição inicial.

u[x, t] == 
        InverseFourierTransform[𝒰[k, t] /. 
              FSol /. ℱ -> F, k, x] // Sp

Table[
        Plot[(1/2) (Erf[((Sqrt[(((-1 + x))^2)])/(2 (Sqrt[t])))] Sign[1 - x] + 
                  Erf[((Sqrt[(((1 + x))^2)])/(2 (Sqrt[t])))] Sign[1 + x]) // 
            Release, {x, -5, 5}, PlotRange -> {0, 1.5}], {t, 
          0.0001,  .15,  .05}];

Solução por convolução com o "Heat Kernel" h[x, t] da condição inicial f[x]

h[x_, t_] = 
      InverseFourierTransform[(E^(-4 (k^2) (π^2) t)), 
        k, x]

u[x, 
          t] == "h*f" == (∫_(-∞)^∞) 
              h[x - ξ, t] f[ξ] \[DifferentialD]ξ // 
      Sp

Table[
        Plot[(1/2) (-Erf[((-1 + x)/(2 (Sqrt[t])))] + 
                  Erf[((1 + x)/(2 (Sqrt[t])))]) // Release, {x, -5, 5}, 
          PlotRange -> {0, 1.5}], {t, 0.0001,  .15,  .05}];

Problema 3

Uma corda de guitarra está presa em x = 0 e x = L. Sugira condições iniciais u[x, 0] = f[x] e ((∂u)/(∂t))[x, 0]=g[x] e, usando transformadas de Fourier e/ou de Laplace, resolva a equação  de ondas (((∂^2) u)/(∂(x^2)))==(1/(v^2)) (((∂^2) u)/(∂(t^2))), onde v = (Sqrt[(T/λ)]), T a tensão na corda e λ a sua densidade linear de massa. Mostre a evolução da solução numa animação e use a função Play para ouvir o som gerado pela corda, e compare a intensidade do tom fundamental com os das harmónicas seguintes .

Solução 3

A solução geral de uma equação do tipo (((∂^2) u)/(∂(x^2)))==(1/(v^2)) (((∂^2) u)/(∂(t^2))) é u[x, 
          t] = 𝒻[x - v t] + ℊ[x + v t], onde 𝒻[x]e ℊ[x] são funções arbitrárias, representando perfis que se deslocam para a direita e a esquerda com velocidade v. De facto:

HoldForm[(∂_(x, x))ψ  - (1/(v^2)) (∂_(t, t))ψ == 
            0] /. ψ -> 𝒻[x - v t] + ℊ[x + v t] // 
      ReleaseHold

% // Simplify

Quando se especificam condições iniciais, restringe-se a forma de 𝒻 e ℊ.

u[x, 0] == 𝒻[x] + ℊ[x] == 
      f[x]
((∂_t) u[x, 
                  t] == (∂_t)(𝒻[
                      x - v t] + ℊ[x + v t]) /. 
            t -> 0) == g[x] // Simplify

Solve[{v (-𝒻'[x]+ℊ'[x])==g[x],𝒻'[x]+ℊ'[x]==f'[x]},{𝒻'[x],ℊ'[x]}]/.Rule -> Equal//Flatten

sol=DSolve[{𝒻'[x]==-g[x]/(2 v)+((f'[x])/2),ℊ'[x]==g[x]/(2 v)+((f'[x])/2)},{𝒻, ℊ},x];
sol = 
      sol /. Thread[(Cases[sol, 
                    z_ /; StringMatchQ[
                        ToString[z], "K$*"], {0, ∞}] // 
                  Union) -> {η, ξ, η, ξ}] // 
        Flatten

Vemos assim que a solução que verifica as condições iniciais é, àparte as constantes, u[x, t]==(1/2) (f[x + v t] + f[x - v t])+(1/(2  v))((∫_ξ^(x + v t))\(g[η] dη\)-(∫_ξ^(x - v t))\(g[η] dη\))

\(\(\(\(u[x, t] == 𝒻[x - v\ t] + ℊ[x + v\ t]	 /. 
              sol\)	 /. 
Integrate[z_ + y_, 
                k___] :> 
              Integrate[z, k] + Integrate[y, k] // 
At[2, 
            Map[Apart, #] &]\) // At[2, {3, 4}, Simplify]\) // 
      At[2, \(-{1, 2}\), Simplify]\)

Frequentemente, as condições iniciais indicam que a corda deformada parte do repouso, i.e. ("\[PartialD]u"/∂)t[x, 0] == 0. 
Neste caso g[x] = 0 e u[x, 
          t] == (1/2) (f[x + v t] + f[x - v t]). 
No caso contrário a corda parte do equilíbrio u[x, 0] = 0 com velocidade inicial dada por g[x], e então u[x, 
          t] == (1/2) (G[x + v t] - G[x - v t]), onde G[x]==(1/v)(∫_ξ^x)g[η] dη.

Para além de condições iniciais pede-se ainda que se verifiquem condições fronteira, i.e. no caso u[x, 0] == 0 e u[x, L] == 0.  Vamos considerar o caso de g[x] == 0.

u[x_, t_] := (1/2) (f[x + v t] + f[x - v t])

O método das imagens permite visualizar de forma simples a solução pretendida. Consideremos uma solução que é a sobreposição de três destas formas de onda: para 0 < x < L, (Subscript[u, o])[x, 0] == f[x], onde f[x ≤ 0] == 
        f[x ≥ L] == 0 .  Esta forma propaga-se para a esquerda como f[x + v t] e para a direita como f[x - v t], e para garantir que u[0, t] == u[L, t] == 0 é necessário adicionar à esquerda uma solução partindo de (Subscript[u, -1])[x, 0] == -f[-x] para -L < x < 0,  e à direita (Subscript[u, +1])[x, 0] == -f[2  L - x]para L < x < 2  L,  . Estas condições iniciais para segmentos virtuais da corda são imagens virtuais, com focos em x == 0 e x == L,  da forma inicial da corda no intervalo [0, L].

SetOptions[Plot, Frame -> True, Axes -> False, 
        FrameTicks -> None, PlotPoints -> 60];

Vamos  exemplificar com L = 4 e v = 0.3.

L = 4;
v =  .3;

Clear[f]
f[x_] := 
      If[Abs[x - (L/4)] < (L/8), (1/2) (Cos[(x - (L/2)) 2  π] + 1), 
        0]
f[x_ /; x < 0 && x ≥ -L] := -f[-x]
f[x_ /; x > L && x ≤ 2  L] := -f[2  L - x]

Plot[f[x], {x, -4 L, 4  L}, PlotRange -> {All, {-1, 1}}, 
        FrameTicks -> {{{0, 
                0}, {L, "L"}, {-L, "-L"}, {2  L, "2L"}}, 
            None}, AspectRatio -> (2/(8  L)), 
        PlotPoints -> 8 L (2^10)];

Infelizmente esta não é a solução final porque, embora se verifiquem as condições fronteira, apenas se observam duas reflexões.   Para se observar um movimento periódico tem que se  estender a definição de f da seguinte forma

Clear[f]
f[x_] := f[x] = 
        If[Abs[x - (L/4)] < (L/8), 
          (1/2) (Cos[(x - (L/2)) 2  π] + 1), 0]
f[x_ /; x < 0] := f[x] = -f[-x]
f[x_ /; x > L] := f[x] = -f[2  L - x]

Obtemos assim uma condição inicial periódica de período 2  L sobre toda a recta ℝ, e s sua evolução temporal como solução da equação da corda vibrante dá-nos a solução pretendida no intervalo [0, L]

Plot[f[x], {x, -6 L, 6  L}, PlotRange -> {All, {-1, 1}}, 
        FrameTicks -> {{{0, 
                0}, {L, "L"}, {-L, "-L"}, {2  L, "2L"}}, 
            None}, AspectRatio -> (2/(12  L)), 
        PlotPoints -> 8 L (2^10)];

A seguinte animação mostra como u[0, t] == 
        (Subscript[u, -1])[0, t] + 
          (Subscript[u, o])[0, t] ==  -(f[v t] + f[-v t]) + (f[
                v t] + f[- v t]) == 0 e u[L, t] == 
        (Subscript[u, o])[L, t] + 
          (Subscript[u, +1])[L, 
            t] == (f[L + v t] + 
              f[L - v t]) - (f[2  L  - L + v t] + 
              f[2  L - L - v t]) == 0.

Table[
              Plot[u[x, t], {x, -L, 2  L}, 
                PlotRange -> {{- .01, L +  .01}, {-1.2, 1.2}}, 
                AspectRatio -> (2.4/L), 
                DisplayFunction -> Identity], {t, 
                0, ((2  L)/v),  .03 (( 2  L)/v)}] // Partition[#, 6] & // 
          GraphicsArray[#, GraphicsSpacing -> 0.1] & // 
        SH[#, Dspl] &;

Module[{nb = EvaluationNotebook[], Target, tar}, 
      Target = ToString[tar] <> 
          ToString[$SessionID]; 
CellPrint[
        Cell["", "Graphics", 
          CellTags -> Target]]; 
Do[
        Plot[u[x, t], {x, -L, 2  L}, 
          PlotRange -> {{- .01, L +  .01}, {-2, 2}}, 
          FrameTicks -> {{{0, 0}, {L, "L"}}, None}, 
          AspectRatio -> (4/L), 
          DisplayFunction -> (NotebookFind[nb, Target, All, 
                CellTags]; 
NotebookWrite[nb, 
                  Cell[GraphicsData["PostScript", 
                      DisplayString[#]], "Graphics", 
                    CellTags -> Target]] &)], {t, 
          0, ((2  L)/v),  .03 (( 2  L)/v)}]]

Olhemos agora para o movimento de um ponto fixo da corda. O som gerado pela corda tem mais a ver com este movimento temporal do que com a forma espacial da corda (o som corresponde a uma oscilação longitudinal do ar, enquanto a deformação da corda é transversal).  Para ouvir um som tem que existir uma oscilação de frequência 20 Hz <  ν  <  15 kHz.  Quando v é conhecido, podemos estimar que ν ≈ (v/L)para formas iniciais muito simples.

L = 4;
v = 831. ;

Clear[f]
f[x_] := (-UnitStep[x - (L/2)] + UnitStep[x]) 
        (1/2) Cos[3 (x - (L/4)) ((2  π)/L)]
f[x_ /; x < 0] := -f[-x]
f[x_ /; x > L] := -f[
          Mod[2  L (IntegerPart[((2  L + x)/(2  L))]) - x, 
            2  L]]

Plot[f[x], {x, -3 L, 5  L}, PlotRange -> {All, {-1, 1}}, 
        FrameTicks -> {{{0, 
                0}, {L, "L"}, {-L, "-L"}, {2  L, "2L"}}, 
            None}, AspectRatio -> (2/(7  L)), PlotPoints -> 1000, 
        PlotRegion -> {{-1.5, 1.5}, {-1.3, 1.3}}];

u[x_, t_] := (1/2) (f[x + v t] + f[x - v t])

Play[u[(L/2), t], {t, 0, 1}, 
        PlayRange ->  .4 {-1, 1}];

Compare este com um tom puro de 220 Hz:

Play[ Sin[2  π 220.  t], {t, 0, 1}];

Problema 4

A evolução de uma partícula num potencial V[x] é descrita em mecânica quântica por uma equação do tipo i ((∂ψ)/(∂t))==(-(1/2) (((∂^2))/(∂(x^2)))+V)ψ ≡ ℍψ. Usando transformadas de Fourier mostre que, para a partícula livre (V = 0) esta equação é fácilmente resolúvel e construa o propagador \(TraditionalForm\`\(e\^\(\(-iℍΔ\) \[RawEscape]t\)\) ψ[x, 
            t] -> ψ[x, t + Δt]\) usando uma discretização adequada para uma forma de onda inicial arbitrária ψ[x, 0]. Represente a sua evolução temporal gráficamente.