4862 Reputation

17 Badges

6 years, 362 days

MaplePrimes Activity

These are replies submitted by mmcdara

@acer @Carl Love 

Great thanks to both of you for this detailed discussion.
To summarize and clarify a few questions :

  • acer/  "but I interpreted the original question as being more about how to write the procedure so that the calls to H could be written with just the short form of the name":
    Absolutely: I wrote "f" first, observed it didn't worked as expected, then tried "g", became happy it worked correctly but remained upset for the syntax orthopoly:-H(...) is obviously heavier than just "uses orthopoly;..... H(...)"
  • acer/  thanks for the solutions you provided in your first answer
  • Carl Love/  about the "semantics of :-" ... I didn't know about this way to access an entry of a table. If I correctly understand what you said, it is some kind of undocummented feature ?
    It has one funny consequence: for instance (attached file), the syntax T:-N (where T is a table and N one if its indices) avoids using eval to visualize T[N] if N is itself a table.
  • Carl Love/  Thanks for having explained "why procedure g works" 
  • both of you/  I understand I should use the "correct coding" acer provided while avoiding the "T:-N" syntax pointed by Carl 


Thanks for your help and your


Thanks for the advice.
I'll do the test since I get back from my holidays.
I'll see you again around 9-10 June


@Carl Love 

I copy-and-pasted the last(?) code.
I get an error.
I think it's because I'm at home with Maple 2015.2 (no Iterator package ; Maple 2016 and Maple 2018 at the office) a

I'm in vacation and I won't be able to do some tests of performances within one week.
I'll see you around 10 June.


