Métodos Numéricos II

Implementação em C++ dos métodos numéricos cobertos na disciplina de Métodos Numéricos II: derivação numérica, integração numérica, autovalores/autovetores, problemas de valor inicial (PVI), problemas de valor de contorno (PVC) e, como módulo extra, detecção de bordas em imagens..

O projeto é uma única CLI interativa em C++, organizada por módulo, cada um com seu próprio submenu.

Sumário

• Como compilar e executar
• Organização do projeto
• Módulo 1 - Derivação Numérica
• Módulo 2 - Integração Numérica
• Módulo 3 - Autovalores e Autovetores
• Módulo 4 - PVI (Problemas de Valor Inicial)
• Módulo 5 - PVC (Problemas de Valor de Contorno)
• Módulo 6 - Processamento de Imagens (Detecção de Bordas)
• Parser de expressões matemáticas
• Bibliotecas de terceiros

Como compilar e executar

Requisitos: CMake ≥ 3.28 e um compilador C++ com suporte a C++17 (GCC, Clang ou MSVC).

mkdir build && cd build
cmake ..
cmake --build .
./methods          # no Windows: methods.exe

O binário methods abre um menu principal com um submenu por módulo; a navegação é toda por número + Enter, e cada submenu tem a opção 0 - Voltar.

Organização do projeto

main.cpp                  # menu principal, despacha para os 6 módulos
include/                  # headers (.h) organizados por módulo
src/                      # implementação (.cpp), mesma estrutura de include/
third_party/              # dependências vendorizadas (ver seção própria)

Cada módulo segue o mesmo padrão de arquivos:

Arquivo	Papel
menuX.cpp	CLI do módulo: lê parâmetros do usuário, chama os algoritmos, imprime resultado
xUtils.cpp	funções de apoio compartilhadas dentro do módulo (álgebra linear, tipos, solvers)
demais .cpp	um arquivo por família de método numérico

Essa separação existe para isolar a CLI (entrada/saída, validação, formatação) da matemática (que não depende de iostream e pode em tese ser testada isoladamente).

Os métodos que não têm uma entrada obrigatória por enunciado (derivação numérica e EDOs de uso geral) aceitam funções digitadas em tempo real via um parser de expressões (seção Parser de expressões matemáticas); os que dependem de uma tarefa específica (integração multidimensional, autovalores, PVC) usam entradas fixas no código, exatamente como definidas nos enunciados.

---

Módulo 1 - Derivação Numérica

Arquivos: src/numericalDiff/, include/numericalDiff/.

Implementa derivadas de 1ª a 4ª ordem por diferenças finitas, todas parametrizadas por passo h:

• Diferenças finitas simples (finiteDiff.h): fórmulas forward, backward e central para as quatro ordens, mais versões recursivas (forwardRec, backwardRec, centralRec) que calculam a derivada de ordem n aplicando a fórmula de 1ª ordem recursivamente sobre a derivada de ordem n-1 - útil para comparar o custo/precisão de aplicar a definição diretamente versus compor diferenças de ordem inferior.
• Fórmulas de Newton com erro de ordem superior (finiteDiffNewton.h): combinações lineares de mais pontos para cancelar termos do erro de truncamento, ex. derivada 1ª com erro O(h²) usando 3 pontos (forward, central ou backward) e derivada 2ª com erro O(h⁴).
• Expansão em série de Taylor (finiteDiffTaylor.h): mesma ideia das fórmulas de Newton, mas cada função retorna também uma estimativa do termo de erro (struct derivate { derivateValue; error; }), permitindo comparar a ordem de convergência teórica com o erro observado ao variar h.

Na CLI (menuNumericalDiff()), o fluxo é: escolher a função (uma das 4 funções de teste em testFunctions1D, ou digitar uma expressão personalizada), escolher o ponto x, o passo h e a família de método.

Módulo 2 - Integração Numérica

Arquivos: src/numericalInt/, include/numericalInt/.

• Newton-Cotes (newtonCotes.h): fórmulas fechadas (trapézio, Simpson 1/3, Simpson 3/8, grau 4) e abertas (trapézio aberto, Milne, graus 3 e 4). newtonCotesParticionado faz refinamento adaptativo: dobra o número de partições sucessivamente até a diferença entre duas estimativas consecutivas ficar abaixo de uma tolerância.
• Gauss-Legendre (gaussLegendre.h): quadratura com 2, 3 ou 4 pontos sobre [xi, xf], via mudança de variável linear para o intervalo de referência [-1, 1]. Tem versão particionada (mesmo esquema adaptativo do Newton-Cotes) e uma versão 2D (gaussLegendre2D) para integrais duplas em domínios retangulares.
• Quadraturas especiais de Gauss (gaussSpecial.h): Gauss-Hermite (peso e^{-s²}, domínio infinito), Gauss-Laguerre (peso e^{-s}, domínio [0, +∞)) e Gauss-Chebyshev (peso 1/√(1-s²), domínio [-1, 1]) - cada uma com nós e pesos tabelados para n = 2, 3, 4.
• Integração com singularidade (singularIntegral.h): trata integrandos com singularidade nos extremos via mudança de variável exponencial, simples (x(s) = tanh(s)) ou dupla (x(s) = tanh(π/2·sinh(s))), que mapeia o intervalo original para toda a reta e concentra pontos de quadratura perto da singularidade. Na prática o integrando transformado converge exponencialmente rápido para 0 nas caudas, então basta truncar a integral em [-c, c] e aumentar c até convergir (singularSolucao2), aplicando Gauss-Legendre no intervalo truncado.
• Integrais multidimensionais (multiDimIntegral.h): área de superfície 3D (∫∫ √((∂f/∂x)² + (∂f/∂y)² + 1) dA) e volume sob uma superfície, sobre domínio retangular ou elíptico. No domínio elíptico usa-se coordenadas polares elípticas para mapear a elipse em um retângulo antes de aplicar Gauss-Legendre 2D.

A função de referência usada nas tarefas de superfície/volume é F_surf(x,y) = 0.2(x² - y²), com derivadas parciais fechadas em function.cpp.

Módulo 3 - Autovalores e Autovetores

Arquivos: src/eigenvalues/, include/eigenvalues/.

• matrixUtils: base de álgebra linear usada por todo o módulo - produto matriz-matriz e matriz-vetor, transposta, norma, decomposição LU sem pivotamento e resolução de sistemas triangulares/gerais via LU (usada, por exemplo, dentro do método da potência inverso).
• Método da Potência (powerMethod.h): três variantes.
  • Regular: converge para o autovalor de maior valor absoluto.
  • Inverso: aplica o método regular sobre A⁻¹ (via LU, sem inverter a matriz explicitamente) para achar o autovalor de menor valor absoluto.
  • Com deslocamento: aplica o inverso sobre A - μI para achar o autovalor mais próximo de um μ escolhido.
• Householder (householderMethod): reduz uma matriz simétrica a uma forma tridiagonal semelhante por uma sequência de reflexões H_i, preservando os autovalores; o produto acumulado H = H1·H2·... permite recuperar os autovetores da matriz original a partir dos autovetores da tridiagonal.
• Jacobi (jacobiMethod): diagonaliza uma matriz simétrica por varreduras sucessivas, cada uma zerando um par (i,j) fora da diagonal via uma rotação G_ij (jacobiRotation); repete até a soma dos quadrados abaixo da diagonal (sumSquaresBelowDiagonal) ficar abaixo de eps.
• QR (qrMethod): itera A → QR → RQ até a matriz convergir para uma forma (quase) diagonal; a decomposição QR é feita por rotações de Givens (givensRotation), que zeram um elemento subdiagonal por vez.

A CLI permite escolher entre as matrizes das tarefas (taskMatrixA1/A2/A3), casos de teste genéricos (simétrica positiva-definida, autovalores próximos, singular, mal-condicionada) ou inserir uma matriz manualmente, e depois roda o método escolhido mostrando os passos intermediários (matrizes após cada varredura/iteração) quando aplicável.

Módulo 4 - PVI (Problemas de Valor Inicial)

Arquivos: src/pvi/, include/pvi/.

