diff --git a/examples/cpp/wave1D_convergence.cpp b/examples/cpp/wave1D_convergence.cpp new file mode 100644 index 00000000..346e6bd7 --- /dev/null +++ b/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 +#include +#include +#include +#include + +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 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 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 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); + + 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; +} diff --git a/examples/matlab_octave/wave1D_convergence/wave1D_convergence.m b/examples/matlab_octave/wave1D_convergence/wave1D_convergence.m new file mode 100644 index 00000000..43da9d8b --- /dev/null +++ b/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; +#### +## + 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..de75f19c --- /dev/null +++ b/examples/matlab_octave/wave1D_convergence/wave1D_solver.m @@ -0,0 +1,106 @@ +function [R_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: +% R _error : Discrete Relative 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.35; % 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 + L(end, :) = 0; + L(end, end) = 1; % Right boundary + + % Define the Force Function (RHS) + % F = c^2 * Laplacian * u + F_op = @(u) (c^2) .* (L * u); + + % ========================================================================= + % 3. Initial Conditions + % ========================================================================= + % 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; + u = u + dt * v + 0.5 * dt^2 * acc; + acc_old = acc; + + % b) 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_old + acc); + end + + % ========================================================================= + % 5. Error Analysis + % ========================================================================= + % Analytical solution + T_simulated = Nt * dt; + + % Analytical solution evaluated at cell centers using T_simulated + u_exact_centers = sin(pi*x_centers)*cos(pi*c*T_simulated) + ... + sin(2*pi*x_centers)*cos(2*pi*c*T_simulated); + u_num_centers = u(2:end-1); + % Compute discrete relative Error Norm + diff = u_num_centers - u_exact_centers; + norm_error = sqrt(sum(diff.^2) * dx); + norm_exact = sqrt(sum(u_exact_centers.^2) * dx); + R_error = norm_error / norm_exact; + +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..409ed58a --- /dev/null +++ b/examples/matlab_octave/wave1D_convergence/wave1D_solver_fd.m @@ -0,0 +1,101 @@ +function [R_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: +% R_error : Discrete relative 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 + c = 2; % Wave speed + T_final = 0.35; % 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; + u = u + dt * v + 0.5 * dt^2 * acc; + acc_old = acc; + + % b) 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_old + acc); + end + u(1) = 0; u(end) = 0; + + % ========================================================================= + % 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; + norm_error = sqrt(sum(diff.^2) * dx); + norm_exact = sqrt(sum(u_exact_all.^2) * dx); + R_error = norm_error / norm_exact; + +end + + +