janhardo

More rules added ..kroneckerproductmodul...

Answered: janhardo 910

21 hours ago

0 0

Simplified commands

@Harry Garst Ja, en hopelijk kan je...

Answered: janhardo 910

Yesterday at 10:40 AM

0 0

@Harry Garst
Ja , en hopelijk kan je er iets mee met de module code :-)
Nog wat verder uitbreiden met nieuwe regels
groet Jan

27-4-2026 ...

Answered: janhardo 910

Yesterday at 8:12 AM

0 0

27-4-2026

>

27-6-2026

>

restart;

>

KroneckerRules := module()
    option package;
    export Kron, `&x`, Simplify, Expand, Shuffle, Check, VisualCheck, Verify, VisualVerify, NumericEvaluate, Info, SimplifyTrace, ToggleKronNotation;

    # ---------- Inert notation ----------
    Kron := (A,B) -> 'KroneckerProduct'(A,B):
    `&x` := (A,B) -> 'KroneckerProduct'(A,B):

    # ---------- Trace mechanism (with improved rule descriptions) ----------
    local _trace_steps := [];

    # Table with mathematical descriptions of each rule
    local _rule_desc := table([
        "KRON7 (chain)" = "(A5B)·(C5D) → (A·C)5(B·D)",
        "Dot left-distribution" = "X·(Y+Z) → X·Y + X·Z",
        "Dot right-distribution" = "(X+Y)·Z → X·Z + Y·Z",
        "KRON5" = "(X+Y)5Z → X5Z + Y5Z",
        "KRON6" = "X5(Y+Z) → X5Y + X5Z",
        "KRON1 (scalar left)" = "(c·A)5B → c·(A5B)",
        "KRON1 (scalar right)" = "A5(c·B) → c·(A5B)",
        "KRON2" = "(A5B)^T → A^T 5 B^T",
        "KRON3" = "(A5B)^H → A^H 5 B^H",
        "KRON4" = "(A5B)5C → A5(B5C)",
        "KRON8" = "Tr(A5B) → Tr(A)·Tr(B)",
        "KRON9" = "det(A5B) → det(A)^(dim B) · det(B)^(dim A)",
        "KRON10" = "(A5B)^(-1) → A^(-1) 5 B^(-1)",
        "Transpose over sum" = "(A+B)^T → A^T + B^T",
        "HermitianTranspose over sum" = "(A+B)^H → A^H + B^H",
        "conjugate over sum" = "conjugate(A+B) → conjugate(A)+conjugate(B)",
        "conjugate over product" = "conjugate(A·B) → conjugate(A)·conjugate(B)",
        "Trace linearity" = "Tr(A+B) → Tr(A)+Tr(B)",
        "Inverse of Shuffle" = "Shuffle(A,B)^(-1) → Shuffle(A^(-1), B^(-1))",
        "Trace(Shuffle)" = "Tr(Shuffle(A,B)) → Tr(A)·Tr(B)"
    ]);

    local _AddTrace := proc(rule, from_expr, to_expr)
    local desc;
    desc := _rule_desc[rule];
    if desc = NULL then desc := rule; end if;
    _trace_steps := [op(_trace_steps),
        sprintf("%s (%s): %a -> %a", rule, desc, 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, outputstyle::string := "default")
    local new, old, i, format_expr;
    new := expr;
    for i from 1 to maxiter do
        old := new;
        new := ApplyRulesNoTrace(new);
        if new = old then break; end if;
    end do;
    format_expr := proc(e)
        if outputstyle = "Kron" then
            return subsindets(e, specfunc(KroneckerProduct),
                proc(k) 'Kron'(op(k)); end proc);
        else
            return e;
        end if;
    end proc;
    return format_expr(new);
end proc;


    # ---------- SimplifyTrace (step-by-step output with optional outputstyle) ----------
    SimplifyTrace := proc(expr, maxiter::integer := 50, outputstyle::string := "default")
    local new, old, i, step, cnt, format_expr;
    _trace_steps := [];
    new := expr;

    format_expr := proc(e)
        if outputstyle = "Kron" then
            return subsindets(e, specfunc(KroneckerProduct),
                proc(k) 'Kron'(op(k)); end proc);
        else
            return e;
        end if;
    end proc;

    printf("Start expression: %a\n\n", format_expr(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;
            if outputstyle = "Kron" then
                step := StringTools:-SubstituteAll(step, "KroneckerProduct", "Kron");
            end if;
            printf(" %d. %s\n", cnt, step);
        end do;
        printf("Result after iteration %d: %a\n\n", i, format_expr(new));
    end do;
    printf("--- Final simplified expression ---\n");
    return ifelse(outputstyle = "Kron", format_expr(new), new);
end proc;

    # ---------- Toggle GLOBAL output notation of KroneckerProduct ----------
     ToggleKronNotation := proc()
    global `print/KroneckerProduct`;
    if assigned(`print/KroneckerProduct`) then
        unassign(`print/KroneckerProduct`);
        printf("KroneckerProduct output notation set to DEFAULT (KroneckerProduct(A,B))\n");
    else
        `print/KroneckerProduct` := (A,B) -> 'Kron'(A,B);
        printf("KroneckerProduct output notation set to PRETTY (Kron(A,B))\n");
    end if;
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;
        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)));
        return eval(e);
    end proc;

    # ---------- Info (educational) ----------
    Info := proc()
    printf("================================================================\n");
    printf("KRONECKERRULES PACKAGE – USER GUIDE\n");
    printf("================================================================\n\n");

    printf("1. PURPOSE\n");
    printf("   Simplify symbolic expressions involving Kronecker products (5),\n");
    printf("   matrix multiplication (.), addition (+), subtraction (-), transpose (^T),\n");
    printf("   Hermitian transpose (^H), inverse, determinant, trace, and shuffle.\n\n");

    printf("2. INERT NOTATION (prevents immediate evaluation)\n");
    printf("   - Kron(A,B)   or   A &x B\n");
    printf("   Both produce an inert 'KroneckerProduct(A,B)' for symbolic simplification.\n");
    printf("   Avoid using LinearAlgebra:-KroneckerProduct directly.\n\n");

    printf("3. MAIN COMMANDS\n");
    printf("   Simplify(expr)                     – apply all rules iteratively (KroneckerProduct output).\n");
    printf("   Simplify(expr, \"Kron\")            – same, but output uses short Kron notation.\n");
    printf("   Simplify(expr, maxiter)            – set maximum number of iterations (default 50).\n");
    printf("   Simplify(expr, maxiter, \"Kron\")   – both maxiter and short notation.\n\n");

    printf("   SimplifyTrace(expr)                – show steps (KroneckerProduct notation).\n");
    printf("   SimplifyTrace(expr, \"Kron\")       – show steps using short Kron notation.\n");
    printf("   SimplifyTrace(expr, maxiter)       – set iterations.\n");
    printf("   SimplifyTrace(expr, maxiter, \"Kron\") – both.\n\n");

    printf("   Expand(expr)                       – fully distribute 5 over sums.\n");
    printf("   Shuffle(A,B)                       – realize B 5 A = S·(A5B)·S^T (perfect shuffle).\n");
    printf("   ToggleKronNotation()               – switch global output between Kron and KroneckerProduct.\n\n");

    printf("4. VERIFICATION OF IDENTITIES\n");
    printf("   Check(expr)                        – test numerically (random SQUARE matrices).\n");
    printf("   VisualCheck(expr)                  – same, but shows matrices and results.\n");
    printf("   Verify(expr)                       – compare Simplify(expr) with direct evaluation.\n");
    printf("       * uses USER-supplied concrete matrices (any dimensions).\n");
    printf("   VisualVerify(expr)                 – same, displays both computations.\n\n");

    printf("5. NUMERIC EVALUATION\n");
    printf("   NumericEvaluate(expr)              – evaluate concrete expressions (inert -> active).\n\n");

    printf("6. STEP-BY-STEP EXPLANATION – SimplifyTrace\n");
    printf("   Each iteration scans the whole expression.\n");
    printf("   All applicable rules are applied (order is not guaranteed).\n");
    printf("   Output: rule name + mathematical description + old -> new subexpression.\n");
    printf("   Example: KRON7 (chain) ((A5B)·(C5D) → (A·C)5(B·D)): ...\n\n");

    printf("7. SUMMARY OF RULES (KRON1–KRON11)\n");
    printf("   KRON1: (cA)5B = A5(cB) = c·(A5B)\n");
    printf("   KRON2: (A5B)^T = A^T 5 B^T\n");
    printf("   KRON3: (A5B)^H = A^H 5 B^H\n");
    printf("   KRON4: (A5B)5C = A5(B5C)\n");
    printf("   KRON5: (A+B)5C = A5C + B5C\n");
    printf("   KRON6: A5(B+C) = A5B + A5C\n");
    printf("   KRON7: (A5B)·(C5D) = (A·C)5(B·D)\n");
    printf("   KRON8: Tr(A5B) = Tr(A)·Tr(B)\n");
    printf("   KRON9: det(A5B) = det(A)^{dim(B)} · det(B)^{dim(A)}\n");
    printf("   KRON10: (A5B)^{-1} = A^{-1} 5 B^{-1}\n");
    printf("   KRON11: Shuffle(A,B) = S·(A5B)·S^T gives B5A\n");
    printf("   Extra: dot distributes over sums, trace linearity, etc.\n\n");

    printf("8. DIMENSIONS AND LIMITATIONS\n");
    printf("   - Kronecker product works for any matrices (even rectangular).\n");
    printf("   - Ordinary multiplication (.), determinant, trace, inverse require square matrices.\n");
    printf("   - Shuffle needs exact dimensions; otherwise it remains inert.\n");
    printf("   - Check and VisualCheck assume all matrix symbols are square.\n");
    printf("   - Use Verify/VisualVerify for non‑square or user‑provided matrices.\n\n");

    printf("9. EXAMPLES (copy and paste into Maple)\n");
    printf("   with(KroneckerRules):\n");
    printf("   expr := (A &x B) . (C &x D);\n");
    printf("   SimplifyTrace(expr);\n");
    printf("   SimplifyTrace(expr, \"Kron\");\n");
    printf("   expr2 := (2*A &x (B+C)) . (Transpose(D) &x MatrixInverse(E));\n");
    printf("   Simplify(expr2, \"Kron\");\n");
    printf("   A := Matrix(2,2,[[1,2],[3,4]]); B := Matrix(2,2,[[0,1],[1,0]]);\n");
    printf("   Verify( (A &x B) . (B &x A) );\n\n");

    printf("10. ADDITIONAL INFORMATION\n");
    printf("    - Use 'ToggleKronNotation()' to change global output style.\n");
    printf("    - Adjust 'maxiter' for deep associativity (default 50).\n");
    printf("    - Shuffle is only numerically evaluable when all dimensions are known.\n");
    printf("================================================================\n");
