Maple Transactions Volume 4 Issue 4 has now been published.

This issue has two Featured Contributions by people who have been plenary speakers at Maple Conferences in the past, namely Veselin Jungić and Juana Sendra. We hope you enjoy both articles.  There is an accompanying video by Professor Sendra, which we will add a link to when it becomes ready.

As usual, there is an article in the Editor's Corner, but this one is a bit different.  In this one, Michelle Hatzell (the new copyeditor for Maple Transactions, who is also a Masters' student working with me at Western) and I have written about a fun use of Maple's colour contour plots to make an image that might be used as the cover of an upcoming book, namely Perturbation Methods using backward error, which I'm just finishing now with Nic Fillion and which SIAM will publish next year.  So, while there's some math in that paper, it's more about Maple's utilities for colour plotting; so you might find it useful.  We also hope you like at least some of the images.  Some are more attractive than others!

We have several Refereed Contributions, not all of which are ready at this time of publication but which will be added as they are revised and sent in.  We have a nice paper on using continued fractions in a high school context, another on code generation, and another on using Digital Signal Processing in Engineering courses.

Finally we have a first publication in French, by Jalale Soussi.  Actually we have the paper also in English: we chose to publish both, in our Communications section, each with links to the other.  It is possible to publish in Maple Transactions solely in French, of course, but the author provided both, so why not?

Happy reading, and best wishes for 2025. 

Download dynamics_of_novel_disease_outbreak.mw

I am very happy to announce the first public release of a project which I have been working on for the last couple of years.

NODEMaple consists of a set of Maple workbooks and a library for structural design based on the Eurocode.

Currently the main development of the workbooks is focused on "Eurocode 5: Design of timber structures" with the Norwegian Annex.

This software has been made public in the hope of that it might be useful for other structural designers, professionals as well as students. Everyone interested is very Welcome to contribute to this project. The code is published under the GPLv3 license.

For more information see https://github.com/Anthrazit68/NODEMaple.

I didn't put it in the title, but of course this is a post about Advent of Code, in particular Days 16 and 18 which feature a perenial favorite type of problem: finding shortest paths in mazes.

Your input for these is always a maze given as an ascii map.  Like so:

###############
#.......#....E#
#.#.###.#.###.#
#.....#.#...#.#
#.###.#####.#.#
#.#....#....#.#
#.#.#####.###.#
#...........#.#
###.#.#####.#.#
#...#.....#.#.#
#.#.#.###.#.#.#
#.....#...#.#.#
#.###.#.#.#.#.#
#S......#.#...#
###############

There's lots of ways to import one of these into Maple and then solve the maze, but I am to highlight how to do it with GraphTheory.  I am going to start with a GridGraph and then remove the walls in order to leave a just the vertices that represent the paths:

with(StringTools): with(GraphTheory):
maze:=
"###############
#.......#....E#
#.#.###.#.###.#
#.....#.#...#.#
#.###.#####.#.#
#.#....#....#.#
#.#.#####.###.#
#...........#.#
###.#.#####.#.#
#...#.....#.#.#
#.#.#.###.#.#.#
#.....#...#.#.#
#.###.#.#.#.#.#
#S......#.#...#
###############
":
mazelines := (Split(Trim(maze), "\n")):
sgrid := ListTools:-Reverse((map(Explode, mazelines)) ):
m,n := nops(sgrid), nops(sgrid[1]);
tgrid := table([seq(seq([i,j]=sgrid[i,j],i=1..m),j=1..n)]):
start := lhs(select(e->rhs(e)="S", [entries(tgrid,'pairs')])[]);
finish := lhs(select(e->rhs(e)="E", [entries(tgrid,'pairs')])[]);

Now the maze is stored in the table tgrid, and it is easy to find the walls and paths.  In a GridGraph the vertices are labeled with their coordinates as "x,y" and so we rewrite our list of paths in that form, so we can create the induced subgraph of the Grid that includes only those vertices.

(walls,paths) := selectremove(e->rhs(e)="#", [entries(tgrid, 'pairs')]):
paths := map(s->sprintf("%d,%d",lhs(s)[]), paths):
H := SpecialGraphs:-GridGraph(m,n);
G := InducedSubgraph(H, paths);

We can use StyleVertex to highlight the start and finish.

StyleVertex(G, sprintf("%d,%d",start[]), color="LimeGreen");
StyleVertex(G, sprintf("%d,%d",finish[]), color="Red");

plots:-display(<
DrawGraph(H, stylesheet=[vertexshape="square", vertexborder=false, vertexcolor="Black"], showlabels=false) | 
DrawGraph(G, stylesheet=[vertexshape="square", vertexborder=false, vertexcolor="Black"], showlabels=false)>);