Unfortunately I'm in vacation for a week and I won't have access to Maple in 12 hours.
I'm going to try to extract the "real time" by doing a few things at home.
A point is: the code Carl Love provided me  is a module including a procedure with a remember option (the [non optimal] procedure I had written myself doesn't have this option and the cpu times assessed by the two procedures are roughly equal).
I will do some supplementary tests and I will keep you informed about their results as soon as possible.




Thanks for the precisions.

By the way, I have an other thread open with Carl Love as the main contributor.
Catrl wrote different procedures to adress the same problem an I did some comparisons of their performances (mainly the cpu time).

Before your previous answer here ("use option iterations")  I used the procedure below to assess the mean cpu time over 100 runs of the same code.

Procedure 1 :
 T := NULL;
for r from 1 to 100 do
   T := T, CodeTools[Usage](Code(...), output=cputime)  # from your answer to the "other" thread
end do:

Once you made me discover the option "iterations" I wrote this second procedure.

Procedure 2 : 
 CodeTools[Usage](Code(...), output=cputime, iterations=10

Here are the performances I got.

  indicator                Proc 1                   Proc2            

 cpu time/run             16.9 s                     32.75 s             

memory size            0.95 GiB                  0.95 GiB

I repeated this twice and obtained the same ratio of 2 between the cpu times.
For proc1 the 10 times range between 15.49 s and 18.08 s

Do you have any idea about where this difference could come from ?
Do you think it is "intrinsic" to the CodeTools[Usage] procedure or it comes from the way Code(...) is written ?

Thanks for answer

PS : this reply is close to the one I sent to Carl Love in the other thread



But the reference provided by Mariusz is more complete and I used it.
I can't obtain a solution for the system described in the reference paper by using pdsolve.

There remain some points to fix:

1/ For NS equations, the pression is only defined up to an arbitrary constant. I think P(0,0)=0 is a reasonable boundary condition (bc).

2/ bc n° 10 and 11 cannot be handle with pdsolve.

3/ The writting P(z) xauses problems: the system has 2 unknowns u and P and each of them has to be defined in terms of r and z : u --> u(r, z) and P --> P(r,z).
Of course this induces some complexification s because we need to account for the equation diff(P(r,z), r) = 0


I think it will not be possible to obtain a solution by using "blindly" pdsolve.
A solution strategy (for instance solve equation 5 and plug it in equation 7) seems necessary.
It will probably take some time to do that and, very honestly, I'm not sure I will be capable of that in a reasonnable time.

I will keep trying to solve this problem, but do not count on me too much.


@Mariusz Iwaniuk 

Be happy !
It works pretty well with Maple 2015.2 and Maple 2016.1 too

@Mariusz Iwaniuk 

Damn !
I'm just going to complete my answer to waseem.

Thanks for the reference.


I can't find a free version of the article but I'm suppose that you study the stationnary laminar flow of a non-newtonian (incompressible) fluid along a duct of radius r(z) ?
In cylindrical coordinates (with adhoc auxiliary assumptions) stationnary Navier Stokes equations for a newtonian fluid 
write for a duct of constant section (Poiseuille's flow)
(1/r) * diff( mu * r * diff(u(r,z), r), r) = diff(P(z), z)

By comparing this equation to your equation (1) it comes to me 
1/ that the right hand side of (1) should be written diff(P(z), z) (pressure gradient along the axis of the "duct" (artery, vein).
    ... or at least P(z)
2/ u should be written as a function of r and z

But ...
If the radius of the duct is a function of z, then the field of the axial velocity u(r, z) cannot be the same in every cross section of the duct. Then :
1/ either the duct has a constant radius, in which case your boundary conditions (3) and (6) should be written "at r = R" (where R is some strictly positive value).
2/or either some terme invoking diff(u(r,z), k) should be present on the left hand side of (1) (unles adhoc assumptions)

Last but not least ...
Isothermal NS (Navier-Stokes) equations consist in two conservations equations, one for mass conservation, the other por impulsion conservation.
Equation (1) is obviously the impulsion condition
Without mass conservation equation you can't obtain a solution and I doubt the authors of the reference you gave were able to ...

On your second reply you say "they got the results in terms of bessels".
Viewed from a distance, if we omit the term "1" in the most inner part of (1), you have a bilaplacien operator of the type encounterd in vibration of circular plates ... which as a solution in terms of the Airy function

Could you provide me the correct set of equations written in Maple ?


@Carl Love 

Your Question mentions "base 10" three times. This problem has nothing whatsoever to do with base 10.
Absolutly, we could imagine the same matrices expressed in any other base, hexadecimal base for instance.
"do c = a+b just as if a and b were numbers in base 10."
Right, it was only to say that we are not doing all the operations in Z/pZ (p being the base) but that some steps are done in an another base (in fact the base in which the elements of the matrices are written)

See also my second Reply to Tom below.
Apologies, I royally missed it.

The main problem with multiple contributions from multiple contributors is to reconstruct correctly the timing of the events and to follow the  thread you think is the most advanced and which is the most promising.
So I recognize that I let go Tom's contributions or the exchanges you had directly with him.



apologies for not having reply sooner to your post.
I was studying the last contributions CarlLove had sent me and, to be honnest, I thought that only Carl was racking his brain on the subject ...

Considering your code: yes,  your answer (906) is correct, Carls' wans't.
But this was the state 8 hours ago.
Meanwhile I went to sleep and this morning I discover the last Carl's post (3 h ago).
(Please have a look to the answer I sent him 2 h ago)

I observed the last Carl's code had corrected an error (probably due to a vagueness from my part about the" carry/not carry"  rule (Carl"'s question 13 h ago, my answer 12 h ago) and maybe some comparison with yours.

In any event, thank you for your contribution.


@Carl Love 

For information.

Please find attached the file containing the performances of the different procedures used to build a "Base=5, K=2" matrix

@Carl Love 


I observe that, at the same time, you have corrected your unitial set of procedures to account for bases greater than 4.
I had observed last night (its 8:30 am here) that you first code didn't returnd the correct Matrix for base 7 but it was two late to warn you about that.

On top of that your module uses a lot of technical stuff I didn't know about. It will probably take me some time to understand it correctly, but I think I will learn a lot from it at the end



@Carl Love 

If you want to go the polynomial route, use PolynomialTools:-FromCoefficientList rather than series.
No, not necessarily, it was just my initial attempt but your procedures suit me perfectly.

If you have foreknowledge of the overall maximum list length needed
I don't have any prior idea of the this value. I want to study those special matrices of dimansion base^K where K is some "smal" positive number, lets's say 3 or 4, no more.
Then you can consider that the overall maximul list length doesn't exceed 4

 It appears that for your work that there is such foreknowledge.
Right. These matrices probably first appeared in game theory. One can also relate them to some cryptography problems: In this later case I guess one can imagine huge matrices (K large), but I'm still in a prospective stage where I'm looking to these matrices as potentially usefull tolls. 
To investigate this the restriction K<=4 is sufficient. 
You can find a situation where the matrix Matrix(3^2, 3^3, (i,j) -> Nim(i-1, j-1, 3)) appears in

I need for you to clarify if a "carry" operation for the addition is needed when the base is greater than 4.
Right again
You raise two points: the one of the carry and, behint it, the value of the base.
Second point first: the base has to be a prime number.
First point: there is effectively an issue when the base is greater or equal to 5... but the additions have to be done withour carry. 
Here is an example with base 7, A=13 and B=27
A = [1, 6] mod 7    (convert(A, base, 7)=[6, 1])
B = [3, 6] mod 7    (convert(A, base, 7)=[6, 3])
A+B = [4, 12]         (bitwise "no carry addition" in base 10)
A+B  mod 7 = [4, 5]
A+B = 33   in base 10

WATCH OUT : I focuse only on bases 2 and 3, then I'm not sure that my procedures return the good matrices for bases greater or equal to 5 !

I hope I have answered to your probing questions.

For information
The alternative ways to generate such matrices start from the "generator" G := Matrix(N, N, (i,j) -> Nim(i-1, j-1, N)) (N being the base, 3 in my initial post) and by using the KroneckerProduct. Their complexities can be easily established.
This product, also termed sometimes "tensor product", makes a (still unformal) connection between these matrices and  tensor product of graphs.

@Carl Love 

(I'm sand15, now out of the office and at home) 

Yes, it's ok.
Here is the matrix I want to construct :

N := 3: # base
K := 3:
M := Matrix(N^K, N^K, (i,j) -> Nim(i-1, j-1, N))

Even if the "plus" operation is the standard way to constructthis matrix M (for different values of K), I was able to derive different alernate ways to build M.
For instance the procedures A and B below 

A := proc(n,k)
  local G, B, un, i:
  uses LinearAlgebra:
  G := Matrix(n, n, (i,j) -> (i-1)+(j-1) mod n):
  B := copy(G):
  un := m -> Matrix(m, m, 1):
  for i from 1 to k do
    B := KroneckerProduct(un(n^i), G) +~ KroneckerProduct(n *~ B, un(n))
  end do:
end proc;

Or again, after having resolved the recurrence relation 
B := proc(n, k)
  local G, un :
  uses LinearAlgebra:
  G  := Matrix(n, n, (i,j) -> (i-1)+(j-1) mod n):
  un := m -> Matrix(m, m, 1):
  add(n^p *~ KroneckerProduct(KroneckerProduct(un(n^(k-p)), G), un(n^p)), p=0..k)
end proc;

A(N,K) = B(N, K) = Matrix(N^K, N^K, (i,j) -> Nim(i-1, j-1, N))

For information, my coding for the direct building of the matrix  was 
C := proc(n, k)
  local a, b, sa, sb, pa, pb, pc, c:
  M := n^k:
  U := Matrix(M, M);
  for p from 0 to M-1 do
    for q from 0 to M-1 do
      a   := convert(p, base, n);
      b   := convert(q, base, n);
      sa  := gfun[listtoseries](a, x, 'ogf');
      sb  := gfun[listtoseries](b, x, 'ogf');
      pa := convert(sa, polynom):
      pb := convert(sb, polynom):
      pc  := pa + pb mod n;
      c   := subs(x=n, pc);
      U[p+1,q+1] := c;
      U[q+1,p+1] := c;
    end do:
  end do:
end proc;

And here again C(N, K) is equal to M

I wasn't happy of the procedure C which seems to make a moutain of a molehill ...  
Your solution seems clever. Nevertheless I'm not sure we can answer the question without making a detour via polynoms (your B10 procedure seems rather close to my "pc := ...; and c:= ...; above).

I will finally adopt your solution.

Thank you for your time



First 95 96 97 98 99 100 101 Last Page 97 of 106