diff --git a/examples/cpp/wave1D_convergence.cpp b/examples/cpp/wave1D_convergence.cpp new file mode 100644 index 00000000..77735781 --- /dev/null +++ b/examples/cpp/wave1D_convergence.cpp @@ -0,0 +1,171 @@ +/* + * ====================================================================================== + * Example: 2nd Order Convergence for 1D Wave Equation (Hyperbolic) + * Language: C++ (using MOLE library & Armadillo) + * ====================================================================================== + * + * Reference: + * Problem based on "Example 10.1" from: + * Mathews, J. H., & Fink, K. D. (2004). Numerical methods using MATLAB (4th ed.). + * Pearson Prentice Hall. + * + * Context: + * This example was presented during the postgraduate course "Introduction + * to Mimetic Difference Methods and Applications", taught by Prof. Jose + * Castillo in October 2023 at the Faculty of Exact Sciences, Engineering + * and Surveying (FCEIA) of the National University of Rosario (UNR), + * Argentina. + * + * Mathematical Formulation: + * ∂²u/∂t² = 4 ∂²u/∂x² on Ω = [0, 1] x [0, 0.5] + * (Comparing with standard form u_tt = c² u_xx, this implies wave speed c = 2) + * + * Domain Description: + * - Spatial: x ∈ [0, 1] + * - Temporal: t ∈ [0, 0.5] + * - Grid: Staggered grid (Mimetic Discretization) + * + * Boundary Conditions (Dirichlet): + * u(0, t) = 0 + * u(1, t) = 0 + * + * Initial Conditions: + * u(x, 0) = sin(πx) + sin(2πx) + * ∂u/∂t(x, 0) = 0 + * + * Analytical Solution (Exact): + * u(x, t) = sin(πx)cos(2πt) + sin(2πx)cos(4πt) + * + * Implementation Details: + * - Spatial Scheme: Mimetic Finite Differences (Order k=2) + * - Time Integration: Verlet Algorithm (Symplectic, 2nd order, Leapfrog equivalent) + * - Library: MOLE (with Armadillo for linear algebra) + * + * Output: + * - Prints a table of L2 errors and convergence rates for successive grid refinements. + * + * Author: Martin S. Armoa + * Programming Assistant: Google Gemini via VS Code + * ====================================================================================== + */ + +#include "mole.h" +#include +#include +#include +#include + +// Use Armadillo namespace +using namespace arma; +using namespace std; + +const double PI = 3.14159265358979323846; +const double WAVE_SPEED_C = 2.0; +const double T_FINAL = 0.5; + +// Analytical Solution +double exact_sol(double x, double t) { + return sin(PI*x)*cos(PI*WAVE_SPEED_C*t) + sin(2*PI*x)*cos(2*PI*WAVE_SPEED_C*t); +} + +// Simulation Runner (Accepts fixed dt) +double run_simulation(int m, double dt) { + int k = 2; // Spatial order + double a = 0.0; // Left boundary + double b = 1.0; // Right boundary + double dx = (b - a) / m; + + // Nt calculated from fixed dt + int Nt = (int)ceil(T_FINAL / dt); + + // 1. Initialize Mimetic Operator + Laplacian L(k, m, dx); + + // 2. State Vectors + vec u(m + 2, fill::zeros); + vec v(m + 2, fill::zeros); + vec acc(m + 2, fill::zeros); + vector x_centers(m + 2); + + // 3. Initial Conditions + for (int i = 1; i <= m; i++) { + double x = a + (i - 0.5) * dx; + x_centers[i] = x; + u(i) = exact_sol(x, 0.0); + } + // Dirichlet BC + u(0) = 0.0; + u(m+1) = 0.0; + + // 4. Time Integration (Velocity Verlet) + // Initial Acc + acc = (WAVE_SPEED_C*WAVE_SPEED_C) * (L * u); + acc(0) = 0.0; acc(m+1) = 0.0; // BC + + for (int t = 0; t < Nt; t++) { + v = v + 0.5 * dt * acc; + u = u + dt * v; + acc = (WAVE_SPEED_C*WAVE_SPEED_C) * (L * u); + acc(0) = 0.0; acc(m+1) = 0.0; // BC + v = v + 0.5 * dt * acc; + } + + // 5. L2 Error + double sum_sq_error = 0.0; + for (int i = 1; i <= m; i++) { + double diff = u(i) - exact_sol(x_centers[i], T_FINAL); + sum_sq_error += diff * diff; + } + + return sqrt(sum_sq_error * dx); +} + +int main() { + vector mesh_sizes = {20, 40, 80, 160, 320}; + + // --- CALC FIXED DT --- + // Match MATLAB logic: based on finest mesh (m=320) + int m_finest = 320; + double dx_min = 1.0 / m_finest; + // CFL Safety factor 0.25 + double dt_fixed = 0.25 * dx_min / WAVE_SPEED_C; + + cout << "Running C++ Convergence Test (Fixed dt)" << endl; // Changed title to verify update + cout << "---------------------------------------" << endl; + cout << "Configuration: Fixed dt = " << scientific << dt_fixed << endl << endl; + + cout << "### Convergence Rate Table (C++)" << endl; + cout << "| Cells (m) | dx | L2 Error | Rate (p) |" << endl; + cout << "| :--- | :--- | :--- | :--- |" << endl; + + vector errors; + vector dx_vals; + + for (size_t i = 0; i < mesh_sizes.size(); i++) { + int m = mesh_sizes[i]; + double dx = 1.0 / m; + + try { + double error = run_simulation(m, dt_fixed); + errors.push_back(error); + dx_vals.push_back(dx); + + string rate_str = "-"; + if (i > 0) { + double rate = log(errors[i-1] / errors[i]) / log(dx_vals[i-1] / dx_vals[i]); + char buffer[50]; + sprintf(buffer, "%.2f", rate); + rate_str = string(buffer); + } + + cout << "| " << setw(9) << left << m + << " | " << scientific << setprecision(4) << dx + << " | " << scientific << setprecision(4) << error + << " | " << setw(8) << rate_str << " |" << endl; + + } catch (const std::exception& e) { + cout << "Error simulation m=" << m << ": " << e.what() << endl; + } + } + return 0; +} diff --git a/examples/matlab_octave/wave1D_convergence/run_convergence_test.m b/examples/matlab_octave/wave1D_convergence/run_convergence_test.m new file mode 100644 index 00000000..bf6eabd6 --- /dev/null +++ b/examples/matlab_octave/wave1D_convergence/run_convergence_test.m @@ -0,0 +1,143 @@ +% ========================================================================= +% Example: 2nd Order Convergence for 1D Wave Equation (Hyperbolic) +% Language: MATLAB / Octave +% ========================================================================= +% +% Reference: +% Problem based on "Example 10.1" from: +% Mathews, J. H., & Fink, K. D. (2004). Numerical methods using MATLAB +% (4th ed.). Pearson Prentice Hall. +% +% Context: +% This example was presented during the postgraduate course "Introduction +% to Mimetic Difference Methods and Applications", taught by Prof. Jose +% Castillo in October 2023 at the Faculty of Exact Sciences, Engineering +% and Surveying (FCEIA) of the National University of Rosario (UNR), +% Argentina. +% +% Mathematical Formulation: +% d2u/dt2 = 4 * d2u/dx2 on Omega = [0, 1] x [0, 0.5] +% (Comparing with standard form u_tt = c^2 u_xx, this implies c = 2) +% +% Domain Description: +% - Spatial: x in [0, 1] +% - Temporal: t in [0, 0.5] +% - Grid: Staggered grid (Mimetic Discretization) +% +% Boundary Conditions (Dirichlet): +% u(0, t) = 0 +% u(1, t) = 0 +% +% Initial Conditions: +% u(x, 0) = sin(pi*x) + sin(2*pi*x) +% du/dt(x, 0) = 0 +% +% Analytical Solution (Exact): +% u(x, t) = sin(pi*x)cos(2*pi*t) + sin(2*pi*x)cos(4*pi*t) +% +% Implementation Details: +% - Spatial Scheme: Mimetic Finite Differences (Order k=2) vs Standard FD +% - Time Integration: Verlet Algorithm (Symplectic, 2nd order, Leapfrog equivalent) +% - Library: MOLE (MATLAB/Octave implementation) +% +% Output: +% - Console: Table of L2 errors and convergence rates. +% - Figure 1: Error vs Grid Spacing (dx) [Comparison] +% - Figure 2: Error vs Number of Cells (m) [Comparison] +% - Figure 3: Wave Profile Comparison (Coarse grid m=20 to visualize error) +% +% Author: Martin S. Armoa +% Programming Assistant: Google Gemini 3 PRO via VS Code +% ========================================================================= +clear; clc; close all; + +% --- 1. Path Configuration --- + +addpath('../../../src/matlab_octave'); +mole_lib_path = fullfile(current_folder, '..', '..', '..', 'src', 'matlab_octave'); +if exist(mole_lib_path, 'dir'), addpath(mole_lib_path); end + +fprintf('Running Comparative Convergence Test (Mimetic vs FD)\n'); +fprintf('--------------------------------------------------\n'); + +mesh_sizes = [20, 40, 80, 160, 320]; +n_sims = length(mesh_sizes); +c = 2; +m_finest = max(mesh_sizes); +dx_finest_mimetic = 1.0 / m_finest; +% CFL condition based on finest mesh: dt <= dx / c. +% We use a safety factor of 0.25 +dt_fixed = 0.25 * dx_finest_mimetic / c; + +fprintf('Configuration: Fixed dt = %.5e (Ensures stability for all meshes)\n', dt_fixed); + +results = struct('dx_mim', zeros(n_sims,1), 'dx_fd', zeros(n_sims,1), ... + 'm', zeros(n_sims,1), ... + 'err_mim', zeros(n_sims,1), 'rate_mim', zeros(n_sims,1), ... + 'err_fd', zeros(n_sims,1), 'rate_fd', zeros(n_sims,1)); + +for i = 1:n_sims + m = mesh_sizes(i); + results.m(i) = m; + + % --- MIMETIC SOLVER (Staggered Grid) --- + [results.err_mim(i), results.dx_mim(i)] = wave1d_solver(m, dt_fixed); + + % --- FD SOLVER (Nodal Grid, m+2 points) --- + [results.err_fd(i), results.dx_fd(i)] = wave1d_solver_fd(m, dt_fixed); + + % Compute Convergence Rates + if i > 1 + % Rate = log(E_prev/E_curr) / log(dx_prev/dx_curr) + log_dx_mim = log(results.dx_mim(i-1) / results.dx_mim(i)); + results.rate_mim(i) = log(results.err_mim(i-1) / results.err_mim(i)) / log_dx_mim; + + log_dx_fd = log(results.dx_fd(i-1) / results.dx_fd(i)); + results.rate_fd(i) = log(results.err_fd(i-1) / results.err_fd(i)) / log_dx_fd; + end +end + +% --- Display Table --- +fprintf('\n### Convergence Comparison: Mimetic vs Standard FD\n'); +fprintf('Note: Time step dt is fixed. Comparison targets spatial order ~2.0.\n\n'); +fprintf('| Cells(m) | Mimetic Err | Rate | FD Error | Rate |\n'); +fprintf('| :--- | :--- | :--- | :--- | :--- |\n'); + +for i = 1:n_sims + r_mim = '-'; r_fd = '-'; + if i > 1 + r_mim = sprintf('%.2f', results.rate_mim(i)); + r_fd = sprintf('%.2f', results.rate_fd(i)); + end + fprintf('| %-8d | %.3e | %s | %.3e | %s |\n', ... + mesh_sizes(i), results.err_mim(i), r_mim, results.err_fd(i), r_fd); +end + +% --- Optional Plotting +%{ +if ~exist('OCTAVE_VERSION', 'builtin') || ~isempty(getenv('DISPLAY')) + % Figure 1: Error vs Grid Spacing (dx) + figure(1); + loglog(results.dx_mim, results.err_mim, '-o', 'DisplayName', 'Mimetic'); + hold on; + loglog(results.dx_fd, results.err_fd, '-x', 'DisplayName', 'FD'); + % Reference line O(h^2) using Mimetic dx + ref_y = results.err_mim(end) * (results.dx_mim / results.dx_mim(end)).^2; + loglog(results.dx_mim, ref_y, '--k', 'DisplayName', 'O(h^2)'); + grid on; xlabel('dx'); ylabel('L2 Error'); legend('Location', 'best'); + title('Convergence Comparison: Error vs dx'); + + % Figure 2: Error vs Number of Cells (m) + figure(2); + loglog(results.m, results.err_mim, '-o', 'DisplayName', 'Mimetic'); + hold on; + loglog(results.m, results.err_fd, '-x', 'DisplayName', 'FD'); + % Reference line O(m^-2) + ref_y_m = results.err_mim(1) * (results.m / results.m(1)).^(-2); + loglog(results.m, ref_y_m, '--k', 'DisplayName', 'O(m^{-2})'); + grid on; xlabel('Cells (m)'); ylabel('L2 Error'); legend('Location', 'best'); + title('Convergence Comparison: Error vs Cells'); + + fprintf('\nPlots are commented out by default. Uncomment in script to view.\n'); +end +%} diff --git a/examples/matlab_octave/wave1D_convergence/wave1d_solver.m b/examples/matlab_octave/wave1D_convergence/wave1d_solver.m new file mode 100644 index 00000000..7b385fcb --- /dev/null +++ b/examples/matlab_octave/wave1D_convergence/wave1d_solver.m @@ -0,0 +1,103 @@ +function [L2_error, dx, u_num_centers, x_centers] = wave1d_solver(m,dt) +% WAVE1D_SOLVER Solves the 1D Wave Equation using Mimetic Finite Differences. +% +% formulation: d2u/dt2 = c^2 * div(grad(u)) +% scheme: Spatial: 2nd Order Mimetic Operators (MOLE library) +% Temporal: Velocity Verlet (Leapfrog) +% +% Input: +% m : Number of cells in the spatial grid. +% +% Outputs: +% L2_error : Discrete L2 norm error against exact solution. +% dx : Grid spacing. +% u_num_centers : Numerical solution vector (interior points). +% x_centers : Coordinate vector of cell centers. + + % ========================================================================= + % 1. Physical Parameters and Grid Configuration + % ========================================================================= + k = 2; % Spatial order of accuracy + a = 0; b = 1; % Domain boundaries [0, 1] + dx = (b-a)/m; % Grid spacing + c = 2; % Wave speed + T_final = 0.5; % Simulation duration + Nt = ceil(T_final/dt); + + % Coordinate system: Cell centers (staggered grid) + % Used for evaluating initial functions on the interior nodes. + x_centers = (a + dx/2 : dx : b - dx/2)'; % Column vector of size m + + % ========================================================================= + % 2. Mimetic Operators Setup + % ========================================================================= + % Initialize Mimetic Laplacian operator + % Size is (m+2) x (m+2) to accommodate boundary values (ghost points) + L = lap(k, m, dx); + + % Impose Boundary Conditions on the operator matrix (Dirichlet u=0) + % This ensures the operator behaves correctly at the edges. + L(1, :) = 0; L(1, 1) = 1; % Left boundary (Identity row) + L(end, :) = 0; L(end, end) = 1; % Right boundary (Identity row) + + % Define the Force Function (RHS) + % F = c^2 * Laplacian * u + F_op = @(u) (c^2) * (L * u); + + % ========================================================================= + % 3. Initial Conditions (Augmented Vector) + % ========================================================================= + % Initialize state vector 'u' with size m+2 (includes boundaries) + u = zeros(m+2, 1); + + % Populate interior nodes (indices 2 to m+1) + ICU = @(x) sin(pi*x) + sin(2*pi*x); + u(2:end-1) = ICU(x_centers); + + % Enforce Boundary Conditions (indices 1 and m+2) + u(1) = 0; + u(end) = 0; + + % Initialize velocity vector 'v' (du/dt = 0 initially) + v = zeros(m+2, 1); + + % ========================================================================= + % 4. Time Integration (Velocity Verlet / Leapfrog) + % ========================================================================= + + % Calculate initial acceleration + acc = F_op(u); + + % Enforce zero acceleration at boundaries to maintain strict Dirichlet BCs + acc(1) = 0; acc(end) = 0; + + for t = 1:Nt + % a) Half-step velocity update + v = v + 0.5 * dt * acc; + + % b) Full-step position update + u = u + dt * v; + + % c) Recalculate forces (acceleration) based on new position + acc = F_op(u); + + acc(1) = 0; acc(end) = 0; + + % d) Full-step velocity update + v = v + 0.5 * dt * acc; + end + + % ========================================================================= + % 5. Error Analysis + % ========================================================================= + % Analytical solution evaluated at cell centers + u_exact_centers = sin(pi*x_centers)*cos(pi*c*T_final) + ... + sin(2*pi*x_centers)*cos(2*pi*c*T_final); + + % Extract the interior of the numerical solution for comparison + u_num_centers = u(2:end-1); + + % Compute discrete L2 Error Norm + diff = u_num_centers - u_exact_centers; + L2_error = sqrt(sum(diff.^2) * dx); +end diff --git a/examples/matlab_octave/wave1D_convergence/wave1d_solver_fd.m b/examples/matlab_octave/wave1D_convergence/wave1d_solver_fd.m new file mode 100644 index 00000000..70a21ed6 --- /dev/null +++ b/examples/matlab_octave/wave1D_convergence/wave1d_solver_fd.m @@ -0,0 +1,95 @@ +function [L2_error, dx, u_num_centers, x_centers] = wave1d_solver_fd(m,dt) +% WAVE1D_SOLVER_FD Solves the 1D Wave Equation using Standard Finite Differences. +% +% formulation: d2u/dt2 = c^2 * d2u/dx2 +% scheme: Spatial: Standard Central Differences (3-point stencil) +% Temporal: Velocity Verlet (Leapfrog) +% context: Acts as a baseline to validate/compare Mimetic methods. +% +% Input: +% m : Number of cells in the spatial grid. +% +% Outputs: +% L2_error : Discrete L2 norm error against exact solution. +% dx : Grid spacing. +% u_num : Numerical solution vector +% + + % ========================================================================= + % 1. Physical Parameters and Grid Configuration + % ========================================================================= + a = 0; b = 1; % Domain boundaries + dx = (b-a)/m; % Grid spacing + c = 2; % Wave speed + T_final = 0.5; % Simulation duration + Nt = ceil(T_final/dt); + % Comparison: Mimetic uses m cells. FD uses m+2 nodes to cover same domain. + N = m + 2; + x_nodes = linspace(a, b, N)'; + dx = x_nodes(2) - x_nodes(1); % Exact nodal spacing + + % ========================================================================= + % 2. Standard Finite Difference Operator (Tridiagonal Matrix) + % ========================================================================= + + e = ones(N, 1); + + % Construct sparse tridiagonal matrix with diagonals [1, -2, 1] + % -1: Lower diagonal + % 0: Main diagonal + % +1: Upper diagonal + L = spdiags([e -2*e e], -1:1, N, N); + L = L / dx^2; + + % Boundary Conditions (Explicit Dirichlet u=0) + % We zero out the first and last rows to prevent the stencil from + % evolving the boundary nodes. They remain fixed at 0. + L(1, :) = 0; L(1, 1) = 0; + L(end, :) = 0; L(end, end) = 0; + + % Define the Force Function (RHS) + F_op = @(u) (c^2) * (L * u); + + % ========================================================================= + % 3. Initial Conditions + % ========================================================================= + u = zeros(N, 1); + + % Initialize interior points + ICU = @(x) sin(pi*x) + sin(2*pi*x); + u = ICU(x_nodes); + + % Initialize velocity (du/dt = 0) + v = zeros(N, 1); + + % ========================================================================= + % 4. Time Integration (Velocity Verlet) + % ========================================================================= + acc = F_op(u); + + for t = 1:Nt + % a) Half-step velocity update + v = v + 0.5 * dt * acc; + + % b) Full-step position update + u = u + dt * v; + + % c) Recalculate acceleration + acc = F_op(u); + + % d) Full-step velocity update + v = v + 0.5 * dt * acc; + end + + % ========================================================================= + % 5. Error Analysis + % ========================================================================= + % Analytical solution evaluated at nodes + u_exact_all = sin(pi*x_nodes)*cos(pi*c*T_final) + sin(2*pi*x_nodes)*cos(2*pi*c*T_final); + + u_num_internal = u(2:end-1); + u_exact_internal = u_exact_all(2:end-1); + + diff = u_num_internal - u_exact_internal; + L2_error = sqrt(sum(diff.^2) * dx); +end