end proc;

end module:

>

>	with(KroneckerRules);

(1.1)

>	ToggleKronNotation();

KroneckerProduct output notation set to PRETTY (Kron(A,B))

>	expr := (A &x B) . (C &x d).(E &x F); # 3 kronecker products

(1.2)

>	expr := Kron(A,B) . Kron(C,d) . Kron(E,F);

(1.3)

>	expr :=((KroneckerProduct(A, B)) . (KroneckerProduct(C, d))) . (KroneckerProduct(E, F));

(1.4)

>	expr :=((KroneckerProduct(A, B)) . (KroneckerProduct(C, d))) . (KroneckerProduct(E, F))=Simplify(expr);

(1.5)

>	expr :=((KroneckerProduct(A, B)) . (KroneckerProduct(C, d))) . (KroneckerProduct(E, F))=Simplify(expr,"Kron");

(1.6)

>

>	Simplify(expr);

(1.7)

>	Simplify(expr, "Kron");

(1.8)

>	SimplifyTrace(expr);# start over with only expr := (A &x B) . (C &x d).(E &x F); # 3 kronecker products

Start expression: KroneckerProduct(A,B) . KroneckerProduct(C,d) . KroneckerProduct(E,F)

Iteration 1:
1. KRON7 (chain) ((A5B)·(C5D) → (A·C)5(B·D)): KroneckerProduct(A,B) . KroneckerProduct(C,d) . KroneckerProduct(E,F) -> KroneckerProduct(A . C . E,B . d . F)
Result after iteration 1: KroneckerProduct(A . C . E,B . d . F)

