option : show_navy_triangles := false;
Spirale_Suites_complexe_animatie_gemaaktDEF_2-5-2026.mw

KroneckerRules := module()
    option package;
    export Kron, `&x`, Simplify, Expand, Shuffle, Check, VisualCheck, Verify, VisualVerify, NumericEvaluate, Info, SimplifyTrace;
    
    # ---------- Inert notation ----------
    Kron := (A,B) -> 'KroneckerProduct'(A,B):
    `&x` := (A,B) -> 'KroneckerProduct'(A,B):
    
    # ---------- Trace mechanism ----------
    local _trace_steps := []:
    local _AddTrace := proc(rule, from_expr, to_expr)
        _trace_steps := [op(_trace_steps),
            sprintf("%s: **%a**  ->  **%a**", rule, from_expr, to_expr)];
    end proc;
    
    # ---------- Perfect shuffle matrix ----------
    local _PerfectShuffleMatrix := proc(p::posint, r::posint)
        local M, i, j;
        M := Matrix(p*r, p*r, 0);
        for i from 1 to p do
            for j from 1 to r do
                M[(j-1)*p + i, (i-1)*r + j] := 1;
            end do;
        end do;
        return M;
    end proc;
    
    # ---------- Helper: flatten a dot product into a list of operands ----------
    local _FlattenOperands := proc(p)
        local terms, rec;
        terms := [];
        rec := proc(x)
            if type(x, `.`) then
                rec(op(1,x));
                rec(op(2,x));
            else
                terms := [op(terms), x];
            end if;
        end proc;
        rec(p);
        return terms;
    end proc;
    
    # ---------- Distribute dot over sums (recursively, fully) ----------
    local _DistributeDot := proc(expr, with_trace::boolean := false)
        local e, changed, term;
        e := expr;
        changed := true;
        while changed do
            changed := false;
            e := subsindets(e, `.`, proc(p)
                local a,b,terms,res;
                a := op(1,p); b := op(2,p);
                if type(a, `+`) then
                    changed := true;
                    terms := [op(a)];
                    res := `+`(seq(term . b, term in terms));
                    if with_trace then _AddTrace("Dot left-distribution", p, res); end if;
                    return res;
                elif type(b, `+`) then
                    changed := true;
                    terms := [op(b)];
                    res := `+`(seq(a . term, term in terms));
                    if with_trace then _AddTrace("Dot right-distribution", p, res); end if;
                    return res;
                else
                    return p;
                end if;
            end proc);
        end do;
        return e;
    end proc;
    
    # ---------- Combine consecutive Kronecker products in a dot chain (full flatten) ----------
    local _FlattenChain := proc(expr, with_trace::boolean := false)
        local e, changed;
        e := expr;
        changed := true;
        while changed do
            changed := false;
            e := subsindets(e, `.`, proc(p)
                local ops, all_kron, left, right, i, to_pair;
                ops := [op(p)];
                all_kron := true;
                for i in ops do
                    if not type(i, specfunc(KroneckerProduct)) then
                        all_kron := false;
                        break;
                    end if;
                end do;
                if all_kron then
                    left := op(1, ops[1]);
                    right := op(2, ops[1]);
                    for i from 2 to nops(ops) do
                        left := left . op(1, ops[i]);
                        right := right . op(2, ops[i]);
                    end do;
                    to_pair := Kron(left, right);
                    if with_trace then _AddTrace("KRON7 (chain)", p, to_pair); end if;
                    changed := true;
                    return to_pair;
                else
                    return p;
                end if;
            end proc);
        end do;
        return e;
    end proc;
    
    # ---------- Core simplification rules (without trace) ----------
    local ApplyRulesNoTrace := proc(expr)
        local e, e_old, changed;
        local sub, a, b, terms, term, new_sub;
        e := expr;
        changed := true;
        while changed do
            e_old := e;
            e := _DistributeDot(e, false);
            e := _FlattenChain(e, false);
            if e = e_old then changed := false; end if;
        end do;
        
        for sub in indets(e, specfunc(Transpose)) do
            a := op(1,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(Transpose(term), term in terms));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(HermitianTranspose)) do
            a := op(1,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(HermitianTranspose(term), term in terms));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(conjugate)) do
            a := op(1,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(conjugate(term), term in terms));
                e := subs(sub = new_sub, e);
            elif type(a, `.`) then
                new_sub := conjugate(op(1,a)) . conjugate(op(2,a));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(Determinant)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Determinant(op(1,a))^(dim(op(2,a))) * Determinant(op(2,a))^(dim(op(1,a)));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            new_sub := sub;
            if type(a, `*`) and type(op(1,a), numeric) then
                new_sub := op(1,a) * Kron(op(2,a), b);
            elif type(b, `*`) and type(op(1,b), numeric) then
                new_sub := op(1,b) * Kron(a, op(2,b));
            end if;
            if new_sub <> sub then e := subs(sub = new_sub, e); end if;
        end do;
        
        for sub in indets(e, specfunc(Transpose)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(Transpose(op(1,a)), Transpose(op(2,a)));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(HermitianTranspose)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(HermitianTranspose(op(1,a)), HermitianTranspose(op(2,a)));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(conjugate)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(conjugate(op(1,a)), conjugate(op(2,a)));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(op(1,a), Kron(op(2,a), b));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(Kron(term, b), term in terms));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            if type(b, `+`) then
                terms := [op(b)];
                new_sub := `+`(seq(Kron(a, term), term in terms));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(Trace)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Trace(op(1,a)) * Trace(op(2,a));
                e := subs(sub = new_sub, e);
            elif type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(Trace(term), term in terms));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(MatrixInverse)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(MatrixInverse(op(1,a)), MatrixInverse(op(2,a)));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(MatrixInverse)) do
            a := op(1,sub);
            if type(a, specfunc(Shuffle)) then
                new_sub := Shuffle(MatrixInverse(op(1,a)), MatrixInverse(op(2,a)));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(Trace)) do
            a := op(1,sub);
            if type(a, specfunc(Shuffle)) then
                new_sub := Trace(op(1,a)) * Trace(op(2,a));
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        return e;
    end proc;
    
    # ---------- Same with trace (for SimplifyTrace) ----------
    local ApplyRulesWithTrace := proc(expr)
        local e, e_old, changed;
        local sub, a, b, terms, term, new_sub;
        e := expr;
        changed := true;
        while changed do
            e_old := e;
            e := _DistributeDot(e, true);
            e := _FlattenChain(e, true);
            if e = e_old then changed := false; end if;
        end do;
        
        for sub in indets(e, specfunc(Transpose)) do
            a := op(1,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(Transpose(term), term in terms));
                _AddTrace("Transpose over sum", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(HermitianTranspose)) do
            a := op(1,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(HermitianTranspose(term), term in terms));
                _AddTrace("HermitianTranspose over sum", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(conjugate)) do
            a := op(1,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(conjugate(term), term in terms));
                _AddTrace("conjugate over sum", sub, new_sub);
                e := subs(sub = new_sub, e);
            elif type(a, `.`) then
                new_sub := conjugate(op(1,a)) . conjugate(op(2,a));
                _AddTrace("conjugate over product", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(Determinant)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Determinant(op(1,a))^(dim(op(2,a))) * Determinant(op(2,a))^(dim(op(1,a)));
                _AddTrace("KRON9", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            new_sub := sub;
            if type(a, `*`) and type(op(1,a), numeric) then
                new_sub := op(1,a) * Kron(op(2,a), b);
                _AddTrace("KRON1 (scalar left)", sub, new_sub);
            elif type(b, `*`) and type(op(1,b), numeric) then
                new_sub := op(1,b) * Kron(a, op(2,b));
                _AddTrace("KRON1 (scalar right)", sub, new_sub);
            end if;
            if new_sub <> sub then e := subs(sub = new_sub, e); end if;
        end do;
        
        for sub in indets(e, specfunc(Transpose)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(Transpose(op(1,a)), Transpose(op(2,a)));
                _AddTrace("KRON2", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(HermitianTranspose)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(HermitianTranspose(op(1,a)), HermitianTranspose(op(2,a)));
                _AddTrace("KRON3", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(conjugate)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(conjugate(op(1,a)), conjugate(op(2,a)));
                _AddTrace("conjugate over KroneckerProduct", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(op(1,a), Kron(op(2,a), b));
                _AddTrace("KRON4", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            if type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(Kron(term, b), term in terms));
                _AddTrace("KRON5", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(KroneckerProduct)) do
            a := op(1,sub); b := op(2,sub);
            if type(b, `+`) then
                terms := [op(b)];
                new_sub := `+`(seq(Kron(a, term), term in terms));
                _AddTrace("KRON6", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(Trace)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Trace(op(1,a)) * Trace(op(2,a));
                _AddTrace("KRON8", sub, new_sub);
                e := subs(sub = new_sub, e);
            elif type(a, `+`) then
                terms := [op(a)];
                new_sub := `+`(seq(Trace(term), term in terms));
                _AddTrace("Trace linearity", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(MatrixInverse)) do
            a := op(1,sub);
            if type(a, specfunc(KroneckerProduct)) then
                new_sub := Kron(MatrixInverse(op(1,a)), MatrixInverse(op(2,a)));
                _AddTrace("KRON10", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(MatrixInverse)) do
            a := op(1,sub);
            if type(a, specfunc(Shuffle)) then
                new_sub := Shuffle(MatrixInverse(op(1,a)), MatrixInverse(op(2,a)));
                _AddTrace("Inverse of Shuffle", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        for sub in indets(e, specfunc(Trace)) do
            a := op(1,sub);
            if type(a, specfunc(Shuffle)) then
                new_sub := Trace(op(1,a)) * Trace(op(2,a));
                _AddTrace("Trace(Shuffle)", sub, new_sub);
                e := subs(sub = new_sub, e);
            end if;
        end do;
        
        return e;
    end proc;
    
    # ---------- Simplify (without trace) ----------
    Simplify := proc(expr, maxiter::integer := 50)
        local new, old, i;
        new := expr;
        for i from 1 to maxiter do
            old := new;
            new := ApplyRulesNoTrace(new);
            if new = old then break; end if;
        end do;
        return new;
    end proc;
    
    # ---------- SimplifyTrace (step-by-step output) ----------
    SimplifyTrace := proc(expr, maxiter::integer := 50)
        local new, old, i, step, cnt;
        _trace_steps := [];
        new := expr;
        printf("Start expression: %a\n\n", new);
        for i from 1 to maxiter do
            old := new;
            _trace_steps := [];
            new := ApplyRulesWithTrace(new);
            if new = old then
                printf("No more rules apply after iteration %d.\n", i-1);
                break;
            end if;
            printf("Iteration %d:\n", i);
            cnt := 0;
            for step in _trace_steps do
                cnt := cnt + 1;
                printf("  %d. %s\n", cnt, step);
            end do;
            printf("Result after iteration %d: %a\n\n", i, new);
        end do;
        printf("--- Final simplified expression ---\n");
        return new;
    end proc;
    
    # ---------- Expand ----------
    local _ExpandKron := proc(expr)
        local a, b, new_a, new_b, term;
        if not type(expr, specfunc(KroneckerProduct)) then
            return expr;
        end if;
        a := op(1, expr);
        b := op(2, expr);
        new_a := _ExpandKron(a);
        new_b := _ExpandKron(b);
        if type(new_a, `+`) then
            return `+`(seq(Kron(term, new_b), term in [op(new_a)]));
        elif type(new_b, `+`) then
            return `+`(seq(Kron(new_a, term), term in [op(new_b)]));
        else
            return Kron(new_a, new_b);
        end if;
    end proc;
    
    Expand := proc(expr)
        local e, new_e;
        e := expr;
        while true do
            new_e := subsindets(e, specfunc(KroneckerProduct), _ExpandKron);
            if new_e = e then break; end if;
            e := new_e;
        end do;
        return e;
    end proc;
    
    # ---------- Shuffle (KRON11) ----------
    Shuffle := proc(A, B)
        local p, q, r, s, S1, S2;
        try
            p := LinearAlgebra:-RowDimension(A);
            q := LinearAlgebra:-ColumnDimension(A);
            r := LinearAlgebra:-RowDimension(B);
            s := LinearAlgebra:-ColumnDimension(B);
        catch:
            return 'Shuffle'(A, B);
        end try;
        if not (type(p, posint) and type(q, posint) and type(r, posint) and type(s, posint)) then
            return 'Shuffle'(A, B);
        end if;
        S1 := _PerfectShuffleMatrix(p, r);
        S2 := _PerfectShuffleMatrix(q, s);
        return S1 . LinearAlgebra:-KroneckerProduct(A, B) . LinearAlgebra:-Transpose(S2);
    end proc;
    
    # ---------- Activation for numeric evaluation (used in Check/VisualCheck) ----------
    local Activate := proc(expr)
        subs([ 'Determinant' = LinearAlgebra:-Determinant,
               'Trace' = LinearAlgebra:-Trace,
               'MatrixInverse' = LinearAlgebra:-MatrixInverse,
               'KroneckerProduct' = LinearAlgebra:-KroneckerProduct,
               'Transpose' = LinearAlgebra:-Transpose,
               'HermitianTranspose' = LinearAlgebra:-HermitianTranspose,
               dim = LinearAlgebra:-RowDimension
             ], expr);
    end proc;
    
    # ---------- Check (symbolic test, square matrices only) ----------
    Check := proc(expr::uneval, {n::posint := 2, samples::posint := 5})
        printf("Check: Use inert Kron(A,B) or A &x B. Avoid direct LinearAlgebra:-KroneckerProduct.\n");
        printf("   Note: This command assumes all matrix symbols are SQUARE. For non‑square matrices, use Simplify and Verify.\n");
        local e, syms, i, subs_set, gen, simp_expr, val_simp, orig_expr, val_orig, ok;
        e := eval(expr);
        syms := indets(e, symbol) minus {Kron, Determinant, Trace, MatrixInverse, KroneckerProduct,
                                         'Kron', 'Determinant', 'Trace', 'MatrixInverse', 'KroneckerProduct',
                                         dim, 'dim', `&x`, Shuffle, Expand, Transpose, HermitianTranspose, conjugate};
        if syms = {} then error "No symbolic variables found. Use symbolic names like A,B,..."; end if;
        gen := proc() local M; do M := LinearAlgebra:-RandomMatrix(n, n, generator=-5..5); until LinearAlgebra:-Determinant(M) <> 0; return M; end proc;
        ok := true;
        for i from 1 to samples do
            subs_set := map(s -> s = gen(), syms);
            simp_expr := Simplify(e);
            simp_expr := Activate(simp_expr);
            val_simp := eval(eval(simp_expr, subs_set));
            orig_expr := subs('Kron' = KroneckerProduct, `&x` = KroneckerProduct, e);
            orig_expr := Activate(orig_expr);
            val_orig := eval(eval(orig_expr, subs_set));
            if type(val_orig, Matrix) and val_simp = 0 then
                if LinearAlgebra:-Norm(val_orig, infinity) = 0 then
                    # ok
                else
                    ok := false; break;
                end if;
            elif type(val_simp, Matrix) and val_orig = 0 then
                if LinearAlgebra:-Norm(val_simp, infinity) = 0 then
                    # ok
                else
                    ok := false; break;
                end if;
            elif type(val_orig, Matrix) and type(val_simp, Matrix) then
                if not LinearAlgebra:-Equal(val_orig, val_simp) then ok := false; break; end if;
            else
                if evalb(val_orig <> val_simp) then ok := false; break; end if;
            end if;
        end do;
        return ok;
    end proc;
    
    # ---------- VisualCheck (same as Check, shows matrices) ----------
    VisualCheck := proc(expr::uneval, {n::posint := 2})
        printf("VisualCheck: Use inert Kron(A,B) or A &x B. Avoid direct LinearAlgebra:-KroneckerProduct.\n");
        printf("   Note: This command assumes all matrix symbols are SQUARE. For non‑square matrices, use Simplify and Verify.\n");
        local e, syms, subs_set, gen, simp_expr, val_simp, orig_expr, val_orig, a;
        e := eval(expr);
        syms := indets(e, symbol) minus {Kron, Determinant, Trace, MatrixInverse, KroneckerProduct,
                                         'Kron', 'Determinant', 'Trace', 'MatrixInverse', 'KroneckerProduct',
                                         dim, 'dim', `&x`, Shuffle, Expand, Transpose, HermitianTranspose, conjugate};
        if syms = {} then error "No symbolic variables found. Use symbolic names like A,B,..."; end if;
        gen := proc() local M; do M := LinearAlgebra:-RandomMatrix(n, n, generator=-5..5); until LinearAlgebra:-Determinant(M) <> 0; return M; end proc;
        subs_set := map(s -> s = gen(), syms);
        printf("Substitutions used:\n");
        for a in subs_set do printf("%a\n", a); end do;
        printf("\n");
        
        simp_expr := Simplify(e);
        simp_expr := Activate(simp_expr);
        val_simp := eval(eval(simp_expr, subs_set));
        
        orig_expr := subs('Kron' = KroneckerProduct, `&x` = KroneckerProduct, e);
        orig_expr := Activate(orig_expr);
        val_orig := eval(eval(orig_expr, subs_set));
        
        printf("========== Direct calculation ==========\n");
        print(val_orig);
        printf("\n========== Simplified calculation (via module) ==========\n");
        print(val_simp);
        printf("\n");
        
        local identical;
        if type(val_orig, Matrix) and val_simp = 0 then
            identical := LinearAlgebra:-Norm(val_orig, infinity) = 0;
        elif type(val_simp, Matrix) and val_orig = 0 then
            identical := LinearAlgebra:-Norm(val_simp, infinity) = 0;
        elif type(val_orig, Matrix) and type(val_simp, Matrix) then
            identical := LinearAlgebra:-Equal(val_orig, val_simp);
        else
            identical := evalb(val_orig = val_simp);
        end if;
        
        if identical then
            printf("Results are IDENTICAL (visual proof).\n");
            return true;
        else
            printf("ERROR – results differ.\n");
            return false;
        end if;
    end proc;
    
    # ---------- Verify (for concrete matrices, returns boolean) ----------
    Verify := proc(expr)
        printf("Verify: Comparing Simplify(expr) with direct evaluation using concrete matrices.\n");
        local simp, direct;
        simp := Simplify(expr);
        simp := Activate(simp);
        direct := subs(`&x` = LinearAlgebra:-KroneckerProduct, expr);
        direct := Activate(direct);
        if type(simp, Matrix) and type(direct, Matrix) then
            if LinearAlgebra:-Equal(simp, direct) then
                printf("Results are IDENTICAL.\n");
                return true;
            else
                printf("ERROR – results differ.\n");
                return false;
            end if;
        else
            if evalb(simp = direct) then
                printf("Results are IDENTICAL.\n");
                return true;
            else
                printf("ERROR – results differ.\n");
                return false;
            end if;
        end if;
    end proc;
    
    # ---------- VisualVerify (visual comparison for concrete matrices) ----------
    VisualVerify := proc(expr)
        printf("VisualVerify: Comparing Simplify(expr) with direct evaluation using concrete matrices.\n");
        local simp, direct;
        simp := Simplify(expr);
        simp := Activate(simp);
        direct := subs(`&x` = LinearAlgebra:-KroneckerProduct, expr);
        direct := Activate(direct);
        printf("========== Direct calculation ==========\n");
        print(direct);
        printf("\n========== Simplified calculation (via module) ==========\n");
        print(simp);
        printf("\n");
        if type(simp, Matrix) and type(direct, Matrix) then
            if LinearAlgebra:-Equal(simp, direct) then
                printf("Results are IDENTICAL (visual proof).\n");
                return true;
            else
                printf("ERROR – results differ.\n");
                return false;
            end if;
        else
            if evalb(simp = direct) then
                printf("Results are IDENTICAL (visual proof).\n");
                return true;
            else
                printf("ERROR – results differ.\n");
                return false;
            end if;
        end if;
    end proc;
    
    # ---------- NumericEvaluate (robust numeric evaluation for concrete matrices) ----------
    NumericEvaluate := proc(expr)
        local e;
        # First replace all non-Kronecker functions that can appear inside KroneckerProduct arguments
        # Order: innermost first
        e := expr;
        e := subsindets(e, specfunc(MatrixInverse), m -> LinearAlgebra:-MatrixInverse(op(m)));
        e := subsindets(e, specfunc(Transpose), t -> LinearAlgebra:-Transpose(op(t)));
        e := subsindets(e, specfunc(HermitianTranspose), h -> LinearAlgebra:-HermitianTranspose(op(h)));
        e := subsindets(e, specfunc(conjugate), c -> conjugate(op(c)));
        e := subsindets(e, specfunc(Shuffle), s -> LinearAlgebra:-KroneckerProduct(op(2,s), op(1,s)));
        e := subsindets(e, specfunc(dim), d -> LinearAlgebra:-RowDimension(op(d)));
        e := subsindets(e, specfunc(KroneckerProduct), k -> LinearAlgebra:-KroneckerProduct(op(k)));
        e := subsindets(e, specfunc(Determinant), d -> LinearAlgebra:-Determinant(op(d)));
        e := subsindets(e, specfunc(Trace), t -> LinearAlgebra:-Trace(op(t)));
        # Now evaluate the whole expression (including sums)
        return eval(e);
    end proc;
    
    # ---------- Info (educational) ----------
    Info := proc()
        printf("================================================================\n");
        printf("KRONECKERRULES PACKAGE - COMPREHENSIVE EDUCATIONAL GUIDE\n");
        printf("================================================================\n\n");
        printf("1. PURPOSE\n");
        printf("   Simplify symbolic expressions involving Kronecker products,\n");
        printf("   matrix multiplication, determinants, traces, inverses, transposes,\n");
        printf("   conjugates, and shuffling.\n\n");
        printf("2. INERT NOTATION (to avoid immediate evaluation)\n");
        printf("   Use either:   Kron(A,B)   or   A &x B\n");
        printf("   Both produce an inert 'KroneckerProduct(A,B)' that the module can simplify.\n");
        printf("   Avoid using LinearAlgebra:-KroneckerProduct directly in expressions\n");
        printf("   you want to simplify symbolically.\n\n");
        printf("3. MAIN COMMANDS\n");
        printf("   Simplify(expr)  - apply all algebraic rules (KRON1‑KRON11) iteratively.\n");
        printf("   Expand(expr)    - fully distribute Kronecker products over sums.\n");
        printf("   Shuffle(A,B)    - implement perfect shuffle: B ⊗ A = S (A⊗B) S^T.\n");
        printf("   Check(expr)     - test symbolic identities using random square matrices.\n");
        printf("   VisualCheck(expr) - same as Check but shows random matrices and results.\n");
        printf("   Verify(expr)    - compare Simplify(expr) with direct evaluation for user‑supplied concrete matrices.\n");
        printf("   VisualVerify(expr) - same as Verify but displays both results (for any matrix dimensions).\n");
        printf("   NumericEvaluate(expr) - directly evaluate concrete expressions to numeric matrices (robust).\n");
        printf("   SimplifyTrace(expr) - same as Simplify but prints each applied rule (with markers).\n");
        printf("   Info()          - this help text.\n\n");
        printf("4. CHECK AND VISUALCHECK – VERIFY SYMBOLIC IDENTITIES (SQUARE MATRICES ONLY)\n");
        printf("   These commands generate random square matrices (size n x n) for each symbol.\n");
        printf("   They are useful for testing algebraic identities, but assume all matrices are square.\n");
        printf("   For non‑square matrices, use Simplify and test with your own concrete matrices via Verify or VisualVerify.\n\n");
        printf("5. NUMERIC EVALUATION (FOR CONCRETE MATRICES)\n");
        printf("   NumericEvaluate(expr) is the recommended way to compute a numeric result from an expression\n");
        printf("   built with &x, Shuffle, Transpose, etc., after you have assigned concrete matrices.\n");
        printf("   It replaces all inert functions in the correct order and returns the final matrix or scalar.\n\n");
        printf("6. VERIFY AND VISUALVERIFY – FOR CONCRETE MATRICES (ANY DIMENSIONS)\n");
        printf("   Verify(expr) returns true/false after comparing simplified and direct results.\n");
        printf("   VisualVerify(expr) shows the direct and simplified results for visual inspection.\n");
        printf("   Use these after you have defined your matrices (e.g., A := Matrix(...)).\n\n");
        printf("7. SIMPLIFYTRACE – EDUCATIONAL STEP-BY-STEP OUTPUT\n");
        printf("   SimplifyTrace(expr) works like Simplify but prints each rule application.\n");
        printf("   The changing subexpression is marked with ** **. Iterations and steps are numbered.\n");
        printf("   This helps you understand how the simplification proceeds.\n\n");
        printf("8. MATRIX DIMENSION REQUIREMENTS\n");
        printf("   - Kronecker product works for any matrices (even non‑square).\n");
        printf("   - Ordinary multiplication (A . B), determinant, trace, inverse require square matrices.\n");
        printf("   - For determinant rule, dim(X) returns the row dimension (assumed square).\n");
        printf("   - Shuffle requires known numeric dimensions; otherwise returns inert.\n\n");
        printf("9. LIST OF SIMPLIFICATION RULES (KRON1‑KRON11)\n");
        printf("   KRON1: (cA)⊗B = A⊗(cB) = c (A⊗B)\n");
        printf("   KRON2: (A⊗B)^T = A^T ⊗ B^T\n");
        printf("   KRON3: (A⊗B)^H = A^H ⊗ B^H\n");
        printf("   KRON4: (A⊗B)⊗C = A⊗(B⊗C)\n");
        printf("   KRON5: (A+B)⊗C = A⊗C + B⊗C\n");
        printf("   KRON6: A⊗(B+C) = A⊗B + A⊗C\n");
        printf("   KRON7: (A⊗B)(C⊗D) = (A·C)⊗(B·D)\n");
        printf("   KRON8: trace(A⊗B) = trace(A)·trace(B), linear over sums\n");
        printf("   KRON9: det(A⊗B) = det(A)^(dim B) · det(B)^(dim A)\n");
        printf("   KRON10: (A⊗B)^(-1) = A^(-1) ⊗ B^(-1)\n");
        printf("   KRON11: Shuffle(A,B) = S (A⊗B) S^T gives B⊗A.\n");
        printf("   Extra: transpose, conjugate, Hermitian transpose distribute over sums and products.\n\n");
        printf("10. NOTES AND LIMITATIONS\n");
        printf("   - No automatic dimension checking; ensure compatibility.\n");
        printf("   - Shuffle requires known dimensions.\n");
        printf("   - Avoid using LinearAlgebra:-KroneckerProduct directly.\n");
        printf("   - Check and VisualCheck are only for symbolic square matrices; for concrete or non‑square, use NumericEvaluate or Verify/VisualVerify.\n");
        printf("================================================================\n");
    end proc;
    
end module:

with(KroneckerRules):
Info();
   
MAIN COMMANDS
   Simplify(expr)  - apply all algebraic rules (KRON1‑KRON11) iteratively.
   Expand(expr)    - fully distribute Kronecker products over sums.
   Shuffle(A,B)    - implement perfect shuffle: B ⊗ A = S (A⊗B) S^T.
   Check(expr)     - test symbolic identities using random square matrices.
   VisualCheck(expr) - same as Check but shows random matrices and results.
   Verify(expr)    - compare Simplify(expr) with direct evaluation for user‑supplied concrete matrices.
   VisualVerify(expr) - same as Verify but displays both results (for any matrix dimensions).
   NumericEvaluate(expr) - directly evaluate concrete expressions to numeric matrices (robust).
   SimplifyTrace(expr) - same as Simplify but prints each applied rule (with markers).
   Info()          - this help text.

delete

pde_vrgstuk_Can_pdsolve_handle_a_procedure_in_the_pde_mprimesDEF_18-4-2026.mw

Download weideproblee_nu_voor_mpromesDEF_13-4-2026.mw

restart;
with(plots):
with(Optimization):

# Parameters
u := 0.3:
v := -0.1:
N := 3:

# Trial function: Fourier series up to 3rd harmonic
x_trial := a0 + a1*cos(t) + b1*sin(t) + a2*cos(2*t) + b2*sin(2*t) + a3*cos(3*t) + b3*sin(3*t):

# Lagrangian
L_expr := 1/2*diff(x_trial, t)^2 + u^2*cos(x_trial) + x_trial*v*sin(t):

# Objective function: numerically integrate the Lagrangian over t = 0..2π
obj := proc(C::Vector)
    local integrand, val;
    integrand := eval(L_expr, [a0=C[1], a1=C[2], b1=C[3], a2=C[4], b2=C[5], a3=C[6], b3=C[7]]);
    val := evalf(Int(integrand, t=0..2*Pi, method=_d01ajc));
    return val;
end proc:

# Initial guess (slightly perturbed to avoid zero-gradient warning)
init := Vector([0.3, 0.1, 0.001, 0.001, -0.001, 0.001, -0.001], datatype=float[8]):

# Minimize the action using the SQP method
opt := NLPSolve(7, obj, initialpoint=init, method=sqp):
optimal_coeffs := opt[2];
J_min := opt[1];

printf("Optimal coefficients:\n");
printf("a0 = %.6f\n", optimal_coeffs[1]);
printf("a1 = %.6f, b1 = %.6f\n", optimal_coeffs[2], optimal_coeffs[3]);
printf("a2 = %.6f, b2 = %.6f\n", optimal_coeffs[4], optimal_coeffs[5]);
printf("a3 = %.6f, b3 = %.6f\n", optimal_coeffs[6], optimal_coeffs[7]);

# Ritz solution function
x_Ritz := t -> optimal_coeffs[1] + optimal_coeffs[2]*cos(t) + optimal_coeffs[3]*sin(t)
                + optimal_coeffs[4]*cos(2*t) + optimal_coeffs[5]*sin(2*t)
                + optimal_coeffs[6]*cos(3*t) + optimal_coeffs[7]*sin(3*t):

# Numerical reference solution from the Euler-Lagrange ODE
ode := diff(x(t), t$2) = -u^2*sin(x(t)) + v*sin(t):
ics := x(0)=0.3, D(x)(0)=0.1:
sol_numeric := dsolve({ode, ics}, numeric, range=0..2*Pi):

# Plot comparison
p1 := odeplot(sol_numeric, [t, x(t)], 0..2*Pi, color=blue, thickness=2, legend="numerical (dsolve)");
p2 := plot(x_Ritz(t), t=0..2*Pi, color=red, thickness=2, linestyle=dash, legend="Ritz (N=3)");
display(p1, p2, title="Comparison of solutions", labels=["t", "x(t)"]);

Download integral_en_minimum_berekend_mprimesDEF_9-4-2026.mw

This explanation does not cover all derivations; it simply gives an idea of how a rotation matrix is derived

Download rotatienmatrices_mprimes31-3-2026.mw

@Mapleliquid 
Hello , yes Heineken a fouvorite beer of  mine :-)

There is a complete guide  in Maple, so thi s make sense :
 

A Complete Guide for performing Tensors computations using Physics via maple help

A Complete Guide for Tensor computations using Physics - MaplePrimes

Tensors in Maple – definition and index actions (PDF)

with(Physics):
Setup();

@salim-barzani 

restart;
with(plots):

# Definieer de functie voor de 3D plot met parameters en componentkeuze
energie_3d_keuze := proc(comp, m, p, q, r, theta, w, v, t, alpha1, alpha2, alpha3, alpha4)
    local Delta, xi, Knum, Pnum, Hnum, plot_expressie, plot_kleur, plot_titel;
    # Bereken Delta
    Delta := sqrt(4*alpha1*w^2 + (4*theta*alpha3 - 4)*w + 3*alpha4*m^2 + 4*theta^2*alpha2 - 4*r) / 
             (2*sqrt(p^2*alpha1 + p*q*alpha3 + q^2*alpha2));
    # Definieer xi
    xi := p*x + q*y - t*v;
    # Bereken energie componenten
    Knum := - (m^2/(16*p^2*alpha1 + 16*p*q*alpha3 + 16*q^2*alpha2)) * 
            (3*m^2*alpha4 + 4*theta^2*alpha2 + 4*theta*w*alpha3 + 4*w^2*alpha1 - 4*r - 4*w) * 
            (-1 + cos(2*Delta*xi));
     Pnum := (1/(4*p^2*alpha1 + 4*p*q*alpha3 + 4*q^2*alpha2)) * cos(Delta*xi)^2 * 
            (alpha4*m^2*cos(Delta*xi)^2 + 2*alpha1*w^2 + (2*theta*alpha3 - 2)*w + 2*theta^2*alpha2 - 2*r) * m^2;
     Hnum := (m^2/(8*p^2*alpha1 + 8*p*q*alpha3 + 8*q^2*alpha2)) * 
            (2*alpha4*m^2*cos(Delta*xi)^4 - 3*alpha4*m^2*cos(Delta*xi)^2 + 
             3*alpha4*m^2 + 4*theta*w*alpha3 + 4*alpha1*w^2 + 4*theta^2*alpha2 - 4*r - 4*w);
    
    # Kies welke component geplot wordt
    if comp = 1 then
        plot_expressie := Knum;
        plot_kleur := red;
        plot_titel := "Kinetische Energie (Knum)";
    elif comp = 2 then
        plot_expressie := Pnum;
        plot_kleur := blue;
        plot_titel := "Potentiële Energie (Pnum)";
    else
        plot_expressie := Hnum;
        plot_kleur := green;
        plot_titel := "Totale Energie (Hnum)";
    end if;
    
    # Maak de 3D plot
    plot3d(plot_expressie, x = -5..5, y = -5..5, 
           view = [-5..5, -5..5, 0..0.35], color = plot_kleur,
           title = sprintf("%s (m=%.2f, p=%.2f, q=%.2f, t=%.2f)", plot_titel, m, p, q, t),
           labels = ["x", "y", "energie"],
           style = patchcontour, axes = boxed);
end proc:

# Maak de Explore versie met componentkeuze (zonder controller optie)
Explore(energie_3d_keuze(comp, m, p, q, r, theta, w, v, t, alpha1, alpha2, alpha3, alpha4),      
        parameters = [
            comp = [1, 2, 3],
            m = 0.1..5.0, p = 0.1..5.0, q = 0.1..5.0,
            r = 0.1..5.0, theta = 0.1..5.0, w = 0.1..5.0,
            v = 0.1..5.0, t = 0.1..5.0,
            alpha1 = 0.1..5.0, alpha2 = 0.1..5.0,
            alpha3 = 0.1..5.0, alpha4 = 0.1..5.0
        ],      
        initialvalues = [
            comp = 1, m = 1, p = 1, q = 1, r = 1, theta = 1, w = 1, v = 1, t = 1,
            alpha1 = 1, alpha2 = 1, alpha3 = 1, alpha4 = 1
        ],    
        placement = right);

Maybe this gui ?

Beta cannot be 0 , as a quick test, not symbolic , but numerically 

# Choose random numerical values for the parameters
#params := {beta = 0.5, eta1 = 1.2, eta2 = 0.8, gamma1 = 0.3, gamma2 = -0.4, N = 2};

# Evaluate both expressions with these values
eval(G0, params);
eval(G01, params);

# Check if they are (approximately) equal
evalb(eval(G0 - G01, params) = 0);
# Or use a tolerance
is(abs(eval(G0 - G01, params)) < 1e-10);

Maybe useful?

nanofluid_figures_mprimesDEFA_19-2-2026.mw

t := 5*Pi/9:
L := [-tan(4*Pi/9) + 4*sin(4*Pi/9)];
evalf~(L,20)[];

identify(-1.7320508075688772935);

A more practical approach ,well done..Maple, for recognizing a root form from the decimal form.

oneliner..( a long line ) 

toPrefix := e -> `if`(type(e,atomic), e,[op(0,e), seq(toPrefix(op(i,e)), i=1..nops(e))]):

ex := (a+b*c)^2:
toPrefix(ex);


exUltra :=
subs(a=b+c,
    expand(
        int(
            diff((x+y)^5*sin(exp(x^2+y^2)),x)
            /(1+sum(k*z^k,k=1..4)),
            x=0..1
        )
    )
):

toPrefix(exUltra );

exInsane :=
((f@g)@@3)(a+b*c)
+
(exp@@2)(sin@@3(x+y))
+
((a&*b)&+c)^5;

 exInsane := f(g(f(g(f(g(b c + a)))))) + @@(exp, 2)(@@(sin, 3))

                   5
    + (a &* b) &+ c 

toPrefix(exInsane );

A name is  a variable when there is a valued assigned to it. 

Used seq command here instead  a loop construction
Note: how fast the trees are growing :-) 

>      >  restart; with(Fractals[LSystem]): with(plots): # Define the L-system generator LGenerator := proc(n, Init, Angle0, Param, Rules, Color)     local Iteration;          # Use try-catch to avoid errors     try         if n = 0 then             Iteration := Init;         else             Iteration := Iterate(Init, Rules, n);         end if;     catch:         # If there's an error, use the initial string         Iteration := Init;     end try;          # Draw the L-system     LSystemPlot(         Iteration,         Param,         initialangle = Angle0,         color = Color,         scaling = constrained,         axes = none,         thickness = 2     ); end proc: # First L-system (number 17) Init17 := "X": Angle17 := 90: Param17 := [   "F" = "draw:1",   "+" = "turn:-25.7",   "-" = "turn:25.7",   "[" = "push:position;push:angle",   "]" = "pop:position;pop:angle" ]: Rules17 := [   "F" = "FF",   "X" = "F[+X][-X]FX" ]: # Second L-system (number 18) Init18 := "X": Angle18 := 90: Param18 := [   "F" = "draw:1",   "+" = "turn:20",   "-" = "turn:-20",   "[" = "push:position;push:angle",   "]" = "pop:position;pop:angle" ]: Rules18 := [   "F" = "FF",   "X" = "F[+X]F[-X]+X" ]: # Test both L-systems without animation print("L-system 17, iteration 4:"); display(LGenerator(4, Init17, Angle17, Param17, Rules17, "Red")); print("L-system 18, iteration 4:"); display(LGenerator(4, Init18, Angle18, Param18, Rules18, "DarkGreen")); # Create smoother animations with more frames print("Smooth animation of L-system 17 (iterations 0-7):"); display([seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)],         insequence = true,         frames = 64);  # More frames for smoother animation print("Smooth animation of L-system 18 (iterations 0-7):"); display([seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)],         insequence = true,         frames = 64);  # More frames for smoother animation # Optional: Show all iterations side by side in a grid print("All iterations of L-system 17 (0-7):"); display(Matrix(2, 4, [seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)])); print("All iterations of L-system 18 (0-7):"); display(Matrix(2, 4, [seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)])); # Even smoother version with interpolation (more iterations shown) print("Extra smooth animation of L-system 17 (more frames):"); display([seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)],         insequence = true,         frames = 128);  # Even more frames for extra smoothness print("Extra smooth animation of L-system 18 (more frames):"); display([seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)],         insequence = true,         frames = 128); >

Download fractal_animatie_mprimes_26-1-2026.mw