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

examples/cpp/wave1D_convergence.cpp

@@ -0,0 +1,192 @@
+/*
+ * ======================================================================================
+ * 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 relative L2 errors and convergence rates for successive grid refinements.
+ *
+ * Author: Martin S. Armoa
+ * ======================================================================================
+ */
+
+#include "mole.h"
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <iomanip>
+#include <sstream>
+
+using namespace arma;
+using namespace std;
+
+const double PI = 3.14159265358979323846;
+const double WAVE_SPEED_C = 2.0;
+const double T_TARGET = 0.35;
+
+// 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;
+    double t_eval = T_TARGET;
+    int Nt = (int)ceil(t_eval / dt);
+    double t_simulated = Nt * 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);
+    vec acc_old(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 Acceleration
+    acc = (WAVE_SPEED_C*WAVE_SPEED_C) * (L * u);
+    acc(0) = 0.0;
+    acc(m+1) = 0.0;
+    for (int t = 0; t < Nt; t++) {
+        // a) Full-step position update
+        u = u + dt * v + 0.5 * dt * dt * acc;
+
+        // Store previous acceleration
+        acc_old = acc;
+
+        // b) Recalculate acceleration at new position
+        acc = (WAVE_SPEED_C*WAVE_SPEED_C) * (L * u);
+        acc(0) = 0.0;
+        acc(m+1) = 0.0; // BC
+
+        // c) Full-step velocity update
+        v = v + 0.5 * dt * (acc_old + acc);
+
+        // Ensure boundaries remain strictly zero
+        u(0) = 0.0;
+        u(m+1) = 0.0;
+    }
+    // 5. Relative Error Calculation
+    double sum_sq_error = 0.0;
+    double sum_sq_exact = 0.0;
+
+    for (int i = 1; i <= m; i++) {
+        double u_ext = exact_sol(x_centers[i], t_simulated);
+        double diff = u(i) - u_ext;
+
+        sum_sq_error += diff * diff;
+        sum_sq_exact += u_ext * u_ext;
+    }
+
+    double norm_error = sqrt(sum_sq_error * dx);
+    double norm_exact = sqrt(sum_sq_exact * dx);
+
+    return norm_error / norm_exact;
+}
+int main() {
+    vector<int> mesh_sizes = {20, 40, 80, 160, 300, 400, 500};
+    double dt_fixed = 0.0001;
+
+    cout << "Running C++ Convergence Test (Fixed dt)" << endl;
+    cout << "---------------------------------------" << endl;
+    cout << "Configuration: Fixed dt = " << dt_fixed << endl << endl;
+
+    cout << "### Convergence Rate Table (C++)" << endl;
+    cout << "| Cells (m) | dx         | Rel 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);
+
+            ostringstream dx_stream;
+            dx_stream << scientific << setprecision(4) << dx;
+            string dx_str = dx_stream.str();
+
+            ostringstream err_stream;
+            err_stream << scientific << setprecision(4) << error;
+            string err_str = err_stream.str();
+
+            string rate_str = "----";
+            if (i > 0) {
+                double rate = log(errors[i-1] / errors[i]) / log(dx_vals[i-1] / dx_vals[i]);
+                ostringstream rate_stream;
+                rate_stream << fixed << setprecision(2) << rate;
+                rate_str = rate_stream.str();
+            }
+
+            cout << "| " << setw(9) << right << m
+                 << " | " << setw(10) << right << dx_str
+                 << " | " << setw(12) << right << err_str
+                 << " | " << setw(8) << right << rate_str << " |" << endl;
+
+        } catch (const std::exception& e) {
+            cout << "Error simulation m=" << m << ": " << e.what() << endl;
+        }
+    }
+    return 0;
+}

examples/matlab_octave/wave1D_convergence/wave1D_convergence.m

@@ -0,0 +1,157 @@
+% =========================================================================
+% 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 Relative 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;
+
+addpath('../../../src/matlab_octave');
+
+fprintf('Running Comparative Convergence Test (Mimetic vs FD)\n');
+fprintf('--------------------------------------------------\n');
+
+mesh_sizes = [20, 40, 80, 160, 300, 400, 500];
+n_sims = length(mesh_sizes);
+c = 2;
+dt_fixed = 0.0001;
+
+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
+
+##
+##    % 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('R 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('R Error'); legend('Location', 'best');
+##    title('Convergence Comparison: Error vs Cells');
+##
+##    fprintf('\nPlots are commented out by default. Uncomment in script to view.\n');
+##
+##
+##
+##
+##    % Figure 1: Error vs dx
+##    figure(1);
+##    loglog(results.dx_mim, results.err_mim, '-o', ...
+##        'LineWidth', 1.2, 'MarkerSize', 6, 'DisplayName', 'Mimetic');
+##    hold on;
+##    loglog(results.dx_fd, results.err_fd, '-x', ...
+##        'LineWidth', 1.2, 'MarkerSize', 6, 'DisplayName', 'FD');
+##    ref_y = results.err_mim(1) * (results.dx_mim / results.dx_mim(1)).^2;
+##    loglog(results.dx_mim, ref_y, '--k', ...
+##        'LineWidth', 1.2, 'DisplayName', 'O(h^2)');
+##    grid on;
+##    grid minor;
+##    xlabel('dx', 'FontSize', 11, 'FontWeight', 'bold');
+##    ylabel('R Error', 'FontSize', 11, 'FontWeight', 'bold');
+##    title('Convergence Comparison: Error vs dx', 'FontSize', 12, 'FontWeight', 'bold');
+##    legend('Location', 'southeast', 'FontSize', 10);
+##    hold off;
+####
+##
+