Stage-relabeling invariance for Runge–Kutta methods

In normal Runge–Kutta / numerical methods language, the idea is usually described like this:

1) Runge–Kutta stage reordering does not change the method

An explicit Runge–Kutta method works by doing several internal “sub-steps” (called stages) inside one time step.
The theorem says:

If you relabel or reorder those stages (stage 1 becomes stage 2, stage 2 becomes stage 1, etc.)

and you also rearrange the method’s coefficient tables to match that new ordering then the method still produces the same final result for the step.

So even though the internal bookkeeping changes, the actual computed step result does not change (except for tiny rounding differences).

This is commonly called: Stage permutation invariance or invariance under stage relabeling

2) This still stays true when you chain methods

Your code doesn’t just run one RK method; it chains two of them: Run method 2 for one step, then run method 1 for one step

The theorem also says:

If you replace method 1 with an equivalent “stage-reordered version” and/or replace method 2 with an equivalent “stage-reordered version”
then the chained result is still the same.

This is often described as:

Composition is compatible with stage reordering or more formally: composition is well-defined under the equivalence

If you want one “research-style umbrella description” (still readable):

Reordering stages does not change RK results, and it also does not change the results of chaining RK methods.

Short introduction to theorem_382A (in terms of your code)

Your code is testing this exact claim:

You build two normal RK methods:

Heun RK2 (2 stages) Classical RK4 (4 stages)

You build reordered versions of the same methods by shuffling their stage order, while also shuffling the tables correctly.

You run both versions on the same ODE step and compare:

Original chain: RK4 step then Heun step

Reordered chain: permuted RK4 step then permuted Heun step

The theorem predicts:

The two final answers should match (tiny floating-point rounding differences are okay)

So the code shows you can reorder internal stages for implementation reasons (like memory layout or simpler storage) without changing what the method computes.

How demo_theorem_382A_5cases_main.c shows this, step by step
Step 1 — Build the original RK methods

Your code creates two well-known methods:

make_heun_rk2() fills M1_s, M1_A, M1_b, M1_c
This represents Heun’s RK2.

make_rk4() fills M2_s, M2_A, M2_b, M2_c
This represents classical RK4.

Each method is stored in the standard RK “tableau” parts:

s = how many stages the method has

c[] = where each stage evaluates time inside the step

b[] = how the final result mixes the stage derivatives

A = how later stages use earlier stage derivatives

Important detail in your implementation:
You store only the lower triangle of A (Atri), because explicit RK never needs “future stage depends on future stage” terms.

That’s why you use:

Aidx(i,j), getA(), setA() for triangular indexing.

Step 2 — Build equivalent “stage-reordered” versions

Inside buildAll() you define permutations:

p1 reorders Heun’s stages

p2 reorders RK4’s stages

Then permute_method_into(...) produces:

M1B_* from M1_*

M2B_* from M2_*

What that function does (in plain words):

It rearranges b[] and c[] so stage i in the new method corresponds to the correct original stage.

It also rearranges the stored A entries so that “stage i depends on stage j” is preserved under the new labels.

So the permuted versions are not “different methods.”
They are the same methods written with a different stage numbering.

That is exactly what theorem_382A is about.

Step 3 — Define composition exactly as the theorem describes

The theorem talks about chaining methods.

Your code implements chaining here:

var compose_step(...)
{
var y_mid = rk_step(method2, ...);
var y_out = rk_step(method1, ..., y_mid, ...);
return y_out;
}


So composition means:

Apply method 2 for one step ? get intermediate result

Apply method 1 for one step ? get final result

Then your comparison is:

yA = composition using original methods (M1 and M2)

yB = composition using permuted methods (M1B and M2B)

The theorem predicts: yA and yB should match.

Step 4 — Test it on 5 different ODE examples

run_case(k, ...) sets G_CASE = k, which changes what f_ode(t,y) returns.

You test five different behaviors: simple decay, simple growth, faster decay (stiffer), nonlinear growth with saturation (logistic)
forcing term that depends on time.

This shows the theorem isn’t only working for one “easy” equation — it works across different ODE styles.

Step 5 — What the printed output means

Each test prints: yA result (original chain), yB result (permuted chain), diff = how far apart they are
If theorem_382A is correct (and your permutation is applied correctly): diff should be extremely small (often near machine precision)