(I omitted a step where I set the vertex locations of the maze grid, you can see that in the attached worksheet)

Now finding a path through the maze is as easy as calling GraphTheory:-ShortestPath

sp := ShortestPath(G, sprintf("%d,%d",start[]), sprintf("%d,%d",finish[]) ):

StyleVertex(G, sp[2..-2], color="Orange");
StyleEdge(G, [seq({sp[i],sp[i+1]}, i=1..nops(sp)-1)], color="Orange");
DrawGraph(G, stylesheet=[vertexshape="square", vertexpadding=10, vertexborder=false,
             vertexcolor="Black"],  showlabels=false, size=[800,800]);

Now, Advent of Code seldom gives you a completely simple maze like this, often these is a twist like having to calculate the costs of turns seperately from the cost of steps, or each direction or position has a seperate cost associated with it.

For example, Day 16 has us starting facing east, and then turns cost 1000, while moving forward costs 1. That sort of problem is no longer exactly a maze, instead of the vertices being representing an "x,y" position, instead you increase the number of vertices by a factor of 4, so that you have a vertex for every position and orientation "x,y,o" with edges of weight 1 between adjacent vertices with the same orientation and edges of wieght 1000 to connect "x,y,N" to "x,y,E" and "x,y,W" e.g.  In that sort of weighted graph, we can use GraphTheory:-DijkstrasAlgorithm to find the shortest path and it's weighted cost.

In this code, we expand our list of maze locations with directions, and the use the grid table to generate a list of weighted edges:

dtable := table([0=[0,1], 1=[1,0], 2=[0,-1], 3=[-1,0]]):
dname := table([0="N",1="E",2="S",3="W"]):
dpaths := map(s->local d;seq(cat(s,",",d), d in ["N","E","S","W"]), paths):

edges := NULL:
for i from 1 to m do for j from 1 to n do
    if tgrid[[i,j]] = "#" then next; end if;
    for d from 0 to 3 do
        dir := dtable[d];
        if tgrid[[i,j]+dir] <> "#" then
            edges := edges, [{cat("",i,",",j,",",dname[d]), cat("",i+dir[1],",",j+dir[2],",",dname[d])},1];
        end if;
        edges := edges, [{cat("",i,",",j,",",dname[d]), cat("",i,",",j,",",dname[d+1 mod 4])}, 1000],
                 [{cat("",i,",",j,",",dname[d]), cat("",i,",",j,",",dname[d-1 mod 4])}, 1000];
    end do;
end do; end do:

Gd := Graph(dpaths,weighted,{edges});

Once that is done, it's a simple matter of calling Dijkstra's Algorithm on the graph, but notice that we can reach the finsh while traveling north or east, so we need to find the sortest path to both (you can pass a list of vertices to Dijkstra, and it will efficiently calculate paths to all of them), and select the smaller of the two:

spds := DijkstrasAlgorithm(Gd, cat("",start[1],",",start[2],",E"), 
    [cat("",finish[1],",",finish[2],",N"), cat("",finish[1],",",finish[2],",E")] , 
    distance):
i := min[index](map2(op,2,spds)):
spd := spds[i];

spd := [["2,2,E", "3,2,E", "4,2,E", "4,2,N", "4,3,N", "4,4,N", "4,5,N", "4,6,N", "4,6,E", "5,6,E",
 "6,6,E", "7,6,E", "8,6,E", "8,6,N", "8,7,N", "8,8,N", "8,9,N", "8,10,N", "8,11,N", "8,12,N", 
"8,12,W", "7,12,W", "6,12,W", "5,12,W", "4,12,W", "3,12,W", "2,12,W", "2,12,N", "2,13,N", 
"2,14,N", "2,14,E", "3,14,E", "4,14,E", "5,14,E", "6,14,E", "7,14,E", "8,14,E", "9,14,E", 
"10,14,E", "11,14,E", "12,14,E", "13,14,E", "14,14,E"], 6036]

We can then plot to compare this to the unweighted shortest path:

dsp := ListTools:-MakeUnique( map(s->s[1..-3], spd[1]) );
StyleVertex(G, dsp[2..-2], color="DarkBlue");
StyleEdge(G, [seq({dsp[i],dsp[i+1]}, i=1..nops(dsp)-1)], color="DarkBlue");

DrawGraph(G, stylesheet=[vertexshape="square", vertexpadding=10,
             vertexborder=false, vertexcolor="Black"],  showlabels=false,
          size=[800,800]);

And you can see it's a path that requires more steps, but definitely uses fewer turns if we start facing east/right (6 vs. 9):

I hope this has given you a little bit of a flavor of how to use GraphTheory commands to solve path finding problems.  Like with the second part here, usually the biggest challenge is figuring out how to encode and construct a graph that represents your problem.  Then the actual commands to solve it, are easy. You can see all the code, and a couple steps I left out from above in this worksheet: Mazeblog.mw

And just for fun, here's a Maple workbook that imports a maze from an image and solves it: MazeFromImage.maple

About once a year Advent of Code give you a problem that is a gift if you are using a computer algebra system. This year that day was Day 13. The day 13 problem was one of crazy claw machines.  Each machine has two buttons that can be pressed to move the claw a given number of X and Y positions and a prize at a given position. A buttons cost 3 tokens to press, and B buttons cost 1 token to press, and we are asked to find the minimum number of tokens needed, IF it is possible to reach the prize. The data presented like so:

Button A: X+94, Y+34
Button B: X+22, Y+67
Prize: X=8400, Y=5400

Button A: X+26, Y+66
Button B: X+67, Y+21
Prize: X=12748, Y=12176

...

Some times the input is harder to parse than the problem is so solve, and this might be one of those cases.  I tend to reach for StringTools to take the input apart, but today the right tool is the old school C-style sscanf (after using StringSplit to split at every double linebreak).

machinesL := StringTools:-StringSplit(Trim(input), "\n\n");
machines := map(m->sscanf(m,"Button A: X+%d, Y+%d\nButton B: X+%d, Y+%d\nPrize: X=%d, Y=%d"), machinesL):

Now we have a list of claw machine parameters in the form [A_x, A_y, B_x, B_y, P_x, P_y] and we need to turn those into equations that we can solve. We want the number of A presses a, and B presses b to get the claw to the P_x, P_y position of the claw, it is simple to just write them down:

for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
end do;

Now because of the discrete nature of this problem, we need our variables a and b to be non-negative integers.  When solving this, I first reached for isolve like this:

tokens := 0;
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   sol := isolve(eqn);
   if sol <> NULL then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;

Now, sometimes Advent of Code inputs contain a lot of hidden structure.  I wrote the code above, it worked on the sample input, so I tried it immediately on my real input (about 300 claw machines like the above) and IT WORKED.  But, you might notice that this code does not deal with a couple cases that could have appeared.  In particular, it doesn't check that the solutions are positive.  It also doesn't handle cases where there is more than one possible solution.  The former is easy to check

if sol <> NULL and eval(a,sol) >= 0 and eval(b,sol) >= 0 then

Unfortunately isolve does not handle inequalities, but you could try with solve, but it doesn't save us any checking, because we'd still have to check if the solutions are integers, so we might as well have just solved the equation and then checked if it were a nonnegative integer.

tokens := 0;
for m in machines do
   eqn := {m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6], a>=0, b>=0};
   sol := solve(eqn);
   if type(eval(a,sol), integer) and type(eval(b,sol),integer) then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;
ans1 = tokens;

In the multiple solution case we get something like {a = 3 - 2*b, 0 <= b, b <= 3/2} which has some great information in it but might be hard to handle programmatically, so let's see what isolve does with those cases to see if it's easier to deal with

> eqn := { 17*a + 84*b = 7870, 34*a + 168*b = 15740 }:
> constr := { a >= 0, b>= 0 }:
> sol := isolve(eqn);
               sol := {a = 458 - 84 _Z1, b = 1 + 17 _Z1}

> constr := eval({ a >= 0, b>= 0 }, sol):
            const := {0 <= 1 + 17 _Z1, 0 <= 458 - 84 _Z1}

> obj := eval(3*a+b, sol);
                         obj := 1375 - 235 _Z1

You can see it's easy to tell if these show up in your input, since your "token" total will have the _Zn variables in it.  Now, since everything is simple and linear here, it seems like you could use solve to find the rational value of _Z1 that makes obj=0 and then take the closest integer but it's not so simple, we actually have to deal with the contraints that a and b be positive too. So, it really just makes sense to bring out the big hammer of Optimization:-Minimize which allows us to directly optimize over just the integers.  So a full solution looks like this:

tokens := 0:
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   sol := isolve(eqn);
   if sol = NULL then
      next;
   end if;
   constr := eval({ a >= 0, b>= 0 }, sol);
   obj := eval(3*a+b, sol);
   if not type(obj, constant) then
      tokens := tokens + Optimization:-Minimize( obj, constr, assume=integer )[1];
   elif andmap(evalb, constr) then
      tokens := tokens + obj; 
   end if;
end do;

But since we're bringing out the big hammer, why not just use Optimization in the first place.  The main reason is that Minimize doesn't simply return NULL when it doesn't work, instead it throws an exception, so we need to find all the exceptions that can occur and handle then with a try-catch, thus:

tokens := 0;
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   try 
       sol := Optimization:-Minimize(3*a+b, eqn, 'assume'='nonnegint')[1];
       tokens := tokens + sol;
   catch "no feasible":
   end try; 
end do;
tokens;

(you can in fact omit the string in the catch: statement, but I can tell you from long experience that that is an excellent way to make your code much much harder to debug)

Alright, so how did people not using Maple solve this problem?  The easiest way to solve it, and the one used by all the cheaters scraping the website and using LLM-based code generators that auto-submit solutions to get into the Top 100, was to just check all possible a, b values in 0..100 and take the values than minimize 3*a+b when reaching the prize coordinates.  That's only feasible because the problem states the 100 is an upperbound for a and b, but it's also very fast (about 1/10 second in Maple):

tokens := 0:
for m in machines do;
sol := infinity;
for i from 0 to 100 do for j from 0 to 100 do
    if i*m[1]+j*m[3]=m[5] and i*m[2]+j*m[4]=m[6] and 3*i+j < sol then
        sol := 3*i+j;
    end if;
end do; end do;
tokens := tokens + ifelse(sol=infinity,0,sol);
end do;

It does not scale at all to part 2 (which modified everything to be bigger by about 10 trillion), and it seems that foiled all the LLM solvers. So, what solutions scaled in languages without integer equation solvers?  Well, the easiest solution is just to solve the general equation using paper and pencil

And you can just hard code that formula in, check that it gives integer values and compute the tokens. As long as you get unique solutions, that looks something like this

solveit := proc(m)
local asol := m[4]*m[5] - m[3]*m[6];
local bsol := m[1]*m[6] - m[2]*m[5];
local deno := m[1]*m[4] - m[2]*m[3];
if deno = 0 then return -2^63; end if; # multiple solution case - not handled
if asol mod deno = 0 and bsol mod deno = 0
   and (   ( deno>=0 and asol>=0 and bsol>=0 ) 
        or ( deno<=0 and asol<=0 and bsol<=0 ) )
then

    return 3*asol/deno+bsol/deno;
else
    return 0;
end if;
end proc:

Which if you have this in Maple, you can impress your friends by auto generating solutions in other languages. Here, for your FORTRAN friends

> CodeGeneration:-Fortran(solveit);

Warning, the following variable name replacements were made: solveit -> cg
       integer function cg (m)
        doubleprecision m(*)
        integer asol
        integer bsol
        integer deno
        asol = int(-m(3) * m(6) + m(4) * m(5))
        bsol = int(m(1) * m(6) - m(2) * m(5))
        deno = int(m(1) * m(4) - m(2) * m(3))
        if (deno .eq. 0) then
          cg = -9223372036854775808
          return
        end if
        if (mod(asol, deno) .eq. 0 .and. mod(bsol, deno) .eq. 0 .and. (0
     # .le. deno .and. 0 .le. asol .and. 0 .le. bsol .or. deno .le. 0 .a
     #nd. asol .le. 0 .and. bsol .le. 0)) then
          cg = 3 * asol / deno + bsol / deno
          return
        else
          cg = 0
          return
        end if
      end

Another way that you might solve this without solve is to use a linear algebra library to solve the linear system.  It works even if you only have a numeric solver, but you have to be careful about checking for integers:

tokens := 0:
for m in machines do
    sol := LinearAlgebra:-LinearSolve(
               Matrix(1..2,1..2,[m[[1,3]],m[[2,4]]], datatype=float), 
               Vector(m[5..6], datatype=float));
    if abs(sol[1]-round(sol[1])) < 10^(-8) and abs(sol[2]-round(sol[2])) < 10^(-8)
       and sol[1] >= 0 and sol[2] >= 0
    then
       tokens := tokens + 3*sol[1]+sol[2];
    end if;
end do;

Finally, a lot of people solved this sort of thing with the Z3 Theorem prover from Microsoft research which is also way more than you need, but it mostly just uses SMTLIB, which we also have in a library for in Maple, and it can just be used in place of solve

tokens := 0;
for m in machines do
   eqn := {m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6], a>=0, b>=0};
   sol := SMTLIB:-Satisfy(eqn) assuming a::nonnegint, b::nonnegint;
   if sol <> NULL and type(eval(a,sol), integer) and type(eval(b,sol),integer) then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;

Notice that Satisfy handled the multiple solution case just by choosing one of the many solutions. It is possible to get SMTLib to optimize but it is slightly more involved, and this post is already too long. This time, I've put all this work in worksheet: Day13-Primes.mw

The Advent of Code 2024 Day 11 problem featured labeled rocks which duplicate according to simple rules:

• Rocks labeled 0 turn into 1's
• Rocks with a even number of digits in the label split into two rocks with the first and last half of the digits
• Other rocks have their labels multiplied by 2024

Your puzzle task is to count the number of rocks after 25, then 75 iterations of these rules.

For 25 iterations, a very simple recursive count suffices

count := proc(n, s)
local f, l;
    if n = 0 then return 1;
    elif s=0 then return count(n-1,1);
    end if;
    l := length(s);
    if l mod 2 = 0 then
         f := floor( s/10^(l/2) );
         return count(n-1,f) + count(n-1,s-f*10^(l/2));
    else return count(n-1, 2024*s);
    end if;
end proc:

For the puzzle input, a collection of 8 numbers (the stone labels) the largest with 7 digits, you can call count(25,s) on each s in the input and add them together in about 1 second, but trying that for 75 "hangs" and clearly will not work.

Generally for something like this, the first trick in my puzzle solving toolbox is memoization, so that the results of recursive calls are stored so they don't have to be recomputed.  In Maple, there are two main ways to add memoization to recursive procedures: option cache and option remember.  The remember option is older and uses a simple table that gets populated with (input,return) pairs whenever the procedure returns. The cache option is newer and uses size limited LRU tables that are more memory friendly.  Typically in the Maple math library we use cache tables since they won't fill up your working memory with cached results, and will remove old cached elements in a smart way (the double option option remember, system also limits table growth by allowing entries to be garbage collected, but it is not as flexible as a cache).

So this preable is to say, that the first memoization option I tend to reach for is option cache but best practice for Maple math library development is not always the thing that solves a code puzzle.  In this case, adding option cache to our procedure count also "hangs".  If you are smart, you will immediately try again with option remember, and you'll find that works easily:

count := proc(n, s)
option remember;
    if n = 0 then return 1;
    elif s=0 then return count(n-1,1);
    end if;
    local l := length(s);
    if l mod 2 = 0 then
         local f := floor( s/10^(l/2) );
         return count(n-1,f) + count(n-1,s-f*10^(l/2));
    else return count(n-1, 2024*s);
    end if;
end proc:
ans1 := CodeTools:-Usage( add( map2(count, 25, input) ) ); # 17ms
ans2 := CodeTools:-Usage( add( map2(count, 75, input) ) ); # 825ms

So why does remember work, when cache fails so spectacularly? Well, you can create a function with a cache, call it once create it and the examine the Cache object with op:

> cacheFunc := proc() option cache; args; end proc:
> cacheFunc(1): op( 4, eval(cacheFunc) );

                       Cache(512, 'temporary' = [1 = 1])

You can see here that the default cache is 512 entries which is a pretty good default for general use, but way way too small for this ridiculous puzzle problem.  Fortunately, the cache option allows you to specify a size to create a much bigger cache for a highly recursive procedure like this one.  I tested caches of various sizes (cache always rounds up the next power of 2) and compared to option remember and remember, system.

1. option remember  - 707ms
2. option remember, system - 4.22s
3. option cache  - too slow, killed after hours
4. option cache(1024) - 30.21s
5. option cache(16384) - 5.60s
6. option cache(131072)- 726ms
7. option cache(1048576)- 806ms

So, for this problem, it looks like a 2^17=~100,000 entry cache (#6) gives about the same performance as a remember table, so it's not surpring a larger cache does no better. 

In fact, since results are not removed from an option remember table, you can actually check exactly how many results were cached and see it is just less than 2^17.

> add( map2(count, 75, input) ): numelems( op(4, eval(count) ) );
                                    120076

And the take away is, for this sort of puzzle solving you should try option remember even if option cache is usually a better choice when implementing a more serious mathematical algorithmn that wants to occasionally needs to reuse results.

Every December since 2015, software engineer and puzzle enthusiast Eric Wastl has created an Advent Calender of code challenges called Advent of Code. These puzzles can be solved in any programming language, since you are only required to submit a string or integer solution, and tend to increase in difficulty over the course of the 25 days. Each puzzle has two parts each using the same input (which is different for each user), the second part being revealed only after the first part is solved.

I started participating in this puzzle challenge in 2021 and I like to solve them using Maple right when they are unlocked at 12am EST each day. I post my my solutions to my github as Maple worksheets, and (usually) later as cleaned up Maple text files: https://github.com/johnpmay/AdventOfCode2024

I typically work in a Maple worksheet since I find the worksheet interface ideal to quickly play with an evolving piece of code to solve a puzzle (it is kind of a race if you want it to be, Advent of Code tells you your place every day and gives you "points" if you are one of the first 100). Also I find the nature of the Maple language and its variety of libraries pretty ideal for approaching these problems.

Today, I am writing about this because 2025 Day 5 was a cool puzzle that really asked to be analyzed further in Maple. (Note, if you are participating in Advent of Code and haven't solved Day 5, there are heavy spoilers ahead).

Once you decypher the silly story about the elves, the problem comes down to a set of pairwise ordering rules, and a collection of lists (called updates) that need to be checked if they are sorted according to the rule. This seems straightforward, but as always you have a (big!) input file to handle.  In this case the input looks like this (the actual input not included here, but it's much much longer)

75|13
53|13
[...elided...]


75,47,61,53,29
97,61,53,29,13
[...elided...]

For this step, StringTools is your best friend (I typically read in the whole input as a string using FileTools:-Text:-ReadFile).  I first use StringSplit on the double line break "\n\n" to seperate the rules from the updates.  Then I map an ordinary Split over those to make everything into nested lists of integer.

(rules, updates) := StringSplit(Trim(input), "\n\n")[]:
rules := map(s->map(s2i,Split(s,"|")), Split(rules,"\n")): 

        rules := [[75, 13], [53, 13], ...]

updates := map(s->(map(s2i,Split(s,","))), Split(updates, "\n")):

        updates := [[75, 47, 61, 53, 29], [97, 61, 53, 29, 13], ...]

(the function s2i is a "string to integer" helper function since it's considered bad practice to just call parse on input handed to you by a stranger on the internet.  In older version of maple I used: s2i := s->sscanf(s,"%d")[1] in newever versions you can use the way less cryptic s2i := s->convert(s,'integer') )

Okay, now that we have two nice lists, it's easy to check which of the updates is sorted according to the rules with a simple nested loop.

goodups := DEQueue():

for u in updates do
    for r in rules do
        if member(r[1],u,'i') and member(r[2],u,'j') then
           if j < i then next 2; end if;
        end if;
     end do;
     push_back(goodups, u);
end do:

For each update, check that if a rule applies to it then that rule is obeyed. If the rule r is violated, advance the outer loop with the relatively new and useful next 2; syntax (it also works with break #) and if not check the next rule. If no rules were violated then push u into our stack of Good Updates.  This is a perfectly fine solution to part 1 of the puzzle, but if you thing about it a little why not just use sort with a custom comparision function built from the rules?

V := ListTools:-MakeUnique(map2(op,1,rules)): # all the numbers that can appear
N := [seq(map2(op,2,select(n->m=n[1], rules)), m in V)]: # all their neighbors in the graph induced by rules
T := table([seq(V[i]={N[i][]}, i=1..nops(V))]): # table of neighbors
comp := (x,y)->evalb(y in T[x]): # "less than" according the rules

Now with that comp function we can just do

goodups := select(u->sort(u, comp)=u, updates);

to find all the correctly sorted updates, this will seem like a better idea when you get the [SPOILER] part 2 of the problem which is to sort all the incorrectly sorted updates.

NOW WAIT A MINUTE.  If you look at the fine help page for ?sort you'll see it very clearly says about the comp function F (emphasis added):

custom boolean procedure

 When F does not match one of the symbols listed above, it must be a Boolean-
 valued function of two arguments.  Specifically, F(a,b) returns false if and  
 only if b must precede a in the sorted output.  That is F(a,b) is a non-
 strict less than comparison function.  In addition, F(a,b) must be defined  
 for all pairs a for a and b in the input structure and F(a,b) must be  
 transitive, that is, if F(a,b) = true and F(b,c) = true then F(a,c) = true.

Dear reader, I guess we should check if our rules from the puzzle create a relation that is transitive.  The good news is that the sample input given in the problem description is indeed transitive. The easiest way I could think to check that was to use GraphTheory and see if there were any cycles in the directed graph given by the rules. So, let's make a graph with an edge for every comparison rule a<b, b<c, etc.  If there is a cycle in the graph, that will be a case where a < b < c < a in our rules and therefore non-transitivity in the comparison.  When we do that, we see it's cycle-free.

with(GraphTheory):
G := Graph({rules[]});
FindCycle(G); # returns []

That graph is really nice, it looks like this:

If you count the in-degree of each vertex, you can see it actually induced a total order on the numbers going linearly from smallest to largest numbers.

So, is that also true of the actual input you are given? Of course not. When you create the graph of your actual puzzle input you find there is a 3-cycle.

with(GraphTheory):
G := Graph({rules[]});
FindCycle(G); # returns a 3-cycle
seq(FindCycle(G,v),v in Vertices(V)); # returns a 3-cycle for every v!

In fact, there is a 3-cycle from every single vertex.  But at least the graph is pleasing to contemplate. Let it calm you:

So, does this mean sort with out comp function doesn't work? Is all lost? In fact, I solved both parts of the puzzle using a custom comparision for sort before stopping to consider that maybe the comparison wasn't transitive. And... it just works! So, WHY does it work? (the comparison is very very non-transitive). The answer is of course that the lists you need to sort don't contain all of the numbers covered by the rules.  In fact, your clever/cruel puzzlemaster has in fact constructed each of those (200 or so) lists so that it misses every one of the 3-cycles in the above graph. You can carefully check that by restricting the relation graph to the subset of numbers in each list and check for cycles:

G := Graph({rules[]});
H := InducedSubgraph(G, updates[i]); # try for i=1..nops(updates)
FindCycles(H); # returns [] for every element of update

In fact, each of those induced subgraphs looks like this

and you can maybe observe that this induced a total order on the numbers in the problem. Or if you are not careful, and you didn't bother to check, you would get the right answer anyway. Because it's only Day 5, and Eric Wastl hasn't decided to make thing really challenging yet.

I hope you've enjoyed my deep dive into this problem. And I hope it's inspired you to go try out Advent of Code. It's a fun challenge. Let me know in a comment if you're participating! Maybe we can start a leaderboard for MaplePrimes folks working on the puzzles.

I was working in my living room.  My computer was upstairs, but I had my phone and tablet.  I'm working on The Book ("Perturbation methods using backward error", with Nic Fillion, which will be published by SIAM next year some time).

I've discovered something quite cool, historically, about the WKB method and George Green's original invention of the idea (that bears other people's names, or, well, initials, anyway).  (As usual.)  Green had written down a PDE modelling waves in a long narrow canal of slowly varying breadth 2*beta(x) and slowly varying depth 2*gamma(x).  Turns out his "approximate" solution is actually an exact solution to an equation of a very similar kind, with an extra term E(x)*phi(x,t).  The extra term depends in a funny way on beta(x) and gamma(x), and only on those.  So a natural kind of question is, "is there a canal shape for which Green's solution is the exact solution with E(x)==0?"  Can we find beta(x) and gamma(x) for which this works?

Yes.  Lots of cases.  In particular, if the breadth beta(x) is constant, you can write down a differential equation for gamma(x).  I wrote it in my notebook using y and not gamma.  I wrote it pretty neatly.  Then I fired up the Maple Calculator on my little tablet, opened the camera, and pow!  Solved.

I wrote the solution down underneath the equation.  It checks out, too.  See the attached image.

Now, after the fact, I figured out how to solve it myself (using Ricatti's trick: put y' = v, then y'' = v*dv/dy, and the resulting first order equation is separable).  But that whole "take a picture of the equation and write down the solution" thing is pretty impressive.

So: kudos to the designers and implementers of the Maple Calculator.  Three cheers!

Hello Maple enthusiasts,

I am excited to share a sample worksheet on Ordinary Differential Equations (ODEs), created as part of my ongoing project—a book I am writing for undergraduate students. This book is designed to teach ODEs using Maple, offering an interactive and intuitive approach to solving differential equations.

As far as I know, there aren’t many books available in the Greek language that combine ODEs with Maple. In fact, I believe there’s only one other such resource, which highlights the lack of materials in this niche. My goal is to fill this gap by providing students and educators with a resource that is both practical and accessible, leveraging Maple's powerful capabilities to deepen understanding and simplify complex concepts.

The worksheet I’m sharing includes:

• Step-by-step solutions to ODEs using Maple.
• Graphical representations to visualize solutions, which I believe are invaluable for fostering comprehension.

I hope this preview sparks your interest and provides insight into the teaching style and structure of the upcoming book. I would love to hear your thoughts, feedback, or suggestions for topics you think should be included.

Download separable_diff_equations.mw

After more than 25 years of leading research in areas such as differential equations, special functions, and computational physics, Edgardo's role with Maplesoft will shift at the end of 2024 as he returns to academic research. At Maplesoft, he will transition into the position of Research Fellow Emeritus. In this role, Edgardo will remain engaged with many of his cherished projects, though he will not have as much time to maintain the intense level of activity that characterized his work for so many years.

Many of you know Edgardo personally or have interacted with him here or on the Maple Beta Forum. I hope you'll join me in wishing him the very best as he begins this new chapter.

Over the last few days, I’ve been creating worksheets on oscillators to support my class’s understanding of these fundamental physics concepts. I wanted to share one of these worksheets that I found particularly useful for illustrating energy exchange and motion dynamics.

Download Simple_Pendulum.mw

Maple Transactions has just published the Autumn 2024 issue at mapletransactions.org

From the header:

This Autumn Issue contains a "Puzzles" section, with some recherché questions, which we hope you will find to be fun to think about.  The Borwein integral (not the Borwein integral of XKCD fame, another one) set out in that section is, so far as we know, open: we "know" the value of the integral because how could the identity be true for thousands of digits but yet not be really true? Even if there is no proof.  But, Jon and Peter Borwein had this wonderful paper on Strange Series and High Precision Fraud showing examples of just that kind of trickery.  So, we don't know.  Maybe you will be the one to prove it! (Or prove it false.)

We also have some historical papers (one by a student, discussing the work of his great grandfather), and another paper describing what I think is a fun use of Maple not only to compute integrals (and to compute them very rapidly) but which actually required us to make an improvement to a well-known tool in asymptotic evaluation of integrals, namely Watson's Lemma, just to explain why Maple is so successful here.

Finally, we have an important paper on rational interpolation, which tells you how to deal well with interpolation points that are not so well distributed.

Enjoy the issue, and keep your contributions coming.

With the 2024 Maple Conference coming up this week, I imagine one of two of you have noticed something missing. We chose not to have a Conference Art Gallery this year, because we have been working to launch new section of MaplePrimes:  The MaplePrimes Art Gallery. This new section of MaplePrimes is designed for showing off your Maple related images, in a gallery format that puts the images up front, like Instagram but for Math.

To get the ball rolling on the gallery, I have populated it with many of the works from past years' Maple Art Galleries, so you can now browse them all in one place, and maybe give "Thumbs Ups" to art pieces that you felt should have won the coveted "People's Choice Award".

Once you are inspired by past entries, we welcome you to submit your new artwork!  Just like the conference galleries, we are looking for all sorts of mathematical art ranging from computer graphics and animations, to photos of needlework, geometrical sculptures, or almost anything!  Art Gallery posts are very similar to regular MaplePrimes posts except that you are asked to supply an Image File and a Caption that will displayed on the main Gallery Page, the post itself can include a longer description, Maple commands, additional images, and whatever else you like.  See for example this Art Gallery post describing Paul DeMarco's sculpture from the 2021 Maple Conference Gallery - it includes an embedded worksheet as well as additional images.

I can't wait to see what new works of art our MaplePrimes community creates!

I recently prepared a worksheet to teach vector fundamentals in one of my classes, and I wanted to share it with you all. It's nothing special, but I found Maple really helpful in demonstrating the concepts visually. Below is a breakdown of what the worksheet covers, with some Maple code examples included.

Feel free to take a look and use it if you find it useful! Any feedback or suggestions on how to improve it would be appreciated.

Download vectors.mw

Just finished an exciting lecture for undergraduate students on Euler methods for solving ordinary differential equations (ODEs)! 📝 We explored the Euler method and the Improved Euler method (a.k.a. Heun's method) and discussed how these fundamental techniques work for approximating solutions to ODEs.

To make things practical and hands-on, I demonstrated how to implement both methods using Maple—a fantastic tool for experimenting with numerical methods. Seeing these concepts come to life with visualizations helps understand the pros and cons of each method.

The Euler method is the simplest numerical approach, where we approximate the solution by stepping along the slope given by the ODE. But there's a catch: using a large step size can lead to large errors that accumulate quickly, causing significant deviation from the real solution.

Plot Results:

• With a larger step size (h=0.5h = 0.5h=0.5), the solution tends to drift away significantly from the actual curve.
• Reducing the step size improves accuracy, but it also means more steps and more computation.

Next, we discussed the Improved Euler method, which is like Euler's smarter sibling. Instead of blindly following the initial slope, it:

1. Estimates the slope at the start.
2. Uses that slope to predict an intermediate value.
3. Then, it takes another slope at the intermediate point.
4. Averages these two slopes for a better approximation.

This technique makes a big difference in accuracy and stability, especially with larger step sizes.

Plot Results:

• Using the Improved Euler method, we found that even with a larger step size, the solution is smoother and closer to the true path.
• The average slope helps mitigate the inaccuracies that arise from only using the beginning point's derivative, effectively reducing the local error.

• The standard Euler method can produce solutions that oscillate or diverge, especially for larger step sizes.
• The Improved Euler method follows the actual trend of the solution more faithfully, even with fewer steps. This makes it a more efficient choice for balancing computational effort and accuracy.

• The Euler method is great for its simplicity but often requires very small step sizes to be accurate, leading to more computational effort.
• The Improved Euler method—which we implemented and visualized on Maple—proved to be more reliable and accurate, especially for larger step sizes.

It was amazing to see the students engage with these foundational methods, and implementing them on Maple brought a deeper understanding of numerical analysis and the challenges of solving differential equations.

Download Diff_equations.mw