O estado do sistema é genérico (State = std::vector<double>), então os mesmos métodos resolvem tanto uma EDO escalar quanto um sistema (como o par velocidade/posição da tarefa de queda com resistência do ar).

• Euler Explícito (explicitEuler.h): S_{i+1} = S_i + dt·F(S_i, t_i).
• Euler Implícito (implicitEuler.h): S_{i+1} = S_i + dt·F(S_{i+1}, t_{i+1}). Como S_{i+1} aparece nos dois lados, o módulo aceita um ImplicitSolver fornecido por fora - pode ser a solução analítica fechada (caso da EDO de queda livre com arrasto, resolvida em edo.cpp) ou um solver genérico por ponto fixo (makeFixedPointImplicitSolver) para qualquer EDO sem forma fechada conhecida.
• Runge-Kutta 2ª, 3ª e 4ª ordem (rungeKutta.h): estágios intermediários calculados a partir de avaliações de F em pontos entre t_i e t_i + dt.
• Preditor-Corretor / Adams (predictorCorrector.h): métodos multi-passo lineares - AB2, ABM3 e ABM4 - que reaproveitam avaliações de F de passos anteriores para reduzir o número de chamadas à função por passo. Como precisam de vários pontos iniciais que ainda não existem no início da integração, cada um se auto-inicializa com um Runge-Kutta de ordem compatível (RK2 para o AB2, RK3 para o ABM3, RK4 para o ABM4).

A EDO fixa pelas tarefas (taskDerivPVI2) é a queda com resistência do ar proporcional à velocidade: dv/dt = -g - (k/m)v, dy/dt = v, com solução exata fechada usada para calcular erro. Além dela, há um catálogo de EDOs de uso geral (testEdos, em edo.cpp) - crescimento exponencial, decaimento, logística com saturação e um caso rígido (stiff) que serve para observar instabilidade de métodos explícitos quando dt não é pequeno o suficiente - e a opção de digitar uma EDO escalar personalizada pelo parser.

Módulo 5 - PVC (Problemas de Valor de Contorno)

Arquivos: src/pvc/, include/pvc/.

Dois problemas de referência são resolvidos por dois métodos diferentes, para comparação:

• PVC1: y'' - y = 0, y(0) = 0, y(1) = 1 (solução exata fechada disponível em exactSolutionPVC1).
• PVC2: ∂²u/∂x² + ∂²u/∂y² = 4 num domínio quadrado com u = 0 no contorno (equação de Poisson 2D).
• Problema final: u'' + 7u' - u = 2, u(0) = 10, u(2) = 1 (solução exata em exactSolutionFinalODE).

Diferenças Finitas (finiteDiff.h): discretiza as derivadas por diferenças centrais de 2ª ordem nos nós internos de uma malha, montando um sistema linear que é resolvido por eliminação de Gauss com pivotamento (solveGauss, PVC2, matriz geral) ou pelo algoritmo de Thomas (solveTridiagonal, PVC1 e problema final, matriz tridiagonal - muito mais barato que Gauss genérico).

