Computational Physics in Julia 2024-12-09 PDF

Computational Physics in Julia Noel Araujo Moreira 2024-12-09 Table of contents Preface 3 I Conceitos Básicos 4...

Computational Physics in Julia Noel Araujo Moreira 2024-12-09 Table of contents Preface 3 I Conceitos Básicos 4 1 Diferenciação Numérica 6 1.1 Diferenças Finitas.................................. 6 1.2 Melhorar a Precisão da Derivada.......................... 9 1.3 Acabou?........................................ 12 2 Os Detalhes da Diferenciação 14 2.1 As Bordas....................................... 14 2.2 O Espaçamento ℎ................................... 15 2.3 Derivada em 2D................................... 17 2 Preface Course under construction 3  Part I Conceitos Básicos 4 A modelagem computacional física depende de três operações essenciais: Diferenciação Quadratura (integração definida) Busca de raízes Essas três operações numéricas fundamentais são essenciais porque permitem resolver numeri- camente problemas físicos quando soluções analíticas não existem. A diferenciação numérica calcula taxas de variação e gradientes em sistemas onde a função é conhecida apenas em pontos discretos. A quadratura determina áreas, volumes, trabalho e outras quantidades físicas que requerem somas contínuas. A busca de raízes localiza pontos de equilíbrio, zeros de funções e soluções de equações transcendentais. 5 1 Diferenciação Numérica A diferenciação numérica calcula derivadas de uma função 𝑓(𝑥) em um ponto x quando não é possível obter uma expressão analítica. Dois casos requerem diferenciação numérica: 1. Quando 𝑓(𝑥) é calculada através de um procedimento numérico que não permite derivação direta Exemplo: Função implícita definida por equação transcendental 𝑥𝑒𝑥 = 𝑐𝑜𝑠 (𝑥) 2. Quando 𝑓(𝑥) é conhecida apenas em pontos discretos 𝑥1.𝑥2 ,..., 𝑥𝑛 Exemplo: Medições de pressão atmosférica em diferentes altitudes 1.1 Diferenças Finitas Para tratar estes casos numericamente, consideramos uma malha de pontos igualmente espaça- dos por um intervalo ℎ - também conhecido como step size. Figure 1.1: Supomos uma função qualquer 𝑓(𝑥), aonde você conhece os valores de 𝑓(𝑥 ± ℎ). Começamos expandindo 𝑓(𝑥 + ℎ) e 𝑓(𝑥 − ℎ) em série de Taylor em torno de x: 6  ℎ2 ″ ℎ3 ℎ4 𝑓(𝑥 + ℎ) = 𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 ‴ (𝑥) + 𝑓 (𝑖𝑣) (𝑥) + 𝒪(ℎ5 ) 2! 3! 4! ℎ2 ″ ℎ3 ℎ4 𝑓(𝑥 − ℎ) = 𝑓(𝑥) − ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) − 𝑓 ‴ (𝑥) + 𝑓 (𝑖𝑣) (𝑥) + 𝒪(ℎ5 ) 2! 3! 4! Subtraindo estas equações: ℎ3 ‴ 𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ) = 2ℎ𝑓 ′ (𝑥) + 2 𝑓 (𝑥) + 𝒪(ℎ5 ) 3! Dividindo por 2ℎ: 𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ) ℎ2 = 𝑓 ′ (𝑥) + 𝑓 ‴ (𝑥) + 𝒪(ℎ4 ) 2ℎ 6 Rearranjando: 𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ) ℎ2 ‴ 𝑓 ′ (𝑥) = − 𝑓 (𝑥) + 𝒪(ℎ4 ) 2ℎ 6 tem um termo de erro proporcional a ℎ2 que depende da terceira derivada 𝑓 ‴ (𝑥). Este termo de erro pode ser minimizado por duas razões: 1. O fator ℎ2 é pequeno quando escolhemos um passo h suficientemente pequeno 2. Para muitas funções físicas relevantes, a terceira derivada 𝑓 ‴ (𝑥) é naturalmente pequena em magnitude quando comparada com 𝑓 ′ (𝑥) Por exemplo, para funções polinomiais de grau 2 ou menor, 𝑓 ‴ (𝑥) = 0 e a fórmula se torna exata. Para funções suaves como 𝑠𝑖𝑛(𝑥), a terceira derivada (−𝑐𝑜𝑠(𝑥)) tem magnitude similar à função original, então o erro é efetivamente da ordem de ℎ2 - por favor, não ignore essa informação. Finalmente podemos escrever A fórmula de diferença central de 2 pontos: 𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ) 𝑓 ′ (𝑥) ≈ 2ℎ Vejamos um exemplo de código concreto em Julia. using PrettyTables 7  function numerical_derivative(f, x; h=1e-6) return (f(x + h) - f(x - h)) / (2h) end x = 1.0 exact = cos(x) # Create data arrays h_values = Float64[] errors = Float64[] # Calculate values for i in 1:10 h = 10.0^(-i) fprime = numerical_derivative(sin, x, h=h) diff = abs(exact - fprime) push!(h_values, h) push!(errors, diff) end # Create data matrix and headers data = hcat(h_values, errors) headers = ["Step Size (h)", "Erro"] # Create formatter for scientific notation fmt = ft_printf("%0.2e") # Print table with proper formatter specification pretty_table(data, header = headers, formatters = (fmt, fmt)) Cujos valores são: Table 1.1: Erro relativo para fórmula de 2 pontos Step Size (h) Erro 1.00e-01 9.00e-04 1.00e-02 9.00e-06 1.00e-03 9.01e-08 1.00e-04 9.00e-10 8  Step Size (h) Erro 1.00e-05 1.11e-11 1.00e-06 2.77e-11 1.00e-07 1.94e-10 1.00e-08 2.58e-09 1.00e-09 2.97e-09 1.00e-10 5.85e-08 Vamos analisar os resultados. Inicialmente, para valores de h entre 10−1 e 10−4 , observamos uma melhoria consistente na precisão, com o erro diminuindo. Entretanto, quando h se torna muito pequeno (ℎ < 10−7 ), o erro paradoxalmente começa a aumentar. A explicação reside na aritmética de precisão finita do computador: quando h é muito pequeno, os valores de 𝑓(𝑥 + ℎ) e 𝑓(𝑥 − ℎ) tornam-se numericamente muito próximos. Isso é consequência de como números são armazenados pelo computador, e em cursos como Arquitetura de Computadores ou Cálculo Numérico você iria estudar isso em muito com mais rigor. A mensagem que você deve lembrar é: não basta diminuir ℎ, existe um valor ótimo. 1.2 Melhorar a Precisão da Derivada Lembra da derivação da fórmula de 2 pontos, nós desprezamos calculamos a derivada usando 2 pontos. Vamos melhor nossa derivada usando mais pontos. Para entender por que isso funciona, considere que ao usar mais termos na expansão de Taylor, aproximamos melhor a função original por um polinômio de maior grau. Enquanto a fórmula de 2 pontos assume que 𝑓(𝑥) é bem aproximada por um polinômio quadrático no intervalo [−ℎ, +ℎ], uma fórmula de 5 pontos por exemplo, usaria um polinômio de quarto grau no intervalo [−2ℎ, 2ℎ]. Para derivar a fórmula de 5 pontos, expandimos 𝑓(𝑥 ± ℎ) e 𝑓(𝑥 ± 2ℎ) em série de Taylor: ℎ2 ″ ℎ3 ℎ4 𝑓(𝑥 ± ℎ) = 𝑓(𝑥) ± ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) ± 𝑓 ‴ (𝑥) + 𝑓 (𝑖𝑣) (𝑥) + 𝒪(ℎ5 ) 2! 3! 4! 4ℎ2 ″ 8ℎ3 ‴ 16ℎ4 (𝑖𝑣) 𝑓(𝑥 ± 2ℎ) = 𝑓(𝑥) ± 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) ± 𝑓 (𝑥) + 𝑓 (𝑥) + 𝒪(ℎ5 ) 2! 3! 4! 9  Figure 1.2: Fazemos o mesmo processo anterior, mas com mais pontos. Combinando estas equações com coeficientes apropriados para cancelar os termos até ordem h : 1 𝑓 ′ (𝑥) = [−𝑓(𝑥 + 2ℎ) + 8𝑓(𝑥 + ℎ) − 8𝑓(𝑥 − ℎ) + 𝑓(𝑥 − 2ℎ)] + 𝒪(ℎ4 ) 12ℎ Esta fórmula tem erro proporcional a ℎ4 , duas ordens melhor que a fórmula de 2 pontos, que tem erro proporcional a ℎ2. Ĺ Derivação da fórmula de 5 pontos Para derivar a fórmula de 5 pontos, vamos expandir cada termo em série de Taylor em torno de x: Para 𝑓(𝑥 ± ℎ): ℎ2 ″ ℎ3 ℎ4 𝑓(𝑥 ± ℎ) = 𝑓(𝑥) ± ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) ± 𝑓 ‴ (𝑥) + 𝑓 (𝑖𝑣) (𝑥) + 𝒪(ℎ5 ) 2! 3! 4! Para 𝑓(𝑥 ± 2ℎ): 4ℎ2 ″ 8ℎ3 ‴ 16ℎ4 (𝑖𝑣) 𝑓(𝑥 ± 2ℎ) = 𝑓(𝑥) ± 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) ± 𝑓 (𝑥) + 𝑓 (𝑥) + 𝒪(ℎ5 ) 2! 3! 4! Será muuuuuuuuito conveniente afirmar que: 𝑓 ′ (𝑥) = −𝑓(𝑥 + 2ℎ) + 8𝑓(𝑥 + ℎ) − 8𝑓(𝑥 − ℎ) + 𝑓(𝑥 − 2ℎ) 10  Agora, com paciência vamos calcular a seguinte expressão 1 [−𝑓(𝑥 + 2ℎ) + 8𝑓(𝑥 + ℎ) − 8𝑓(𝑥 − ℎ) + 𝑓(𝑥 − 2ℎ)] 12ℎ Substituindo as expansões: 1 4ℎ3 ‴ 2ℎ4 (𝑖𝑣) [−{𝑓(𝑥) + 2ℎ𝑓 ′ (𝑥) + 2ℎ2 𝑓 ″ (𝑥) + 𝑓 (𝑥) + 𝑓 (𝑥)} 12ℎ 3 3 ℎ2 ″ ℎ3 ℎ4 +8{𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 ‴ (𝑥) + 𝑓 (𝑖𝑣) (𝑥)} 2 6 24 ℎ2 ″ ℎ3 ℎ4 −8{𝑓(𝑥) − ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) − 𝑓 ‴ (𝑥) + 𝑓 (𝑖𝑣) (𝑥)} 2 6 24 4ℎ3 ‴ 2ℎ4 (𝑖𝑣) +{𝑓(𝑥) − 2ℎ𝑓 ′ (𝑥) + 2ℎ2 𝑓 ″ (𝑥) − 𝑓 (𝑥) + 𝑓 (𝑥)}] 3 3 Agrupando termos: - Termos em 𝑓(𝑥) se cancelam - Termos em 𝑓 ′ (𝑥): 12ℎ𝑓 ′ (𝑥) - Termos em 𝑓 ″ (𝑥) se cancelam - Termos em 𝑓 ‴ (𝑥) se cancelam - Termos em 𝑓 ( 𝑖𝑣)(𝑥) contribuem para 𝒪(ℎ4 ) Portanto: 1 [12ℎ𝑓 ′ (𝑥) + 𝒪(ℎ5 )] = 𝑓 ′ (𝑥) + 𝒪(ℎ4 ) 12ℎ Alterando o código anterior para criar tabela de erros function numerical_derivative(f, x; h=1e-6) # 5-point formula: f'(x) [-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)]/(12h) return (-f(x + 2h) + 8f(x + h) - 8f(x - h) + f(x - 2h)) / (12h) end Temos a nova tabela: Table 1.2: Erro relativo para fórmula de 5 pontos Step Size (h) Error 1.00e-01 1.80e-06 1.00e-02 1.80e-10 1.00e-03 7.56e-14 1.00e-04 4.24e-13 1.00e-05 3.87e-14 1.00e-06 9.21e-12 11  Step Size (h) Error 1.00e-07 2.87e-10 1.00e-08 2.58e-09 1.00e-09 2.48e-08 1.00e-10 5.85e-08 Os novos resutlados são promissores, para ℎ = 10−5 , o erro é ℎ = 10−14. Nesni assim, dimin- uindo o valor de h, não diminui o erro da derivada. 1.3 Acabou? A escolha da fórmula de diferenciação numérica depende da precisão necessária para seu prob- lema. A seguir você pode consultar algumas fórmulas já conhecidas na literatura. Fórmulas de 2 ou 3 pontos têm erro proporcional a 𝒪(ℎ2 ), fórmulas de 5 pontos têm erro proporcional a 𝒪(ℎ4 ), e existem outras fórmulas de ordem ainda maior. Table 1.3: Derivadas de Primeira Ordem Fórmula Erro 𝑓(𝑥+ℎ)−𝑓(ℎ) Diferença para 𝑓 ′ (𝑥) ≈ ℎ 𝒪(ℎ) frente de 2 pontos 𝑓(𝑥)−𝑓(𝑥−ℎ) Diferença para 𝑓 ′ (𝑥) ≈ ℎ 𝒪(ℎ) trás de 2 pontos 𝑓(𝑥+ℎ)−𝑓(𝑥−ℎ) Diferença central 𝑓 ′ (𝑥) ≈ 2ℎ 𝒪(ℎ2 ) de 2 pontos 𝑓(𝑥−2ℎ)−8𝑓(𝑥−ℎ)+8𝑓(𝑥+ℎ)−𝑓(𝑥+2ℎ) Diferença central 𝑓 ′ (𝑥) ≈ 12ℎ 𝒪(ℎ4 ) de 5 pontos Table 1.4: Derivadas de Segunda e Terceira Ordem Fórmula Erro 𝑓(𝑥−ℎ)−2𝑓(𝑥)+𝑓(𝑥+ℎ) Diferença central 𝑓 ″ (𝑥) ≈ ℎ2 𝒪(ℎ2 ) de 3 pontos (2ª derivada) 12  Fórmula Erro −𝑓(𝑥−2ℎ)+16𝑓(𝑥−ℎ)−30𝑓(𝑥)+16𝑓(𝑥+ℎ)−𝑓(𝑥+2ℎ) Diferença central 𝑓 ″ (𝑥) ≈ 12ℎ2 𝒪(ℎ3 ) de 5 pontos (2ª derivada) −𝑓(𝑥)+3𝑓(𝑥+ℎ)−3𝑓(𝑥+2ℎ)+𝑓(𝑥+3ℎ) Diferença para 𝑓 ‴ (𝑥) ≈ ℎ3 𝒪(ℎ) frente de 4 pontos (3ª derivada) −𝑓(𝑥−2ℎ)+2𝑓(𝑥−ℎ)−2𝑓(𝑥+ℎ)+𝑓(𝑥+2ℎ) Diferença central 𝑓 ‴ (𝑥) ≈ 2ℎ3 𝒪(ℎ2 ) de 5 pontos (3ª derivada) Baseado no problema que você está trabalhando, você será o responsável por decidir que fórmula usar. Felizmente, as fórmulas mais simples com 3 pontos, normalmente são o suficiente para vários problemas clássicos da física. 13 2 Os Detalhes da Diferenciação 2.1 As Bordas Ao calcular derivadas numéricas com diferenças finitas, você enfrentará um problema nas extremidades do seu conjunto de dados. A fórmula da derivada central 𝑓(𝑥+ℎ)−𝑓(𝑥−ℎ) 2ℎ não funciona para o primeiro e o último ponto, pois não há dados disponíveis antes do primeiro ponto ou após o último. Uma possível solução é usar métodos diferentes para as bordas. 𝑥𝑖+1 −𝑥𝑖−1 Para os pontos internos, aplique a fórmula de diferença central 𝑡𝑖+1 −𝑡𝑖−1. Para o primeiro ponto, usa a diferença para frente 𝑥𝑡2 −𝑥 1. 2 −𝑡1 𝑥𝑛 −𝑥𝑛−1 Para o último ponto, empregue a diferença para trás 𝑡 −𝑡. 𝑛 𝑛−1 Vamos vizualizar: using CairoMakie # Dados de exemplo t = [0.0, 1.0, 2.0, 3.0, 4.0] x = [0.0, 2.0, 4.5, 8.0, 12.5] # Função para calcular a derivada numérica function numerical_derivative(t, x) n = length(t) dx = similar(x) # Diferença finita central para pontos internos for i in 2:n-1 dx[i] = (x[i+1] - x[i-1]) / (t[i+1] - t[i-1]) end # Diferença finita para as bordas dx = (x - x) / (t - t) 14  dx[n] = (x[n] - x[n-1]) / (t[n] - t[n-1]) return dx end # Calcular a derivada dx = numerical_derivative(t, x) # Criar o plot fig = Figure(size=(800, 600), fontsize=14) ax1 = Axis(fig[1, 1], xlabel="Tempo", ylabel="Posição") scatter!(ax1, t, x, color=:royalblue, markersize=15, label="Dados Originais") lines!(ax1, t, x, color=:royalblue, linewidth=2) axislegend(ax1, position=:lt, framevisible=false) ax2 = Axis(fig[2, 1], xlabel="Tempo", ylabel="Velocidade") scatter!(ax2, t, dx, color=:firebrick, markersize=15, label="Derivada Numérica") lines!(ax2, t, dx, color=:firebrick, linewidth=2) axislegend(ax2, position=:lt, framevisible=false) # Ajustar espaçamento entre subplots fig[1, 1] = ax1 fig[2, 1] = ax2 rowgap!(fig.layout, 20) # Exibir a figura display(fig) 2.2 O Espaçamento ℎ Em conjuntos de dados reais, onde o intervalo entre medições pode variar, e ao calcular a derivada numérica usando diferenças finitas dados discretos, o valor de ℎ não é fixo como na fórmula teórica, mas varia de acordo com o espaçamento dos pontos. Para um ponto 𝑖, ℎ é metade da distância entre os pontos adjacentes utilizados no cálculo. ℎ = (𝑡𝑖+1 − 𝑡𝑖−1 )/2 Repare que no código anterior, usamos 2ℎ = 𝑡𝑖+1 − 𝑡𝑖−1. 15 Figure 2.1: Posição e Velocidade ao longo do tempo 16 2.3 Derivada em 2D A derivada numérica em 2D, também conhecida como gradiente numérico, é uma extensão do conceito de diferenças finitas para funções de duas variáveis. Para uma função 𝑓(𝑥, 𝑦), calculamos as derivadas parciais em relação a 𝑥 e 𝑦 separadamente. 𝜕𝑓 𝑓(𝑥+ℎ,𝑦)−𝑓(𝑥−ℎ,𝑦) 𝜕𝑓 Utilizamos a fórmula de diferença central para cada direção: 𝜕𝑥 ≈ 2ℎ $ e 𝜕𝑦 ≈ 𝑓(𝑥,𝑦+ℎ)−𝑓(𝑥,𝑦−ℎ) 2ℎ , onde ℎ é um pequeno incremento. Aqui está um exemplo de implementação em Julia que calcula o gradiente numérico de uma função de duas variáveis em um ponto específico function gradiente_numerico_2d(f, x, y; h=1e-5) fx = (f(x+h, y) - f(x-h, y)) / (2h) fy = (f(x, y+h) - f(x, y-h)) / (2h) return [fx, fy] end # Exemplo de uso f(x, y) = x^2 + y^2 # Função de exemplo x, y = 1.0, 2.0 # Quero o gradiente nessas coordenadas grad = gradiente_numerico_2d(f, x, y) println("Gradiente em (1,2): ", grad) 17