Feat: Add 1D Wave convergence example (Matlab & C++) ref Mathews & Fi…

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 <iostream>
+#include <vector>
+#include <cmath>
+#include <iomanip>
+
+// 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<double> 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<int> 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<double> errors;
+    vector<double> 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;
+}

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
+%}