No more rules apply after iteration 1.
--- Final simplified expression ---

(1.9)

>	SimplifyTrace(expr, "Kron");

Start expression: Kron(A,B) . Kron(C,d) . Kron(E,F)

Iteration 1:
1. KRON7 (chain) ((A5B)·(C5D) → (A·C)5(B·D)): Kron(A,B) . Kron(C,d) . Kron(E,F) -> Kron(A . C . E,B . d . F)
Result after iteration 1: Kron(A . C . E,B . d . F)

No more rules apply after iteration 1.
--- Final simplified expression ---

(1.10)

>	#Info(); # remove #

>

Download kroneckerproducts_module__27-4-2026.mw

@janhardo Testing commando's&nb...

Answered: janhardo 910

April 25 2026

0 0

@janhardo
Testing commando's
Performance test.comparison, seems to be quick with the module

with(LinearAlgebra):

# matrix size (adjust as needed!)
n := 40:

# random matrices
A := RandomMatrix(n,n, generator=-5..5):
B := RandomMatrix(n,n, generator=-5..5):
C := RandomMatrix(n,n, generator=-5..5):
d := RandomMatrix(n,n, generator=-5..5):

# ---------- 1. Naive method ----------
t1 := time():

K1 := KroneckerProduct(A,B):
K2 := KroneckerProduct(C,d):
R1 := K1 . K2:

t_naive := time() - t1:

# ---------- 2. Smart method (KRON7) ----------
t2 := time():

R2 := KroneckerProduct(A . C, B . d):

t_smart := time() - t2:

# ---------- Results ----------
printf("Naive time  = %.4f sec\n", t_naive);
printf("Smart time  = %.4f sec\n", t_smart);
printf("Speedup(faster)     = %.2f x\n", t_naive / t_smart);

# ---------- Verification ----------
printf("Equal? %a\n", Equal(R1, R2));
Naive time  = 244.7970 sec
Smart time  = 1.2500 sec
Speedup(faster)     = 195.84 x
Equal? true

> ....

Answered: janhardo 910

April 22 2026

0 0

>

restart;
with(LinearAlgebra):
ZehfussReplace := proc(expr::uneval)
    uses LinearAlgebra:
    local replace;
    replace := proc(det)
        local kr_args, A, B;
        kr_args := op(1, det);
        if not type(kr_args, specfunc(KroneckerProduct)) then return eval(det) end if;
        A := op(1, kr_args);
        B := op(2, kr_args);
        'Determinant'(A)^(RowDimension(B)) * 'Determinant'(B)^(RowDimension(A));
    end proc;
    subsindets(expr, specfunc(specfunc(KroneckerProduct), Determinant), replace);
end proc:

A := Matrix(2,2,symbol=a);
B := Matrix(3,3,symbol=b);

# Important: call the procedure directly with the unevaluated form
result := ZehfussReplace( Determinant(KroneckerProduct(A,B)) + 2*Determinant(KroneckerProduct(A,A)) );
# print(result);
  eval(result);

(a[1, 1]*a[2, 2]-a[1, 2]*a[2, 1])^3*(b[1, 1]*b[2, 2]*b[3, 3]-b[1, 1]*b[2, 3]*b[3, 2]-b[1, 2]*b[2, 1]*b[3, 3]+b[1, 2]*b[2, 3]*b[3, 1]+b[1, 3]*b[2, 1]*b[3, 2]-b[1, 3]*b[2, 2]*b[3, 1])^2+2*(a[1, 1]*a[2, 2]-a[1, 2]*a[2, 1])^4

(1)

Download Zehfuss_-kronecker_determinant__mprimes_22-4-2026.mw

@Harry Garst It doesn't s...

Answered: janhardo 910

April 21 2026

0 0

@Harry Garst

It doesn't seem necessary to use two variables.

These are two ways of describing the same variable..

Apparently, the expression hold for all n =1.. 2
proof by induction further

kronecker_-Zehfuss_mprimes_20-4-2026.mw...

Answered: janhardo 910

April 20 2026

0 0

kronecker_-Zehfuss_mprimes_20-4-2026.mw

For a bounded surface, would that ...

Answered: janhardo 910

April 20 2026

0 0

Of all closed surfaces in space with a fixed surface area A, which surface encloses the largest volume?

But there are two variables in the functional here.

a good example to to do by hand...

> ....

Answered: janhardo 910

April 20 2026

0 0

Via differential geometry

>

restart;
Typesetting:-Settings(typesetdot=true);

printf("================================================================\n");
printf("ISOPERIMETRIC PROBLEM: Maximize area for fixed perimeter\n");
printf("================================================================\n\n");
printf("Goal: Prove that the optimal closed curve is a circle.\n");
printf("We do this by showing that its curvature must be constant.\n\n");

# Step 1: Lagrangian
F := x(t)*diff(y(t),t) + lambda*sqrt(diff(x(t),t)^2 + diff(y(t),t)^2);
printf("Step 1: Lagrangian F = %a\n\n", F);

# Step 2: Euler-Lagrange equations
printf("Step 2: Euler-Lagrange equations\n");
EL1 := diff(Physics:-diff(F, diff(x(t),t)), t) - Physics:-diff(F, x(t)) = 0;
EL2 := diff(Physics:-diff(F, diff(y(t),t)), t) - Physics:-diff(F, y(t)) = 0;
printf("EL1: %a\n", EL1);
printf("EL2: %a\n\n", EL2);

# Step 3: First integrals (integrate once)
printf("Step 3: First integrals\n");
eq1 := lambda*diff(x(t),t)/sqrt(diff(x(t),t)^2+diff(y(t),t)^2) = y(t) + C1;
eq2 := x(t) + lambda*diff(y(t),t)/sqrt(diff(x(t),t)^2+diff(y(t),t)^2) = C2;
printf("First integral I: %a\n", eq1);
printf("First integral II: %a\n\n", eq2);

