Oscilações Acopladas/Problema de Fermi-Pasta-Ulam-Tsingou: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Linha 36: Linha 36:
Uma forma direta de se montar as equações de movimento de um sistema acoplado é pensar em termos dos deslocamentos em relação às posições de equilíbrio (<math>x_{i}</math>, com <math>i=1,2</math>, no sistema considerado aqui). Por exemplo, na primeira das equações acima, a partícula 1 está sujeita à força elástica da mola conectada à parede (termo <math>-kx_{1}</math>) e à força da mola conectada à partícula 2 (termo <math>-k(x_{1}-x_{2}</math>). Esse último termo é definido conforme o seguinte: caso a mola que está conectada às partículas 1 e 2 esteja comprimida, <math>x_{1}</math> deve ser maior que <math>x_{2}</math> (definindo os <math>x_{i}</math> como deslocamentos em relação à posição de equilíbrio positivos para a direita e negativos para a esquerda) e a partícula 1 estará sofrendo uma força que é contrária à compressão, ou seja, com sinal negativo. De modo semelhante para a partícula 2, mas nesse caso, devido à posição ocupada pela partícula 2 na cadeia, com <math>x_{1}</math> maior que <math>x_{2}</math>, o sentido da força contrária à compressão tem sinal positivo.
Uma forma direta de se montar as equações de movimento de um sistema acoplado é pensar em termos dos deslocamentos em relação às posições de equilíbrio (<math>x_{i}</math>, com <math>i=1,2</math>, no sistema considerado aqui). Por exemplo, na primeira das equações acima, a partícula 1 está sujeita à força elástica da mola conectada à parede (termo <math>-kx_{1}</math>) e à força da mola conectada à partícula 2 (termo <math>-k(x_{1}-x_{2}</math>). Esse último termo é definido conforme o seguinte: caso a mola que está conectada às partículas 1 e 2 esteja comprimida, <math>x_{1}</math> deve ser maior que <math>x_{2}</math> (definindo os <math>x_{i}</math> como deslocamentos em relação à posição de equilíbrio positivos para a direita e negativos para a esquerda) e a partícula 1 estará sofrendo uma força que é contrária à compressão, ou seja, com sinal negativo. De modo semelhante para a partícula 2, mas nesse caso, devido à posição ocupada pela partícula 2 na cadeia, com <math>x_{1}</math> maior que <math>x_{2}</math>, o sentido da força contrária à compressão tem sinal positivo.


Uma dificuldade imposta para a resolução do sistema (1) é o fato das equações serem acopladas: note-se que a aceleração da partícula 1 depende da posição da partícula 2, e vice-versa. Vamos supor que esse sistema de equações tenha soluções nas formas (com <math>B_{1}</math> e <math>B_{2}</math> constantes):
Uma dificuldade imposta para a resolução do sistema (1) é o fato das equações serem acopladas: note-se que a aceleração da partícula 1 depende da posição da partícula 2, e vice-versa. Vamos supor que esse sistema de equações tenha soluções nas formas (com <math>A_{1}</math> e <math>A_{2}</math> constantes):


:<math>\begin{align}
:<math>\begin{align}
x_{1}(t) &= B_{1}e^{-i\omega t} \\
x_{1}(t) &= A_{1}e^{-i\omega t} \\
x_{2}(t) &= B_{2}e^{-i\omega t} \quad (2)
x_{2}(t) &= A_{2}e^{-i\omega t} \quad (2)
\end{align}
\end{align}
</math>
</math>
Linha 47: Linha 47:


:<math>\begin{align}
:<math>\begin{align}
-mw^{2}B_{1}e^{-i\omega t} + 2kB_{1}e^{-i\omega t} -kB_{1}e^{-i\omega t} &= 0 \\
-mw^{2}A_{1}e^{-i\omega t} + 2kA_{1}e^{-i\omega t} -kA_{1}e^{-i\omega t} &= 0 \\
-mw^{2}B_{2}e^{-i\omega t} + 2kB_{2}e^{-i\omega t} -kB_{2}e^{-i\omega t} &= 0 \\
-mw^{2}A_{2}e^{-i\omega t} + 2kA_{2}e^{-i\omega t} -kA_{2}e^{-i\omega t} &= 0 \\
\end{align}</math>
\end{align}</math>


Linha 54: Linha 54:


:<math>\begin{align}
:<math>\begin{align}
(2k-mw^{2})&B_{1} &-kB_{2} &= 0 \\
(2k-mw^{2})&A_{1} &-kA_{2} &= 0 \\
-k&B_{1} &+(2k-mw^{2})B_{2} &= 0 \quad (3)\\  
-k&A_{1} &+(2k-mw^{2})A_{2} &= 0 \quad (3)\\  
\end{align}</math>
\end{align}</math>


O sistema de equações (3) apenas terá soluções não triviais se o determinante dos coeficientes dos <math>B_{i}</math> for igual a zero, i.e.:
O sistema de equações (3) apenas terá soluções não triviais se o determinante dos coeficientes dos <math>A_{i}</math> for igual a zero, i.e.:


:<math>\begin{vmatrix}
:<math>\begin{vmatrix}
Linha 68: Linha 68:


:<math>\begin{align}
:<math>\begin{align}
\omega_{1} = \sqrt{\frac{k}{m}} \quad \quad \omega_{2} = \sqrt{\frac{3k}{m}}
\omega_{1} = \sqrt{\frac{k}{m}} \quad \quad \omega_{2} = \sqrt{\frac{3k}{m}} \quad (4)
\end{align} </math>
\end{align} </math>


<math>\omega_{1}</math> e <math>\omega_{2}</math> são as frequências características ou autofrequências do sistema.
<math>\omega_{1}</math> e <math>\omega_{2}</math> são as frequências características ou autofrequências do sistema. As soluções mais gerais do sistema de equações diferenciais lineares (1) vão ser então combinações lineares das soluções (2) com as frequências dadas por (4):
 
:<math>\begin{align}
x_{1}(t) = a_{1}^{1}e^{-i\omega_{1} t} + a_{1}^{2}e^{-i\omega_{2} t} &= \sum_{k}a_{1}^{k}e^{-i\omega_{k}t}\\
x_{2}(t) = a_{2}^{1}e^{-i\omega_{1} t} + a_{2}^{2}e^{-i\omega_{2} t} &= \sum_{k}a_{2}^{k}e^{-i\omega_{k}t}\quad (5)
\end{align}
</math>


== Adição de Termos Não-Lineares: Problema de FPUT ==  
== Adição de Termos Não-Lineares: Problema de FPUT ==  

Edição das 17h53min de 30 de abril de 2022

Grupo: Paula Pandolfo, Ramiro de Souza, Samuel Dieterich e Wallace Carvalho

Objetivo: Este trabalho tem dois objetivos principais: apresentar alguns resultados analíticos de osciladores lineares acoplados, comparando esses resultados com simulações computacionais; e implementar o modelo de osciladores acoplados com a adição de um termo quadrático, conforme inicialmente apresentado pelo artigo original do problema de Fermi-Pasta-Ulam-Tsingou (FPUT), analisando os resultados. Apresentaremos algumas simulações dos casos bidimensionais, mas as análises de resultados serão restritas aos casos unidimensionais, por simplicidade. Inicialmente será introduzido o formalismo de oscilações acopladas lineares. [falta complementar]

Introdução

Os osciladores são talvez os sistemas mais estudados na Física, sendo capazes de modelar uma ampla gama de fenômenos, como, p. ex., pêndulos, circuitos eletrônicos, interações moleculares. O comportamento linear desses sistemas, em particular, possui resultados analíticos bem conhecidos.

O problema de FPUT (Enrico Fermi, John R. Pasta, Stanislaw M. Ulam, Mary Tsingou) resulta da análise computacional de um sistema de partículas que apenas interagem com seus vizinhos, com interações modeladas por oscilações acopladas com a adição de um termo não-linear, que pode ser quadrático ou cúbico. O intuito original da simulação era estudar como esse sistema evolui para o equilíbrio térmico. Se as forças fossem estritamente lineares, a energia alocada em cada modo de vibração não se distribuiria entre os demais modos, ou seja, não se atingiria o equilíbrio térmico. Entretanto, com a adição dos termos não lineares, pelo Teorema da Equipartição da Energia, supunha-se que, após um certo tempo, a energia total do sistema seria distribuída uniformemente entre os modos normais de vibração, o que significaria que o sistema teria atingido o equilíbrio térmico. Entretanto, isso não foi observado.

O caso foi estudado pela primeira vez em Los Alamos, nos Estados Unidos, e implementado no computador MANIAC I (Mathematical Analyzer Numerical Integrator and Automatic Computer Model I). Além dos três participantes coautores do artigo que relatou o caso em 1955, Mary Tsingou implementou o código e resolveu numericamente o sistema. Atualmente, por essa razão, o paradoxo é denominado pela sigla FPUT (Fermi-Pasta-Ulam-Tsingou).

A abordagem adotada no presente trabalho é a seguinte: inicialmente, serão apresentados alguns resultados teóricos bem conhecidos de osciladores lineares acoplados. A seguir, compararemos esses resultados com simulações computacionais. [falta complementar]

Osciladores Lineares Acoplados

Um modelo geral de sistema unidimensional de osciladores lineares acoplados é ilustrado pela Figura 1. Para fins de simplificação do problema, estamos considerando que todas as massas e constantes das molas são iguais, mas esse não precisaria ser o caso.

  • Figura 1. Ilustração de um modelo geral unidimensional de oscilações acopladas. A figura foi criada por Paula Pandolfo, uma das integrantes do grupo.
  • Cada partícula possui duas vizinhas, com as quais interage por meio das molas, exceto as partículas localizadas nos extremos da cadeia, que possuem apenas uma partícula vizinha cada. As interações das partículas dos extremos das cadeias se restringem, portanto, à interação com uma vizinha e com uma mola conectada a uma das paredes externas à cadeia. A posição de cada partícula pode ser descrita por um grau de liberdade associado ao deslocamento em relação à respectiva posição de equilíbrio. No total, um sistema com partículas terá, portanto, graus de liberdade. Vamos tratar aqui o caso em que as forças das molas são lineares, i.e., dadas por .

    N=2

    Vamos inicialmente considerar o caso do oscilador linear acoplado mais simples, com duas partículas (), cada uma com massa , e três molas com os mesmos valores de constantes, .

    As equações de movimento do sistema são:

    Uma forma direta de se montar as equações de movimento de um sistema acoplado é pensar em termos dos deslocamentos em relação às posições de equilíbrio (, com , no sistema considerado aqui). Por exemplo, na primeira das equações acima, a partícula 1 está sujeita à força elástica da mola conectada à parede (termo ) e à força da mola conectada à partícula 2 (termo ). Esse último termo é definido conforme o seguinte: caso a mola que está conectada às partículas 1 e 2 esteja comprimida, deve ser maior que (definindo os como deslocamentos em relação à posição de equilíbrio positivos para a direita e negativos para a esquerda) e a partícula 1 estará sofrendo uma força que é contrária à compressão, ou seja, com sinal negativo. De modo semelhante para a partícula 2, mas nesse caso, devido à posição ocupada pela partícula 2 na cadeia, com maior que , o sentido da força contrária à compressão tem sinal positivo.

    Uma dificuldade imposta para a resolução do sistema (1) é o fato das equações serem acopladas: note-se que a aceleração da partícula 1 depende da posição da partícula 2, e vice-versa. Vamos supor que esse sistema de equações tenha soluções nas formas (com e constantes):

    Essa suposição é fisicamente justificável: sabemos que as soluções são oscilatórias, e exponenciais imaginárias podem ser escritas em termos de senos e cossenos pela fórmula de Euler. Se substituirmos a equação (2) na equação (1) e rearranjarmos os termos, obtemos:

    Ou, eliminando as exponenciais e reagrupando termos:

    O sistema de equações (3) apenas terá soluções não triviais se o determinante dos coeficientes dos for igual a zero, i.e.:

    As soluções da equação acima podem ser facilmente obtidas e resultam de uma equação quadrática simples, apresentaremos apenas o resultado:

    e são as frequências características ou autofrequências do sistema. As soluções mais gerais do sistema de equações diferenciais lineares (1) vão ser então combinações lineares das soluções (2) com as frequências dadas por (4):

    Adição de Termos Não-Lineares: Problema de FPUT

    Implementação numérica

    Resultados e discussão

    Programa

    Código em C para FPUT com deslocamento vertical.

    // Libraries
    
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <assert.h>
    
    
    // Defines
    
    #define PI 3.14159265358979323846
    
    
    // Prototypes
    
    double *zero_array(int);
    double *initial_condition(double, int);
    double acceleration(double *, double, double, double);
    double velocity(double, double, double, double);
    double position(double, double, double, double);
    
    
    // Constants
    
    const double k      = 0.95;
    const double m      = 1.05;
    const double alpha  = 1.2;
    
    const double dt     = 0.5;
    const double tmax   = 5500.0;
    
    const double length = 20.0;
    const int N         = 32;
    
    
    // Main
    
    int main(){
    
        // Variables
    
        int i;
        double t = 0.0;
        double *pos, *vel, *acc;
        double acc_new = 0.0;
    
    
        // Output
    
        FILE *file;
        // char output[50];
        char filename[] = "../data/output.dat";
    
    
        // Generate inicial conditions for
        pos = initial_condition(length, N); // position
        vel = zero_array(N);                // velocity
        acc = zero_array(N);                // acceleration
    
        // Output file
        // sprintf(output, filename);
        // file = fopen(output, "w");
        file = fopen(filename, "w");
    
    
        // Time loop
        while (t < tmax){
    
            // Velocity verlet
            for (i=1; i<N; i++){
                pos[i] = position(pos[i], vel[i], acc[i], dt);
            }
    
            for (i=1; i<N; i++){
                acc_new = acceleration((pos+i), k, m, alpha);
                vel[i] = velocity(vel[i], acc[i], acc_new, dt);
                acc[i] = acc_new;
            }
    
            // Save to file
            for (i=0; i<=N; i++) fprintf(file, "%f %f\n", length/N*i, pos[i]);
            fprintf(file, "\n\n");
    
            // Update the current time
            t += dt;
        }
    
    
        // Close the output file
        fclose(file);
    
    
        // Free the pointers
        free(pos);
        free(vel);
        free(acc);
    
    
        return 0;
    
    }
    
    
    // Functions
    
    // Generate an array of zeros (double) with size N
    double *zero_array(int N){
        assert(N>0);
        int i;
        double *pointer;
    
        // Calloc initializes every element with 0.0
        pointer = (double *)calloc(N, sizeof(double));
    
        if (pointer != NULL){
            return(pointer);
        }
        else {
            fprintf(stderr,"Error in routine zero_array\n");
            exit(8);
        }
    }
    
    // Generate an array with the initial condition determined by the FPUT original problem, N elements with x values ranged from 0 to length.
    double *initial_condition(double length, int N){
        assert(N>0);
        int i;
        double *pointer;
        double step = length/N;
    
        pointer = (double *)calloc(N, sizeof(double));
    
        if (pointer != NULL){
            for (i=0; i<=N; ++i) *(pointer+i) = sin(PI*(step*i)/length);
            return(pointer);
        }
        else {
            fprintf(stderr,"Error in routine initial_condition\n");
            exit(8);
        }
    }
    
    // Velocity-Verlet method
    
    // Updates the acceleration
    double acceleration(double *pos, double k, double m, double alpha){
        return k/m * (pos[1] + pos[-1] - 2*pos[0]) * (1.0 + alpha * (pos[1] - pos[-1]));
    }
    
    // Updates the velocity
    double velocity(double vel_old, double acc, double acc_new, double dt){
        return vel_old + 0.5*(acc+acc_new)*dt;
    }
    
    // Updates the position
    double position(double pos_old, double vel, double acc, double dt){
        return pos_old + vel*dt + 0.5*acc*dt*dt;
    }
    


    Script do Gnuplot para gerar as imagens.

    # Reset previously sessions
    reset
    
    # Data input
    file = "../data/output.dat"
    
    # Statistics about the data
    stats file name "A"
    
    # Set range of plots
    set xrange [A_min_x:A_max_x]    # x-axis
    set yrange [-1.2:1.2]           # y-axis
    
    # Customize the output 
    set term pngcairo size 800, 600 enhanced font "Arial,10" fontscale 1.0
    
    # Title and labels
    set title "FPUT - Descolamento vertical" font "Oswald,18"
    set xlabel "Posição (x)" font ",12"
    set ylabel "Altura (y)" font ",12"
    
    # Total number of frames
    tmax = int(A_blocks-1)
    
    # Text position in y
    text_y = (A_max_x-A_min_x)*0.05
    
    # Loop
    do for [i=0:tmax] {  
        # Print: current_frame total_frame
        print i,tmax
        # Output file format
        outfile = sprintf("../img/".dir."/fput%d.png", i)
        set output outfile
        # Time
        unset label
        set label "Iteração: ".i at text_y, 1 font ",10"
        # Plot
        plot file index i with points pt 7 ps 2 lc rgb 'black' notitle
        }
    

    Script para automatizar a criação dos vídeos.

    #!/bin/bash
    
    # Run all the scripts
    # Compile the C code -> Run the executable -> Run the gnuplot script -> Create the video file from the images on the img folder
    
    name=$(date +%s)
    
    gcc -lm fput.c 
    ./a.out
    mkdir -p ../img/$name
    gnuplot -e "dir='$name'" graph.plt
    mkdir -p video
    ffmpeg -r 120 -i ../img/$name/fput%d.png -vcodec libx264 -crf 25 -pix_fmt yuv420p ../video/output-$name.mp4
    

    Referências


    Bibliografia principal

    • Stephen T. Thornton, Jerry B. Marion, "Classical Dynamics of Particles and Systems". Thomson Learning, Belmont, 2004.
    • Giordano, N.J., Nakanishi, H. "Computational Physics". 2nd Edition. Prentice Hall, Upper Saddle River, 2006.