Elementos Finitos (finiteElements.h): formulação fraca de Galerkin com elementos lineares (PVC1 e problema final) ou bilineares (PVC2), integrando por partes os termos de derivada mais alta antes de montar o sistema. O problema final tem um termo de 1ª ordem (7u') que não é auto-adjunto, então a montagem da forma fraca trata esse termo separadamente do termo de 2ª ordem.

A CLI permite resolver com um N (número de partições/elementos) escolhido pelo usuário, rodar as comparações fixas pelas tarefas (N=4 vs N=8) e exportar a solução do problema final em CSV para plotagem externa.

Módulo 6 - Processamento de Imagens (Detecção de Bordas)

Arquivos: src/imageProc/, include/imageProc/.

Módulo à parte do restante da disciplina, mas construído sobre os mesmos operadores de diferenças finitas: leitura/escrita de imagem em qualquer formato comum é feita pela stb_image/stb_image_write (ver Bibliotecas de terceiros); a imagem carregada vira uma matriz GrayImage de double em escala de cinza, com 0–255 representando intensidade (sem normalização para [0,1]).

• convolution.cpp: convolução 2D genérica usada por todos os filtros abaixo, com tratamento de borda por replicação do pixel mais próximo (atClamped).
• gaussianFilter.cpp: kernel Gaussiano normalizado (soma dos pesos = 1) para suavizar a imagem antes de derivar - reduz a sensibilidade dos filtros de gradiente/Laplaciano a ruído de alta frequência.
• sobelFilter.cpp: kernels de Sobel em x e y, construídos como o produto de uma diferença central (derivada) numa direção por uma suavização triangular [1,2,1] na direção ortogonal.
• laplacianFilter.cpp: kernel de 5 pontos [[0,1,0],[1,-4,1],[0,1,0]], a soma das diferenças centrais de 2ª ordem em x e em y.
• edgeDetection.cpp: implementa os dois algoritmos de detecção de bordas descritos abaixo, compondo os filtros acima.

Algoritmo 1 - Gradiente de Sobel + limiar (threshold)

1. Suaviza a imagem com o filtro Gaussiano.
2. Calcula as derivadas de Sobel em x (matriz A) e em y (matriz B).
3. Calcula a magnitude do gradiente C = √(A² + B²).
4. Gera a imagem final D: D = 0 onde C < threshold, D = 1 (borda) caso contrário.

Algoritmo 2 - Laplaciano + tolerância

1. Suaviza a imagem com o filtro Gaussiano.
2. Convolui com o kernel de Laplace, gerando a matriz A.
3. Gera a imagem final B: B = 1 (borda) onde |A| está fora de uma tolerância de zero, B = 0 caso contrário.

Em ambos, o resultado final é uma matriz binária (0/1) que é escalada para 0/255 apenas no momento de salvar, com a convenção de que borda = pixel preto e fundo = pixel branco.

Como a imagem não é normalizada, o parâmetro de sensibilidade de cada algoritmo precisa ser escolhido na escala 0–255, e não em frações próximas de zero:

Parâmetro	Faixa inicial recomendada
Kernel Gaussiano	5 (3 se a imagem já for limpa, 7–9 se houver bastante ruído)
Sigma da Gaussiana	1.0 – 1.5
Threshold (Algoritmo 1)	50 – 100
Tolerância (Algoritmo 2)	5 – 20

O Laplaciano é uma derivada de 2ª ordem, então é bem mais sensível a ruído do que o gradiente de Sobel; tolerâncias muito pequenas (frações de 1) tendem a marcar como "borda" qualquer variação mínima de textura, deixando o resultado praticamente todo branco antes da inversão de cor.

Parser de expressões matemáticas

Arquivos: src/parser/, include/parser/, biblioteca ExprTk em third_party/.

Os módulos de Derivação Numérica e PVI aceitam, além dos casos de teste pré-programados, digitar a função (f(x)) ou EDO escalar (dy/dt = f(t, y)) diretamente na CLI, sem recompilar. mathExpression.cpp encapsula a integração com o ExprTk atrás de uma interface própria (MathExpression), e functionBuilder.cpp conecta essa interface aos tipos usados pelo resto do projeto (Func1D, StateFunc), incluindo o prompt de leitura e a validação de sintaxe na CLI.

Sintaxe aceita: operadores aritméticos (+ - * / % ^), multiplicação implícita (2x, 3(x+1)), comparadores e lógicos, condicionais (if(cond, a, b) ou cond ? a : b), as funções matemáticas usuais (sin, cos, exp, log, sqrt, etc. - log é logaritmo natural; use log10/logn para outras bases) e as variáveis x (funções 1D) ou t/y (EDO escalar). Atribuições, laços e declaração de variáveis locais são desabilitados de propósito, já que a expressão pode ser avaliada milhares de vezes por segundo dentro dos laços numéricos. Erros de sintaxe são reportados com a posição do problema no texto, sem derrubar a execução.

Bibliotecas de terceiros

Vendorizadas em third_party/, sem necessidade de instalação separada:

• stb_image / stb_image_write - leitura e escrita de imagens em formatos comuns (PNG, JPG, BMP, TGA), usadas pelo módulo de processamento de imagens.
• ExprTk - parser e avaliador de expressões matemáticas em tempo real (header-only, licença MIT), usado pela entrada de funções/EDOs personalizadas.