# Step 4: Switch to arc length parameter s
printf("Step 4: Arc length parametrization (ds/dt = sqrt(x'^2+y'^2))\n");
printf("Then dx/ds = x'/(ds/dt), dy/ds = y'/(ds/dt)\n");
printf("The first integrals become:\n");
printf("   lambda * dx/ds = y + C1      (A)\n");
printf("   lambda * dy/ds = C2 - x      (B)\n\n");

# Step 5: Compute curvature explicitly
printf("Step 5: Curvature calculation\n");
printf("For a plane curve, curvature kappa = (dx/ds)*(d^2y/ds^2) - (dy/ds)*(d^2x/ds^2)\n");
printf("From (A) and (B) we obtain:\n");

dx_ds := (y + C1)/lambda;
dy_ds := (C2 - x)/lambda;
printf("   dx/ds = %a\n", dx_ds);
printf("   dy/ds = %a\n", dy_ds);

d2x_ds2 := (C2 - x)/lambda^2;
d2y_ds2 := -(y + C1)/lambda^2;
printf("   d^2x/ds^2 = %a\n", d2x_ds2);
printf("   d^2y/ds^2 = %a\n\n", d2y_ds2);

kappa := dx_ds * d2y_ds2 - dy_ds * d2x_ds2;
kappa_expanded := expand(kappa);
printf("Substituting into the curvature formula and expanding:\n");
printf("   kappa = %a\n", kappa_expanded);
printf("Thus kappa = -(C2-x)^2/lambda^3 - (y+C1)^2/lambda^3\n\n");

# Step 6: Use the arc length normalization
printf("Step 6: Arc length normalization\n");
printf("Since (dx/ds)^2 + (dy/ds)^2 = 1, from (A) and (B):\n");
printf("   ((y+C1)/lambda)^2 + ((C2-x)/lambda)^2 = 1\n");
printf("   => (y+C1)^2 + (C2-x)^2 = lambda^2\n");
norm_eq := (y+C1)^2 + (C2-x)^2 = lambda^2;
printf("Hence: %a\n\n", norm_eq);

# Compact simplification of kappa_final
printf("Step 6a: Compact simplification of curvature\n");
# Substitute the sum directly
kappa_final := algsubs((y+C1)^2+(C2-x)^2 = lambda^2, kappa_expanded);
printf("   kappa = %a\n", kappa_final);
printf("Therefore the curvature is constant: kappa = -1/lambda.\n");
printf("This proves that the optimal curve has constant curvature.\n\n");

# Step 7: Conclusion – circle
printf("Step 7: Conclusion\n");
printf("A plane curve with constant nonzero curvature is an arc of a circle.\n");
printf("Since the curve is closed and smooth, it must be a full circle.\n");
printf("The radius is R = |lambda|, and the center is (C2, -C1).\n");
printf("Thus the optimal shape is a circle.\n\n");

# Step 8: Example plot
printf("Step 8: Example plot (lambda=1, C1=0, C2=0)\n");
lambda_val := 1; C1_val := 0; C2_val := 0; phi := 0;
x_circ := t -> C2_val + lambda_val*cos(2*Pi*t + phi);
y_circ := t -> -C1_val + lambda_val*sin(2*Pi*t + phi);
plot([x_circ(t), y_circ(t), t=0..1], scaling=constrained,
     title="Optimal closed curve: circle", labels=["x","y"],
     color=blue, thickness=2);
printf("Plot displayed.\n");

================================================================
ISOPERIMETRIC PROBLEM: Maximize area for fixed perimeter
================================================================

Goal: Prove that the optimal closed curve is a circle.
We do this by showing that its curvature must be constant.

Step 1: Lagrangian F = x(t)*diff(y(t),t)+lambda*(diff(x(t),t)^2+diff(y(t),t)^2)^(1/2)

Step 2: Euler-Lagrange equations

-(1/2)*lambda*(diff(x(t), t))*(2*(diff(x(t), t))*(diff(diff(x(t), t), t))+2*(diff(y(t), t))*(diff(diff(y(t), t), t)))/((diff(x(t), t))^2+(diff(y(t), t))^2)^(3/2)+lambda*(diff(diff(x(t), t), t))/((diff(x(t), t))^2+(diff(y(t), t))^2)^(1/2)-(diff(y(t), t)) = 0

diff(x(t), t)-(1/2)*lambda*(diff(y(t), t))*(2*(diff(x(t), t))*(diff(diff(x(t), t), t))+2*(diff(y(t), t))*(diff(diff(y(t), t), t)))/((diff(x(t), t))^2+(diff(y(t), t))^2)^(3/2)+lambda*(diff(diff(y(t), t), t))/((diff(x(t), t))^2+(diff(y(t), t))^2)^(1/2) = 0

EL1: -1/2*lambda/(diff(x(t),t)^2+diff(y(t),t)^2)^(3/2)*diff(x(t),t)*(2*diff(x(t),t)*diff(diff(x(t),t),t)+2*diff(y(t),t)*diff(diff(y(t),t),t))+lambda/(diff(x(t),t)^2+diff(y(t),t)^2)^(1/2)*diff(diff(x(t),t),t)-diff(y(t),t) = 0
EL2: diff(x(t),t)-1/2*lambda/(diff(x(t),t)^2+diff(y(t),t)^2)^(3/2)*diff(y(t),t)*(2*diff(x(t),t)*diff(diff(x(t),t),t)+2*diff(y(t),t)*diff(diff(y(t),t),t))+lambda/(diff(x(t),t)^2+diff(y(t),t)^2)^(1/2)*diff(diff(y(t),t),t) = 0

Step 3: First integrals

lambda*(diff(x(t), t))/((diff(x(t), t))^2+(diff(y(t), t))^2)^(1/2) = y(t)+C1

x(t)+lambda*(diff(y(t), t))/((diff(x(t), t))^2+(diff(y(t), t))^2)^(1/2) = C2

First integral I: lambda/(diff(x(t),t)^2+diff(y(t),t)^2)^(1/2)*diff(x(t),t) = y(t)+C1
First integral II: x(t)+lambda/(diff(x(t),t)^2+diff(y(t),t)^2)^(1/2)*diff(y(t),t) = C2

Step 4: Arc length parametrization (ds/dt = sqrt(x'^2+y'^2))
Then dx/ds = x'/(ds/dt), dy/ds = y'/(ds/dt)
The first integrals become:
lambda * dx/ds = y + C1 (A)
lambda * dy/ds = C2 - x (B)

Step 5: Curvature calculation
For a plane curve, curvature kappa = (dx/ds)*(d^2y/ds^2) - (dy/ds)*(d^2x/ds^2)
From (A) and (B) we obtain:

dx/ds = (y+C1)/lambda
dy/ds = (C2-x)/lambda

(C2-x)/lambda^2

-(y+C1)/lambda^2

d^2x/ds^2 = (C2-x)/lambda^2
d^2y/ds^2 = -(y+C1)/lambda^2

-(C2-x)^2/lambda^3-(y+C1)^2/lambda^3

-C2^2/lambda^3+2*C2*x/lambda^3-x^2/lambda^3-C1^2/lambda^3-2*y*C1/lambda^3-y^2/lambda^3

Substituting into the curvature formula and expanding:
   kappa = -1/lambda^3*C2^2+2/lambda^3*C2*x-1/lambda^3*x^2-1/lambda^3*C1^2-2/lambda^3*y*C1-1/lambda^3*y^2
Thus kappa = -(C2-x)^2/lambda^3 - (y+C1)^2/lambda^3

Step 6: Arc length normalization
Since (dx/ds)^2 + (dy/ds)^2 = 1, from (A) and (B):
   ((y+C1)/lambda)^2 + ((C2-x)/lambda)^2 = 1
   => (y+C1)^2 + (C2-x)^2 = lambda^2

Hence: (y+C1)^2+(C2-x)^2 = lambda^2

Step 6a: Compact simplification of curvature

kappa = -1/lambda
Therefore the curvature is constant: kappa = -1/lambda.
This proves that the optimal curve has constant curvature.

Step 7: Conclusion
A plane curve with constant nonzero curvature is an arc of a circle.
Since the curve is closed and smooth, it must be a full circle.
The radius is R = |lambda|, and the center is (C2, -C1).
Thus the optimal shape is a circle.

Step 8: Example plot (lambda=1, C1=0, C2=0)

Plot displayed.

>

via vectorcalculus
bewijs_cirkel_gesloten_cvorm_viavectot_calculus_mprimes_20-4-2026.mw

Download bewijs_cirkel_gesloten_cvorm_via_diff_geom_mprimes_20-4-2026.mw

@Rouben Rostamian Thanks. Ye...

Answered: janhardo 910

April 20 2026

0 0

@Rouben Rostamian

Thanks for your educational code , yes, this circle is a interesting example . I had actually seen another derivation involving the curvature of the circle, but that leads into vector calculus and differential geometry, which require more advanced knowledge.

It’s also instructive to see different approaches to solving the problem.

Of all closed curves with a fixed length...

Answered: janhardo 910

April 19 2026

0 0

Of all closed curves with a fixed length L, which curve encloses the largest area?

> ....

Answered: janhardo 910

April 18 2026

0 0

>

restart;
with(plots):

# ============================================================
# NONLINEAR WAVE EQUATION SOLVED AS COUPLED ODEs
# ============================================================
#
# ORIGINAL PDE: ∂²u/∂t² = c²(λ(H + ∂u/∂x)) · ∂²u/∂x²
#
# AFTER SPATIAL DISCRETIZATION (FINITE DIFFERENCES):
# We obtain a SYSTEM of COUPLED ORDINARY DIFFERENTIAL EQUATIONS (ODEs):
#
# For each internal point i = 1, 2, ..., Nx-1:
#
#   d²u_i/dt² = c²(λ(H + (u_{i+1}-u_{i-1})/(2Δx))) · (u_{i+1}-2u_i+u_{i-1})/Δx²
#
# This is a system of (Nx-1) coupled, nonlinear ODEs.
# ============================================================

# ============================================================
# STEP 1: PHYSICAL SYSTEM PARAMETERS
# ============================================================
L := 0.04;                         # Length of the rod [m] (4 cm)
A0 := 12*L/2*10^(-6);              # Initial amplitude [m] (~2.4e-7 m)
Nx := 100;                         # Number of spatial points (creates 99 ODEs)
dx := L/Nx;                        # Spatial step size [m]
x_vals := [seq(i*dx, i=0..Nx)];    # Positions along the rod [m]

printf("========================================\n");
printf("     COUPLED ODEs FROM WAVE EQUATION\n");
printf("========================================\n\n");
printf("STEP 1: SYSTEM PARAMETERS\n");
printf(" Rod length L = %.4f m (%.1f cm)\n", L, L*100);
printf(" Initial amplitude A0 = %.2e m\n", A0);
printf(" Number of spatial points Nx+1 = %d\n", Nx+1);
printf(" Number of COUPLED ODEs = %d (for internal points)\n", Nx-1);
printf(" Spatial step dx = %.2e m\n\n", dx);

# ============================================================
# STEP 2: MATERIAL MODEL (PIECEWISE FUNCTIONS)
# ============================================================
# c_sq(H) = squared speed of sound [m²/s²]
c_sq := proc(H)
    if H < 700 then
        return -3143*H + 1.99e7;      # Region 1: decreasing speed
    elif H < 1000 then
        return 7333*H + 1.257e7;      # Region 2: increasing speed
    else
        return 1.99e7;                # Region 3: constant speed
    end if;
end proc:

# lambda(H) = helper function for inverse relation
lambda := proc(H)
    if H < 1000 then
        return (3/250000000)*H;
    else
        return 3/250000;
    end if;
end proc:

printf("STEP 2: MATERIAL MODEL\n");
printf(" c²(H) = piecewise function:\n");
printf("    H < 700:   c² = -3143·H + 1.99e7\n");
printf("    700 ≤ H < 1000: c² = 7333·H + 1.257e7\n");
printf("    H ≥ 1000: c² = 1.99e7 (constant)\n");
printf(" λ(H) = 3/250000000 * H for H < 1000\n\n");

# ============================================================
# STEP 3: NUMERICAL PARAMETERS
# ============================================================
c_max := sqrt(1.99e7);              # Maximum wave speed [m/s]
dt := 0.5 * dx / c_max;             # Time step (CFL = 0.5 for stability)
f0 := sqrt(c_sq(650))/(2*L);        # Resonance frequency [Hz]
T0 := 1/f0;                         # Period [s]
Nt := 1500;                         # Number of time steps

printf("STEP 3: NUMERICAL PARAMETERS\n");
printf(" Maximum wave speed c_max = %.1f m/s\n", c_max);
printf(" Resonance frequency f0 = %.1f Hz (ultrasonic!)\n", f0);
printf(" Period T0 = %.2e s\n", T0);
printf(" Time step dt = %.2e s\n", dt);
printf(" Number of time steps Nt = %d\n", Nt);
printf(" CFL number = %.3f (must be < 1 for stability)\n\n", c_max*dt/dx);

# ============================================================
# STEP 4: INITIALIZE THE SOLUTION (Three time layers)
# ============================================================
# We need three arrays for the time integration:
#   u_old = solution at time t - dt
#   u_curr = solution at time t
#   u_new = solution at time t + dt

u_old := Array(0..Nx);
u_curr := Array(0..Nx);
u_new := Array(0..Nx);

# INITIAL CONDITION: u(x,0) = A0·cos(πx/L)
# This is the FIRST HARMONIC vibration mode
for i from 0 to Nx do
    u_curr[i] := A0 * cos(Pi * x_vals[i+1] / L);
    u_old[i] := u_curr[i];         # Initial velocity = 0 => u_old = u_curr
end do:

printf("STEP 4: INITIAL CONDITION\n");
printf(" u(x,0) = A0·cos(πx/L)\n");
printf(" ∂u/∂t(x,0) = 0\n");
printf(" This creates the first harmonic vibration mode.\n");
printf(" The rod is maximally stretched at the middle (x = %.2f m).\n\n", L/2);

# Display the initial values for the first few points (educational)
printf(" Initial values for the first 5 points (i=0 to 4):\n");
for i from 0 to 4 do
    printf("    u[%d](0) = %.2e m at x = %.4f m\n", i, u_curr[i], x_vals[i+1]);
end do;
printf("    ...\n\n");

# ============================================================
# STEP 5: BOUNDARY CONDITIONS (Neumann - Free Ends)
# ============================================================
# Neumann: ∂u/∂x = 0 at both ends
# This means: u[0] = u[1] and u[Nx] = u[Nx-1]
# Physically: The ends are FREE (no force)
apply_bc := proc()
    u_new[0] := u_new[1];
    u_new[Nx] := u_new[Nx-1];
end proc:

printf("STEP 5: BOUNDARY CONDITIONS\n");
printf(" Neumann: ∂u/∂x = 0 at x = 0 and x = L\n");
printf(" This means: u₀(t) = u₁(t) and u_%d(t) = u_%d(t)\n", Nx, Nx-1);
printf(" Physically: The ends are FREE (can move without force).\n");
printf(" Waves reflect at the ends WITHOUT sign reversal.\n\n");

# ============================================================
# STEP 6: THE COUPLED ODE SYSTEM - EXPLANATION
# ============================================================
printf("STEP 6: THE COUPLED ODE SYSTEM\n");
printf(" After spatial discretization, we get %d coupled ODEs:\n\n", Nx-1);

for i from 1 to 3 do
    printf(" ODE %d: d²u_%d/dt² = c²(λ(H + (u_%d-u_%d)/(2dx))) · (u_%d-2u_%d+u_%d)/dx²\n",
           i, i, i+1, i-1, i+1, i, i-1);
end do;
printf(" ...\n");
for i from Nx-3 to Nx-1 do
    printf(" ODE %d: d²u_%d/dt² = c²(λ(H + (u_%d-u_%d)/(2dx))) · (u_%d-2u_%d+u_%d)/dx²\n",
           i, i, i+1, i-1, i+1, i, i-1);
end do;
printf("\n Each ODE connects point i with its neighbors i-1 and i+1.\n");
printf(" This is why they are called COUPLED ODEs.\n\n");

# ============================================================
# STEP 7: SOLVING THE COUPLED ODEs USING FINITE DIFFERENCES
# ============================================================
# We use an explicit time integration scheme:
#   u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)
#
# This comes from the central difference approximation:
#   d²u_i/dt² ≈ (u_new[i] - 2*u_curr[i] + u_old[i]) / dt²

printf("STEP 7: SOLVING THE COUPLED ODEs\n");
printf(" Using explicit time integration (central differences):\n");
printf("    u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)\n\n");
printf(" Starting numerical integration...\n\n");

# Storage indices for visualization
save_indices := [42, 106, 211, 317, 422];
solutions := table();

# TIME INTEGRATION LOOP
for n from 1 to Nt do

    # --------------------------------------------------------
    # COMPUTE ALL COUPLED ODEs
    # Each iteration of this loop computes ONE ODE for point i
    # The ODEs are COUPLED because u_new[i] depends on u_curr[i-1] and u_curr[i+1]
    # --------------------------------------------------------
    for i from 1 to Nx-1 do

        # ----- COUPLING TERM 1: LOCAL STRAIN (∂u/∂x) -----
        # This connects u_i with its neighbors u_{i-1} and u_{i+1}
        du_dx := (u_curr[i+1] - u_curr[i-1]) / (2 * dx);

        # ----- MATERIAL PARAMETER H DEPENDS ON STRAIN -----
        H_total := 650 + du_dx;

        # ----- SQUARED SPEED OF SOUND c² -----
        if H_total < 700 then
            c2 := -3143*H_total + 1.99e7;
        elif H_total < 1000 then
            c2 := 7333*H_total + 1.257e7;
        else
            c2 := 1.99e7;
        end if;

        # ----- COUPLING TERM 2: SECOND DERIVATIVE (∂²u/∂x²) -----
        # This also connects u_i with its neighbors u_{i-1} and u_{i+1}
        u_xx := (u_curr[i+1] - 2*u_curr[i] + u_curr[i-1]) / dx^2;

        # ----- THE COUPLED ODE FOR POINT i -----
        # d²u_i/dt² = c² · u_xx
        # Discretized in time:
        u_new[i] := 2*u_curr[i] - u_old[i] + dt^2 * c2 * u_xx;

    end do;
    # --------------------------------------------------------
    # END OF COUPLED ODE COMPUTATION
    # --------------------------------------------------------

    # Apply boundary conditions (coupling the end points)
    apply_bc();

    # Store solutions for visualization
    for idx from 1 to 5 do
        if n = save_indices[idx] then
            solutions[n] := copy(u_new);
            printf(" Stored: t = %.3f T0 (time step %d/%d)\n", n*dt/T0, n, Nt);
        end if;
    end do;

    # Shift time layers for next iteration
    for i from 0 to Nx do
        u_old[i] := u_curr[i];
        u_curr[i] := u_new[i];
    end do;

    # Progress indicator
    if n mod 500 = 0 then
        printf(" Progress: %d/%d (%.0f%%)\n", n, Nt, 100*n/Nt);
    end if;
end do:

printf("\nSimulation completed!\n\n");

# ============================================================
# STEP 8: VISUALIZATION OF THE COUPLED ODE SOLUTION
# ============================================================
printf("STEP 8: VISUALIZATION\n");
printf(" Generating plots of the solution u_i(t)...\n");

colors := ["red", "blue", "green", "magenta", "black"];
plot_list := [];

for idx from 1 to 5 do
    n_idx := save_indices[idx];
    t_val := evalf(n_idx * dt / T0, 3);
    points := [seq([x_vals[i+1], solutions[n_idx][i]], i=0..Nx)];
    p := plot(points,
              labels=["x [m]", "u(x,t) [m]"],
              title=sprintf("Wave at t = %.3f T₀", t_val),
              color=colors[idx],
              legend=sprintf("t = %.3f T₀", t_val),
              thickness=2);
    plot_list := [op(plot_list), p];
end do:

display(plot_list, title="Solution of the Coupled ODE System", size=[800,500]);
printf(" Plot shows the wave at 5 different times.\n\n");

# ============================================================
# STEP 9: ANIMATION OF THE COUPLED ODE SOLUTION
# ============================================================
printf("STEP 9: ANIMATION\n");
printf(" Generating animation of the moving wave...\n");

animation_frames := [];
t_labels := [0.1, 0.25, 0.5, 0.75, 1.0];
explanation := [
    "t = 0.10 T₀: Wave begins moving right.\nThe wave deforms because c² depends on strain.",
    "t = 0.25 T₀: Wave approaches the right end.\nAmplitude changes due to nonlinearity.",
    "t = 0.50 T₀: Wave reaches the end and reflects.\n∂u/∂x=0 at the end (free boundary).",
    "t = 0.75 T₀: Wave moves back left.\nWave shape is mirrored.",
    "t = 1.00 T₀: Wave returns to initial condition.\nOne complete period completed."
];

for idx from 1 to 5 do
    n_idx := save_indices[idx];
    t_val := evalf(n_idx * dt / T0, 3);
    points := [seq([x_vals[i+1], solutions[n_idx][i]], i=0..Nx)];

    p := plot(points,
              labels=["x [m]", "u(x,t) [m]"],
              title=sprintf("Longitudinal wave at t = %.2f T₀", t_labels[idx]),
              color="blue",
              thickness=3,
              view=[0..L, -A0*2..A0*2]);

    text := plots:-textplot([0.02, A0*1.5, explanation[idx]],
                            font=["Arial", 10],
                            color="darkred",
                            align="LEFT");

    animation_frames := [op(animation_frames), display(p, text)];
end do:

# Repeat frames for smooth animation
smooth_frames := [];
for rep from 1 to 4 do
    for idx from 1 to 5 do
        smooth_frames := [op(smooth_frames), animation_frames[idx]];
    end do;
end do:

display(smooth_frames, insequence=true, size=[800,600]);
printf(" Animation started! (%d frames)\n\n", nops(smooth_frames));

# ============================================================
# STEP 10: SUMMARY FOR THE STUDENT
# ============================================================
printf("========================================\n");
printf("              SUMMARY\n");
printf("========================================\n\n");
printf("WHAT IS A COUPLED ODE SYSTEM?\n");
printf(" • We have %d unknown functions u_i(t), i = 1..%d\n", Nx-1, Nx-1);
printf(" • Each ODE connects u_i with its neighbors u_{i-1} and u_{i+1}\n");
printf(" • This is why they are called COUPLED ODEs\n\n");
printf("WHERE ARE THE ODEs IN THE CODE?\n");
printf(" • Inside the main time loop (for n from 1 to Nt)\n");
printf(" • In the inner loop (for i from 1 to Nx-1)\n");
printf(" • Each iteration computes ONE ODE for point i\n\n");
printf("WHAT DOES EACH ODE DO?\n");
printf(" • Computes d²u_i/dt² = c² · (u_{i+1}-2u_i+u_{i-1})/dx²\n");
printf(" • Updates u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)\n\n");
printf("ANSWER TO THE ORIGINAL QUESTION:\n");
printf(" Q: Can pdsolve handle c²(λ(H + ∂u/∂x))?\n");
printf(" A: NO, pdsolve CANNOT handle this nonlinear PDE.\n");
printf(" SOLUTION: Convert to COUPLED ODEs and use finite differences.\n");
printf("========================================\n");

========================================
     COUPLED ODEs FROM WAVE EQUATION
========================================

STEP 1: SYSTEM PARAMETERS
  Rod length L = 0.0400 m (4.0 cm)
  Initial amplitude A0 = 2.40e-07 m
  Number of spatial points Nx+1 = 101
  Number of COUPLED ODEs = 99 (for internal points)
  Spatial step dx = 4.00e-04 m

STEP 2: MATERIAL MODEL
  c²(H) = piecewise function:
    H < 700:   c² = -3143·H + 1.99e7
    700 ≤ H < 1000: c² = 7333·H + 1.257e7
    H ≥ 1000: c² = 1.99e7 (constant)
  λ(H) = 3/250000000 * H for H < 1000

STEP 3: NUMERICAL PARAMETERS
  Maximum wave speed c_max = 4460.9 m/s
  Resonance frequency f0 = 52822.0 Hz (ultrasonic!)
  Period T0 = 1.89e-05 s
  Time step dt = 4.48e-08 s
  Number of time steps Nt = 1500
  CFL number = 0.500 (must be < 1 for stability)

STEP 4: INITIAL CONDITION
  u(x,0) = A0·cos(πx/L)
  ∂u/∂t(x,0) = 0
  This creates the first harmonic vibration mode.
  The rod is maximally stretched at the middle (x = 0.02 m).

  Initial values for the first 5 points (i=0 to 4):
    u[0](0) = 2.40e-07 m at x = 0.0000 m
    u[1](0) = 2.40e-07 m at x = 0.0004 m
    u[2](0) = 2.40e-07 m at x = 0.0008 m
    u[3](0) = 2.39e-07 m at x = 0.0012 m
    u[4](0) = 2.38e-07 m at x = 0.0016 m
    ...

Warning, (in apply_bc) `u_new` is implicitly declared local

STEP 5: BOUNDARY CONDITIONS
  Neumann: ∂u/∂x = 0 at x = 0 and x = L
  This means: u₀(t) = u₁(t) and u_100(t) = u_99(t)
  Physically: The ends are FREE (can move without force).
  Waves reflect at the ends WITHOUT sign reversal.

STEP 6: THE COUPLED ODE SYSTEM
  After spatial discretization, we get 99 coupled ODEs:

  ODE 1: d²u_1/dt² = c²(λ(H + (u_2-u_0)/(2dx))) · (u_2-2u_1+u_0)/dx²
  ODE 2: d²u_2/dt² = c²(λ(H + (u_3-u_1)/(2dx))) · (u_3-2u_2+u_1)/dx²
  ODE 3: d²u_3/dt² = c²(λ(H + (u_4-u_2)/(2dx))) · (u_4-2u_3+u_2)/dx²
  ...
  ODE 97: d²u_97/dt² = c²(λ(H + (u_98-u_96)/(2dx))) · (u_98-2u_97+u_96)/dx²
  ODE 98: d²u_98/dt² = c²(λ(H + (u_99-u_97)/(2dx))) · (u_99-2u_98+u_97)/dx²
  ODE 99: d²u_99/dt² = c²(λ(H + (u_100-u_98)/(2dx))) · (u_100-2u_99+u_98)/dx²

  Each ODE connects point i with its neighbors i-1 and i+1.
  This is why they are called COUPLED ODEs.

STEP 7: SOLVING THE COUPLED ODEs
  Using explicit time integration (central differences):
    u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)

  Starting numerical integration...

  Stored: t = 0.099 T0 (time step 42/1500)
  Stored: t = 0.251 T0 (time step 106/1500)
  Stored: t = 0.500 T0 (time step 211/1500)
  Stored: t = 0.751 T0 (time step 317/1500)
  Stored: t = 0.999 T0 (time step 422/1500)
  Progress: 500/1500 (33%)
  Progress: 1000/1500 (67%)
  Progress: 1500/1500 (100%)

Simulation completed!

STEP 8: VISUALIZATION
  Generating plots of the solution u_i(t)...

Plot shows the wave at 5 different times.

STEP 9: ANIMATION
Generating animation of the moving wave...

  Animation started! (20 frames)

========================================
              SUMMARY
========================================

WHAT IS A COUPLED ODE SYSTEM?
  • We have 99 unknown functions u_i(t), i = 1..99
  • Each ODE connects u_i with its neighbors u_{i-1} and u_{i+1}
  • This is why they are called COUPLED ODEs

WHERE ARE THE ODEs IN THE CODE?
  • Inside the main time loop (for n from 1 to Nt)
  • In the inner loop (for i from 1 to Nx-1)
  • Each iteration computes ONE ODE for point i

WHAT DOES EACH ODE DO?
  • Computes d²u_i/dt² = c² · (u_{i+1}-2u_i+u_{i-1})/dx²
  • Updates u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)

ANSWER TO THE ORIGINAL QUESTION:
  Q: Can pdsolve handle c²(λ(H + ∂u/∂x))?
  A: NO, pdsolve CANNOT handle this nonlinear PDE.
  SOLUTION: Convert to COUPLED ODEs and use finite differences.
========================================

>

Download pde_vrgstuk_Can_pdsolve_handle_a_procedure_in_the_pde_mprimesDEFA_18-4-2026.mw

@Rouben Rostamian Thanks for the c...

Answered: janhardo 910

April 17 2026

0 0

@Rouben Rostamian
Thanks for the code! It showed me a different way to use the Euler-Lagrange differential equation. for his first derivative

@Rouben Rostamian Well done.Howeve...

Answered: janhardo 910

April 16 2026

0 0

@Rouben Rostamian

Amzing work too,

However, I get the impression that it’s not quite as straightforward as in dharr’s code.

Take this as an example:

To further simplify, observe that the fraction (diff(y(x), x))/sqrt(1+(diff(y(x), x))^2) takes values between -1 and 1, and therefore

we may take it to be the sine of some angle theta, as in sin(theta) = (diff(y(x), x))/sqrt(1+(diff(y(x), x))^2), `and`(-(1/2)*Pi <= theta, theta <= (1/2)*Pi), is this trivial ?

@dharr Amazing work. One thi...

Answered: janhardo 910

April 14 2026

0 0

@dharr
Amazing work.
One thing about the fence length : 100 m ?

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

910 Reputation

13 Badges

MaplePrimes Activity

These are replies submitted by janhardo

More rules added ..kroneckerproductmodul...

@Harry Garst Ja, en hopelijk kan je...

27-4-2026 ...

@janhardo Testing commando's&nb...

> ....

@Harry Garst It doesn't s...

kronecker_-Zehfuss_mprimes_20-4-2026.mw...

For a bounded surface, would that ...

> ....

@Rouben Rostamian Thanks. Ye...

Of all closed curves with a fixed length...

> ....

@Rouben Rostamian Thanks for the c...

@Rouben Rostamian Well done.Howeve...

@dharr Amazing work. One thi...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again