Skip to content

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

Open
tin7 wants to merge 2 commits intocsrc-sdsu:mainfrom
tin7:feature/second-order-convergence
Open

Feat: Add 1D Wave convergence example (Matlab & C++) ref Mathews & Fi…#270
tin7 wants to merge 2 commits intocsrc-sdsu:mainfrom
tin7:feature/second-order-convergence

Conversation

@tin7
Copy link

@tin7 tin7 commented Jan 22, 2026

…nk Ex 10.1

What type of PR is this? (check all applicable)

  • Refactor
  • Feature
  • Bug Fix
  • Optimization
  • Example
  • Documentation

Description

Academic example comparing Mimetic Finite Differences vs Standard Finite Differences for the 1D Wave Equation (Hyperbolic).

Context & Reference

  • Source Problem: Example 10.1 from Numerical Methods Using MATLAB (Mathews & Fink, 2004).
  • Institutional Context: Example presented in the postgraduate course "Introduction to Mimetic Difference Methods and Applications" (Prof. Jose Castillo, UNR FCEIA, Oct 2023).

Mathematical Formulation

$$\frac{\partial^2 u}{\partial t^2} = 4 \frac{\partial^2 u}{\partial x^2}, \quad x \in [0, 1], \quad t \in [0, 0.5]$$

  • BC: Dirichlet $u(0,t)=0, u(1,t)=0$.
  • IC: $u(x,0) = \sin(\pi x) + \sin(2\pi x)$.

Implementation Details

  1. C++ (examples/cpp/wave1D_convergence.cpp):

    • Uses MOLE library
    • Time integration: Velocity Verlet.
    • Demonstrates global 2nd-order convergence (L2 norm).
  2. MATLAB/Octave (examples/matlab_octave/wave1D_convergence/):

    • run_convergence_test.m: Driver script.
    • Generates comparative plots (see QA Instructions below).

Numerical Verification

Mimetic Scheme Convergence Rates (C++ Output):

Cells (m) dx L2 Error Rate
20 5.0000e-02 4.3367e-04 -
40 2.5000e-02 4.9180e-05 3.14
80 1.2500e-02 5.8741e-06 3.07
160 6.2500e-03 6.4405e-07 3.19
320 3.1250e-03 5.3138e-08 3.60

Note: The observed rates > 2.0 indicate superconvergence effects typical of Mimetic operators on staggered grids with smooth initial conditions.

Related Issues & Documents

  • Closes #

QA Instructions, Screenshots, Recordings

Visual validation of the implementation:

(1) Convergence: Error vs Grid Spacing
convergence_vs_dx

(2) Convergence: Error vs cells
convergence_vs_cells

(3) Wave Profile Comparison (Coarse Grid m=20)
Visual validation of phase and amplitude accuracy on a low-resolution grid. It is noteworthy that the Mimetic method closely tracks the analytical solution even on this coarse mesh ($m=20$), exhibiting significantly less numerical dispersion than the standard Finite Difference scheme.
wave_profile_coarse

How to test:

  1. C++: cd build && make && ./examples/cpp/wave1D_convergence
  2. MATLAB: Run examples/matlab_octave/wave1D_convergence/run_convergence_test.m

Added/updated tests?

_We encourage you to test all code included with MOLE, including examples.

  • Yes
  • No, and this is why: please replace this line with details on why tests
    have not been included
  • I need help with writing tests

Read Contributing Guide and Code of Conduct

[optional] Are there any post deployment tasks we need to perform?

[optional] What gif best describes this PR or how it makes you feel?

Copy link
Collaborator

@jbrzensk jbrzensk left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for contributing @tin7 . It will be good to have this comparison.

For the code, I am going to reference the MATLAB version, but the same suggestion may apply to the C++ code.

In the MATLAB folder, you have div, grad, lap. Please use the operators in the MOLE src/matlab_octave folder.

Also, you do not need the pictures in the PR. Those are generated by the user. You can include them in the pull request, so we can see what the results look like, but leave them out of the PR.

To check the convergence of the operator in space, you FIX the time discretization (dt) and vary the space discretization (dx). You have both varying.

Also, for the FD method, you have a staggered grid (line 34). To do a correct comparison, you need m+2 equally spaced points. The finite difference does not work on a staggered grid. Currently, this introduces additional error in the FD method. When done correctly, both of these measurements (FD and MD methods) should be really close to 2.

@tin7
Copy link
Author

tin7 commented Jan 23, 2026

@jbrzensk

Thank you for the detailed review and the suggestions regarding the grid setup and time stepping.
Changes applied:

  1. Fixed Time Step: Logic updated in both languages. dt is now calculated based on the finest mesh ($m=320$) and kept constant for all simulations to strictly isolate spatial convergence error.

  2. FD Grid: The MATLAB FD solver now uses a nodal grid (linspace) with $m+2$ points to avoid staggered grid interpolation errors.

  3. Removed plotting files and local operators from the PR. Plotting code is preserved but commented out in the script.

Ready for re-review.

Copy link
Collaborator

@jbrzensk jbrzensk left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Getting better. I put some notes about the paths and the error. The mimetic is order 2, you should get order 2, like really close to 2, (1.99->2.01). Same for the finite difference. If you get something different, you did something wrong.

Some other notes about adding paths much simpler, and using the dot operator for scalar times matrix.

Again, thanks for this. I will probably have notes about naming the files and their folder, but thats minor. Key is to get the error to show around order 2. Once it does, change the order of the mimetic operator and ensure that it is order 4. Once yuou fix the matlab, I just assume you will fix the cpp code.

Here is what I got with a little tinkering of your code:
image

clear; clc; close all;

% --- 1. Path Configuration ---
current_file = mfilename('fullpath');
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Files in the current folder are automatically in the path for both Matlab and Octave.
You can add the src files to the path like other matlab examples, with

addpath('../../../src/matlab_octave');

mesh_sizes = [20, 40, 80, 160, 320];
n_sims = length(mesh_sizes);
c = 2;
m_finest = max(mesh_sizes);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just set a small dt. No need to to complicated logic to pick a value. You already picked the numbers of cells, pick an appropriate dt. Something that devides evenly into 0.5

fprintf('| :--- | :--- | :--- | :--- | :--- |\n');

for i = 1:n_sims
r_mim = '-'; r_fd = '-';
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use '----' to make the table print out nicer


% Define the Force Function (RHS)
% F = c^2 * Laplacian * u
F_op = @(u) (c^2) * (L * u);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

scalar times matrix should have '.*' as the operator. MATLAB usually handles this OK, but its also a clue that its a scalar and matrix in the equation


% Compute discrete L2 Error Norm
diff = u_num_centers - u_exact_centers;
L2_error = sqrt(sum(diff.^2) * dx);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

L2 error has N, but I think you want relative error here, for testing convergence. sqrt((Uh-U)^2) / sqrt(U^2).

% 1. Physical Parameters and Grid Configuration
% =========================================================================
a = 0; b = 1; % Domain boundaries
dx = (b-a)/m; % Grid spacing
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

dx set here, and changed later. Just use the later one.

% 0: Main diagonal
% +1: Upper diagonal
L = spdiags([e -2*e e], -1:1, N, N);
L = L / dx^2;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

./ for division. For sanity.

L(end, :) = 0; L(end, end) = 0;

% Define the Force Function (RHS)
F_op = @(u) (c^2) * (L * u);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

.* for multiply

u_exact_internal = u_exact_all(2:end-1);

diff = u_num_internal - u_exact_internal;
L2_error = sqrt(sum(diff.^2) * dx);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here, I think you want relative error.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants