Compute Derivatives: Yes or No

Numerical Evaluation

Two Hankel Transform Definitions

More integral transform results

Mellin and Hankel transform solutions for Partial Differential Equations with boundary conditions

Compute derivatives: Yes or No.

For historical reasons, previous implementations of these integral transform commands did not follow a standard paradigm of computer algebra: "Given a function , the input  should return the derivative of ". The implementation instead worked in the opposite direction: if you were to input the result of the derivative, you would receive the derivative representation. For example, to the input  you would receive . This is particularly useful for the purpose of using integral transforms to solve differential equations but it is counter-intuitive and misleading; Maple knows the differentiation rule of these functions, but that rule was not evident anywhere. It was not clear how to compute the derivative (unless you knew the result in advance).

To solve this issue, a new command, setup, has been added to the package, so that you can set "whether or not" to compute derivatives, and the default has been changed to computederivatives = true while the old behavior is obtained only if you input . For example, after having installed the Physics Updates,

the current settings can be queried via

and so differentiating returns the derivative computed

while changing this setting to work as in previous releases you have this computation reversed: you input the output (1.3) and you get the corresponding input

Reset the value of computederivatives

In summary: by default, derivatives of integral transforms are now computed; if you need to work with these derivatives as in  previous releases, you can input . This setting can be changed any time you want within one and the same Maple session, and changing it does not have any impact on the performance of intsolve, dsolve and pdsolve to solve differential equations using integral transforms.

Numerical Evaluation

In previous releases, integral transforms had no numerical evaluation implemented. This is in the process of changing. So, for example, to numerically evaluate the inverse laplace transform ( invlaplace  command), three different algorithms have been implemented: Gaver-Stehfest, Talbot and Euler, following the presentation by Abate and Whitt, "Unified Framework for Numerically Inverting Laplace Transforms", INFORMS Journal on Computing 18(4), pp. 408–421, 2006.

For example, consider the exact solution to this partial differential equation subject to initial and boundary conditions

Note that these two conditions are not entirely compatible: the solution returned cannot be valid for  and  simultaneously. However, a solution discarding that point does exist and is given by

Verifying the solution, one condition remains to be tested

Since we now have numerical evaluation rules, we can test that what looks different from 0 in the above is actually 0.

Add a small number to the initial value of t to skip the point

The default method used is the method of Euler sums and the numerical evaluation is performed as usual using the evalf command. For example, consider

The Laplace transform of F is given by

and the inverse Laplace transform of LT in inert form is

At  we have

This result is consistent with the one we get if we first compute the exact form of the inverse Laplace transform at t = 1:

In addition to the standard use of evalf to numerically evaluate inverse Laplace transforms, one can invoke each of the three different methods implemented using the MathematicalFunctions:-Evalf  command

Regarding the method we use by default: from a numerical experiment with varied problems we have concluded that our implementation of the Euler (sums) method is faster and more accurate than the other two.

Two Hankel transform definitions

In previous Maple releases, the definition of the Hankel transform was given by

where  is the  function. This definition, sometimes called alternative definition of the Hankel transform, has the inconvenience of the square root  in the integrand, complicating the form of the hankel transform for the Laplacian in cylindrical coordinates. On the other hand, the definition more frequently used in the literature is

With it, the Hankel transform of  is given by the simple ODE form . Not just this transform but several other ones acquire a simpler form with the definition that does not have a square root in the integrand.

So the idea is to align Maple with this simpler definition, while keeping the previous definition as an alternative. Hence, by default, when you load the inttrans package, the new definition in use for the Hankel transform is

You can change this default so that Maple works with the alternative definition as in previous releases.  For that purpose, use the new inttrans:-setup command (which you can also use to query about the definition in use at any moment):

This change in definition is automatically taken into account by other parts of the Maple library using the Hankel transform. For example, the differentiation rule with the new definition is

This differentiation rule resembles (is connected to) the differentiation rule for BesselJ, and this is another advantage of the new definition.

Furthermore, several transforms have acquired a simpler form, as for example:

Let's compare: make the definition be as in previous releases.

The differentiation rule with the previous (alternative) definition was not as simple:

And the transform (3.5) was also not so simple:

Reset to the new default value of the definition.

More integral transform results

The revision of the integral transforms includes also filling gaps: previous transforms that were not being computed are now computed. Still with the Hankel transform, consider the operators

Being able to transform these operators into algebraic expressions or differential equations of lower order is key for solving PDE problems with Boundary Conditions.

Consider, for instance, this ODE

Its Hankel transform is a simple algebraic expression

An example with formula_plus

In the case of hankel , not just differential operators but also several new transforms are now computable

Mellin and Hankel transform solutions for Partial Differential Equations with Boundary Conditions

In previous Maple releases, the Fourier and Laplace transforms were used to compute exact solutions to PDE problems with boundary conditions. Now, Mellin and Hankel transforms are also used for that same purpose.

Example:

As usual, you can let pdsolve choose the solving method, or indicate the method yourself:

It is sometimes preferable to see these solutions in terms of more familiar integrals. For that purpose, use

An example where the hankel transform is computable to the end:

restart;

SerpentinePaths:=proc(S::list, P::list:=[1,1])
local OneStep, A, m, F, B, T, a;

OneStep:=proc(A::listlist)
local s, L, B, T, k, l;

s:=max[index](A);
L:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]];
T:=table(); k:=0;
for l in L do
if l[1]>=1 and l[1]<=S[1] and l[2]>=1 and l[2]<=S[2] and A[op(l)]=0 then k:=k+1; B:=subsop(l=a+1,A);
T[k]:=B fi;
od;
convert(T, list);
end proc;
A:=convert(Matrix(S[], {(P[])=1}), listlist);
m:=S[1]*S[2]-1;
B:=[$ 1..m];
F:=LM->ListTools:-FlattenOnce(map(OneStep, `if`(nops(LM)<=30000,LM,LM[-30000..-1])));
T:=[A];
for a in B do
T:=F(T);
od;
map(convert, T, Matrix);

end proc:
 

NumbrixPuzzle:=proc(A::Matrix)
local A1, L, N, S, MS, OneStepLeft, OneStepRight, F1, F2, m, L1, p, q, a, b, T, k, s1, s, H, n, L2, i, j, i1, j1, R;
uses ListTools;
S:=upperbound(A); N:=nops(op(A)[3]); MS:=`*`(S);
A1:=convert(A, listlist);
for i from 1 to S[1] do
for j from 1 to S[2] do
for i1 from i to S[1] do
for j1 from 1 to S[2] do
if A1[i,j]<>0 and A1[i1,j1]<>0 and abs(A1[i,j]-A1[i1,j1])<abs(i-i1)+abs(j-j1) then return `no solutions` fi;
od; od; od; od;
L:=sort(select(e->e<>0, Flatten(A1)));
L1:=[`if`(L[1]>1,seq(L[1]-k, k=0..L[1]-2),NULL)];
L2:=[seq(seq(`if`(L[i+1]-L[i]>1,L[i]+k,NULL),k=0..L[i+1]-L[i]-2), i=1..nops(L)-1), `if`(L[-1]<MS,seq(L[-1]+k,k=0..MS-L[-1]-1),NULL)];
OneStepLeft:=proc(A1::listlist)
local s, M, m, k, T;
uses ListTools;
s:=Search(a, Matrix(A1));   
M:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]];
T:=table(); k:=0;
for m in M do
if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 then k:=k+1; T[k]:=subsop(m=a-1,A1);
fi;
od;
convert(T, list);
end proc;
OneStepRight:=proc(A1::listlist)
local s, M, m, k, T, s1;
uses ListTools;
s:=Search(a, Matrix(A1));  s1:=Search(a+2, Matrix(A1));  
M:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]];
T:=table(); k:=0;
for m in M do
if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 and `if`(a+2 in L, `if`(is(abs(s1[1]-m[1])+abs(s1[2]-m[2])>1),false,true),true) then k:=k+1; T[k]:=subsop(m=a+1,A1);
fi;
od;
convert(T, list);   
end proc;
F1:=LM->ListTools:-FlattenOnce(map(OneStepLeft, LM));
F2:=LM->ListTools:-FlattenOnce(map(OneStepRight, LM));
T:=[A1];
for a in L1 do
T:=F1(T);
od;
for a in L2 do
T:=F2(T);
od;
R:=map(t->convert(t,Matrix), T);
if nops(R)=0 then return `no solutions` else R fi;
end proc:

SerpentinePaths([3,3]);  # All the serpentine paths for the matrix  3x3, starting with [1,1]-position
SerpentinePaths([3,3],[1,2]);  # No solutions if the start with [1,2]-position
SerpentinePaths([4,4]):  # All the serpentine paths for the matrix  4x4, starting with [1,1]-position
nops(%);
nops(SerpentinePaths([4,4],[1,2]));  # The number of all the serpentine paths for the matrix  4x4, starting with [1,2]-position
nops(SerpentinePaths([4,4],[2,2]));  # The number of all the serpentine paths for the matrix  4x4, starting with [2,2]-position

k:=0:  n:=6:
for i from 1 to n do
for j from i to n do
k:=k+1; S[k]:=SerpentinePaths([n,n],[i,j])[];
od: od:
L1:={seq(S[i][], i=1..k)}:
L2:=map(A->A^%T, L1):
L:=L1 union L2:
nops(L);

T:='T':
N:=20:
M:=[seq(L[i], i=combinat:-randcomb(nops(L),N))]:
for i from 1 to N do
for k from floor(n^2/4) do
T[i]:=Matrix(n,{seq(op(M[i])[3][j], j=combinat:-randcomb(n^2,k))});
if nops(NumbrixPuzzle(T[i]))=1 then break; fi;
od:  od:
T:=convert(T,list):
T1:=[seq([seq(T[i+j],i=1..5)],j=[0,5,10,15])]:
DocumentTools:-Tabulate(Matrix(4,5, (i,j)->T1[i,j]), fillcolor = "LightYellow", width=95):

DocumentTools:-Tabulate(Matrix(4,5, (i,j)->NumbrixPuzzle(T1[i,j])[]), fillcolor = "LightYellow", width=95):

The random number generator needs to be seeded differently in each node. (The reason for this is easy to understand.)

The random variables generated by Statistics:-RandomVariable need to have different names in each node. This one is mind-boggling to me. Afterall, each node has its own kernel and so its own memory It's as if the names of random variables are stored in a disk file which all kernels access. And also the generator has been seeded differently in each node.

restart
:

Digits:= 15
:

#Specify the size of the computation:
(n1,n2,n3):= (100, 100, 1000):
# n1 = size of each random sample;
# n2 = number of samples in a batch;
# n3 = number of batches.

#
#Procedure to initialize needed globals on each node:
Init:= proc(n::posint)
local node:= Grid:-MyNode();
   #This is wrapped in parse so that it'll act globally. Otherwise, an environment
   #variable would be reset when this procedure ends.
   parse("Digits:= 15;", 'statement');

   randomize(randomize()+node); #Initialize independent RNG for this node.
   #If repeatability of results is desired, remove the inner randomize().

   (:-X,:-Y):= Array(1..n, 'datatype'= 'hfloat') $ 2;

   #Perhaps due to some oversight in the design of Statistics, it seems necessary that
   #r.v.s in different nodes **need different names** in order to be independent:
   N||node:= Statistics:-RandomVariable('Normal'(0,1));
   :-TRS:= (X::rtable)-> Statistics:-Sample(N||node, X);
   #To verify that different names are needed, change N||node to N in both lines.
   #Doing so, each node will generate identical samples!

   #Perform some computation. For the pedagogical purpose of this worksheet, all that
   #matters is that it's some numeric computation on some Arrays of random Samples.
   :-GG:= (X::Array, Y::Array)->
      evalhf(
         proc(X::Array, Y::Array, n::posint)
         local s, k, S:= 0, p:= 2*Pi;
            for k to n do
               s:= sin(p*X[k]);  
               S:= S + X[k]^2*cos(p*Y[k])/sqrt(2-sin(s)) + Y[k]^2*s
            od
         end proc
         (X, Y, n)
      )      
   ;
   #Perform a batch of the above computations, and somehow numerically consolidate the
   #results. Once again, pedagogically it doesn't matter how they're consolidated.  
   :-TRX1:= (n::posint)-> add(GG(TRS(X), TRS(Y)), 1..n);
   
   #It doesn't matter much what's returned. Returning `node` lets us verify that we're
   #actually running this on a grid.
   return node
end proc
:

#
#Time the initialization:
st:= time[real]():
   #Send Init to each node, but don't run it yet:
   Grid:-Set(Init)
   ;
   #Run Init on each node:
   Nodes:= Grid:-Run(Init, [n1], 'wait');
time__init_Grid:= time[real]() - st;

num_nodes:= numelems(Nodes);

#Time the actual execution:
st:= time[real]():
   R1:= [Grid:-Seq['tasksize'= iquo(n3, num_nodes)](TRX1(k), k= [n2 $ n3])]:
time__run_Grid:= time[real]() - st

#Just for comparison, run it sequentially:
st:= time[real]():
   Init(n1):
time__init_noGrid:= time[real]() - st;

st:= time[real]():
   R2:= [seq(TRX1(k), k= [n2 $ n3])]:
time__run_noGrid:= time[real]() - st;

plots:-display(
   Statistics:-Histogram~(
      <R1 | R2>, #side-by-side plots
      'title'=~ <<"With Grid\n"> | <"Without Grid\n">>,
      'gridlines'= false
   )
);

speedup_factor:= time__run_noGrid / time__run_Grid;

efficiency:= speedup_factor / num_nodes;

Background

Parade magazine, a filler in the local Sunday newspaper, contains a Numbrix puzzle, the object of which is to find a serpentine path of consecutive integers, 1 through 81, through a nine by nine grid.  The puzzle typically specifies the content of every other border cell.

The Maple Logic  package has a procedure, Satisfy , that can be used to solve this puzzle.  Satisfy is a SAT-solver; given a boolean expression it attempts to find a set of equations of the form , where  are the boolean variables in the given expression and  are boolean values (true or false) that satisfy the expression (cause it to evaluate true).

A general technique to solve this and other puzzles with Satisfy is to select boolean-values variables that encode the state of the puzzle (a trial solution, whether valid or not), generate a boolean-expression of these variables that is satisfied (true) if and only if the variables are given values that correspond to a solution, pass this expression to Satisfy, then translate the returned set of boolean values (if any) to the puzzle solution.

Setup

Assign a matrix that defines the grid and the initial position.  Use zeros to indicate the cells that need values. To make it easy to inspect the expressions, a small 2 x 3 matrix is used for this demo---a full size example is given at the end.

Extract the dimensions of the Matrix

Boolean Variables

Let the boolean variable x[i,j,k] mean that cell (i,j) has value k. For example, x[2,3,6] is true when cell (2,3) contains 6, otherwise it is false. There are  boolean variables.

Initial Position

The initial position can be expressed as the following and-clause.

Adjacent Cells

The requirement that an interior cell with value k is adjacent to the cell with value k+1 can be expressed as the implication
   

   x[i,j,k] &implies &or(x[i-1,j,k+1], x[i+1,j,k+1], x[i,j-1,k+1], x[i,j+1,k+1])

Extending this to handle all cells results in the following boolean expression.

All Values Used

The following expression is true when each integer , from 1 to , is assigned to one or more cells.

Single Value

The following expression is satisfied when each cell has no more than one value.

Solution

Combine the boolean expressions into a a single and-clause and pass it to Satisfy.

Select the equations whose right size is true

Extract the lhs of the true equations

Extract the result from the indices of the vars and assign to a new Matrix

Procedure

We can now combine the manual steps into a procedure that takes an initialized Matrix and fills in a solution.

Example

Create the initial position for a 9 x 9 Numbrix and solve it.

restart;

contplot:=proc(ee, rng1, rng2)
  local clabels, clegend, i, ncrvs, newP, otherdat, others, tcrvs, tempP;
  (clegend,others):=selectremove(type,[_rest],identical(:-legend)=anything);
  (clabels,others):= selectremove(type,others,identical(:-contourlabels)=anything);
  if nops(clegend)>0 then
    tempP:=:-plots:-contourplot(ee,rng1,rng2,others[],
                                ':-contourlabels'=rhs(clegend[-1]));
    tempP:=subsindets(tempP,'specfunc(:-_HOVERCONTENT)',
                      u->`if`(has(u,"null"),NULL,':-LEGEND'(op(u))));
    if nops(clabels)>0 then
      newP:=plots:-contourplot(ee,rng1,rng2,others[],
                              ':-contourlabels'=rhs(clabels[-1]));
      tcrvs:=select(type,[op(tempP)],'specfunc(CURVES)');
      (ncrvs,otherdat):=selectremove(type,[op(newP)],'specfunc(CURVES)');
      return ':-PLOT'(seq(':-CURVES'(op(ncrvs[i]),op(indets(tcrvs[i],'specfunc(:-LEGEND)'))),
                          i=1..nops(ncrvs)),
                      op(otherdat));
    else
      return tempP;
    end if;
  elif nops(clabels)>0 then
    return plots:-contourplot(ee,rng1,rng2,others[],
                              ':-contourlabels'=rhs(clabels[-1]));
  else
    return plots:-contourplot(ee,rng1,rng2,others[]);
  end if;
end proc:
 

contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = 9,
      size=[500,400],
      legendstyle = [location = right],
      legend=true,
      contourlabels=true,
      view=[-2.1..2.1,-2.1..2.1]
);

contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = 17,
      size=[500,400],
      legendstyle = [location = right],
      legend=['contourvalue',$("null",7),'contourvalue',$("null",7),'contourvalue'],
      contourlabels=true,
      view=[-2.1..2.1,-2.1..2.1]
);

# Apparently legend items must be unique, to persist on document re-open.

contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = 11,
      size=[500,400],
      legendstyle = [location = right],
      legend=['contourvalue',seq(cat($(` `,i)),i=2..5),
              'contourvalue',seq(cat($(` `,i)),i=6..9),
              'contourvalue'],
      contourlabels=true,
      view=[-2.1..2.1,-2.1..2.1]
);

contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Green","Red"],
      contours = 8,
      size=[400,450],
      legend=true,
      contourlabels=true
);

contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = 13,
      legend=['contourvalue',$("null",5),'contourvalue',$("null",5),'contourvalue'],
      contourlabels=true
);

(low,high,N):=0.1,7.6,23:
conts:=[seq(low..high*1.01, (high-low)/(N-1))]:
contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = conts,
      legend=['contourvalue',$("null",floor((N-3)/2)),'contourvalue',$("null",ceil((N-3)/2)),'contourvalue'],
      contourlabels=true
);

plots:-display(
  subsindets(contplot((x^2+y^2)^(1/2), x=-2..2, y=-2..2,
                      coloring=["Yellow","Blue"],
                      contours = 7,
                      filledregions),
             specfunc(CURVES),u->NULL),
  contplot((x^2+y^2)^(1/2), x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = 7, #grid=[50,50],
      thickness=0,
      legendstyle = [location=right],
      legend=true),
  size=[600,500],
  view=[-2.1..2.1,-2.1..2.1]
);

plots:-display(
  contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = 5,
      thickness=0, filledregions),
  contplot(x^2+y^2, x=-2..2, y=-2..2,
      coloring=["Yellow","Blue"],
      contours = 5,
      thickness=3,
      legendstyle = [location=right],
      legend=typeset("<=",contourvalue)),
  size=[700,600],
  view=[-2.1..2.1,-2.1..2.1]
);

N:=11:
plots:-display(
  contplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,
      coloring=["Yellow","Blue"],
      contours = [seq(-1+(i-1)*(1-(-1))/(N-1),i=1..N)],
      thickness=3,
      legendstyle = [location=right],
      legend=true),
   plots:-densityplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,
      colorscheme=["zgradient",["Yellow","Blue"],colorspace="RGB"],
      grid=[100,100],
      style=surface, restricttoranges),
   plottools:-line([-2*Pi,-1],[-2*Pi,1],thickness=3,color=white),
   plottools:-line([2*Pi,-1],[2*Pi,1],thickness=3,color=white),
   plottools:-line([-2*Pi,1],[2*Pi,1],thickness=3,color=white),
   plottools:-line([-2*Pi,-1],[2*Pi,-1],thickness=3,color=white),
   size=[600,500]
);

N:=13:
plots:-display(
  contplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,
      coloring=["Yellow","Blue"],
      contours = [seq(-1+(i-1)*(1-(-1))/(N-1),i=1..N)],
      thickness=6,
      legendstyle = [location=right],
      legend=['contourvalue',seq(cat($(` `,i)),i=2..3),
              'contourvalue',seq(cat($(` `,i)),i=5..6),
              'contourvalue',seq(cat($(` `,i)),i=8..9),
              'contourvalue',seq(cat($(` `,i)),i=11..12),
              'contourvalue']),
   plots:-densityplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,
      colorscheme=["zgradient",["Yellow","Blue"],colorspace="RGB"],
      grid=[100,100],
      style=surface, restricttoranges),
   plottools:-line([-2*Pi,-1],[-2*Pi,1],thickness=6,color=white),
   plottools:-line([2*Pi,-1],[2*Pi,1],thickness=6,color=white),
   plottools:-line([-2*Pi,1],[2*Pi,1],thickness=6,color=white),
   plottools:-line([-2*Pi,-1],[2*Pi,-1],thickness=6,color=white),
  size=[600,500]
);