tomleslie

4741 Reputation

15 Badges

9 years, 154 days

MaplePrimes Activity


These are replies submitted by tomleslie

in the Mapleprimes toolbar to upload a worksheet with your calculation

That way the problem can be sensibly invetigated

@Carl Love 

for the wrong question! It was intended for

https://www.mapleprimes.com/questions/227415-How-To-Convert-The-Tekplot-File?reply=answer

because I only have the last six major versions of Maple installed on this machine.

Code I supplied works fine on all of them, but that only goes back to Maple 18!

 

Maple uses various methods to "restrict" the display of "output", amongst which are

  1. If a command is terminated with a colon ( ie ':'), the command will be executed but no output will be displayed
  2. If a command is terminated with a semi-colon ( ie ';'), the command will be executed and output will be displayed
  3. By default Maple will not display of matrices bigger than 10X10. This can be changed by an appropriate setting of interface(rtablesize).
  4. If you want to display the complete contents of a table (called say T), then you need to use T()

In the attached, I have terminated commands with ";" rather than ":" (see 1 above), so the output of more commands will be displayed. I have left colons terminating the 'for' loops - in general you don't want to see the output of every iteration of a loop, becuase this can fill a lot of "screens"

I have also changed set interface(rtablesize=15), so anything with 15 (or fewer) entries in either dimension ought to display completely

If you  download/execute the attached, then it should display (more-or-less) as it displays here

  restart:
  interface(rtablesize = 15):
  with(LinearAlgebra):

  A:= 8;
  B:= 5;
  q:= 0.4;
  p:= 0.2;
  r:= 1 - p - q;
  dimP:= A + B + 1;

8

 

5

 

.4

 

.2

 

.4

 

14

(1)

  P:= Matrix(dimP, dimP, [0 $ dimP*dimP]);
  P[1, 1]:= 1;
  P[1, 2]:= 0;
  P[dimP, dimP]:= 1;
  P[dimP, dimP - 1]:= 0;
  for i from 2 to dimP - 1 do
      P[i, i - 1]:= q;
      P[i, i]:= r;
      P[i, i + 1]:= p;
  end do:

Matrix(14, 14, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (1, 11) = 0, (1, 12) = 0, (1, 13) = 0, (1, 14) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 0, (2, 11) = 0, (2, 12) = 0, (2, 13) = 0, (2, 14) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0, (3, 10) = 0, (3, 11) = 0, (3, 12) = 0, (3, 13) = 0, (3, 14) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0, (4, 11) = 0, (4, 12) = 0, (4, 13) = 0, (4, 14) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (5, 9) = 0, (5, 10) = 0, (5, 11) = 0, (5, 12) = 0, (5, 13) = 0, (5, 14) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (6, 9) = 0, (6, 10) = 0, (6, 11) = 0, (6, 12) = 0, (6, 13) = 0, (6, 14) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (7, 9) = 0, (7, 10) = 0, (7, 11) = 0, (7, 12) = 0, (7, 13) = 0, (7, 14) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0, (8, 9) = 0, (8, 10) = 0, (8, 11) = 0, (8, 12) = 0, (8, 13) = 0, (8, 14) = 0, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = 0, (9, 6) = 0, (9, 7) = 0, (9, 8) = 0, (9, 9) = 0, (9, 10) = 0, (9, 11) = 0, (9, 12) = 0, (9, 13) = 0, (9, 14) = 0, (10, 1) = 0, (10, 2) = 0, (10, 3) = 0, (10, 4) = 0, (10, 5) = 0, (10, 6) = 0, (10, 7) = 0, (10, 8) = 0, (10, 9) = 0, (10, 10) = 0, (10, 11) = 0, (10, 12) = 0, (10, 13) = 0, (10, 14) = 0, (11, 1) = 0, (11, 2) = 0, (11, 3) = 0, (11, 4) = 0, (11, 5) = 0, (11, 6) = 0, (11, 7) = 0, (11, 8) = 0, (11, 9) = 0, (11, 10) = 0, (11, 11) = 0, (11, 12) = 0, (11, 13) = 0, (11, 14) = 0, (12, 1) = 0, (12, 2) = 0, (12, 3) = 0, (12, 4) = 0, (12, 5) = 0, (12, 6) = 0, (12, 7) = 0, (12, 8) = 0, (12, 9) = 0, (12, 10) = 0, (12, 11) = 0, (12, 12) = 0, (12, 13) = 0, (12, 14) = 0, (13, 1) = 0, (13, 2) = 0, (13, 3) = 0, (13, 4) = 0, (13, 5) = 0, (13, 6) = 0, (13, 7) = 0, (13, 8) = 0, (13, 9) = 0, (13, 10) = 0, (13, 11) = 0, (13, 12) = 0, (13, 13) = 0, (13, 14) = 0, (14, 1) = 0, (14, 2) = 0, (14, 3) = 0, (14, 4) = 0, (14, 5) = 0, (14, 6) = 0, (14, 7) = 0, (14, 8) = 0, (14, 9) = 0, (14, 10) = 0, (14, 11) = 0, (14, 12) = 0, (14, 13) = 0, (14, 14) = 0})

 

1

 

0

 

1

 

0

(2)

  p0:= Matrix(dimP, 1, [0 $ dimP]);
  p0[A + 1, 1]:= 1;
  pV[0]:= p0;
  PT:= Transpose(P);

Vector(14, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0})

 

p0[9, 1] := 1

 

Vector(14, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 1, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0})

 

Matrix(%id = 18446744074386657518)

(3)

  for n to 200 do
      pV[n]:= PT . (pV[n - 1]);
  end do:

  map(x -> evalf(x, 3), Transpose(pV[5]));

Vector[row](14, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.102e-1, (5) = 0.512e-1, (6) = .128, (7) = .205, (8) = .230, (9) = .189, (10) = .115, (11) = 0.512e-1, (12) = 0.160e-1, (13) = 0.320e-2, (14) = 0.320e-3})

(4)

 

 

NULL


 

Download mat3.mw

@Rouben Rostamian

As well as the sign of the second derivative, you also flipped the sign of the exponent term in the definition of 'f'. Deliberate?

 

Is this worksheet using 'j' to represent the square root of -1, perhaps set up using an 'init' file?

Maple's default is I

@minhthien2016 

'sol' is a sequence, so you need nops([sol]) to get the number of entries.

Alternatively you can change the solve() command so that it returns a list, rather than a sequence just by writing

sol:=[solve(......)];

in which case nops(sol) will work.

@Rouben Rostamian  

If you incorporate the option 'useint' in the dsolve() command, then Maple will perform the necessary integration for this example.

I am amused by the explanation for this behaviour in the help page at ?dsolve/details, quoted below (with emphasis added)

To request the use of int instead of the default integrator during the solving process. This option makes dsolve attempt computing any integral present in the solution, regardless of the convenience of doing that. This option is of use when the heuristic procedure for deciding about the convenience of performing an integral does not lead to the simpler ODE solution.

I have no idea what the "heuristic procedure for deciding" mentioned in the above actually covers. However I have noted that *sometimes*, when dsolve|() returns an unevaluated integral, the option 'useint' changes this behaviour

An alternative method to "force" evaluation of integrals returned by dsolve() is to use

dsolve({g, v[3](20)=0}, v[3](t));
value(rhs(%));

Since this amounts to evaluating an "inert" integral (ie value(Int(....) ), I assume that this is equivaalint to the 'useint' option

 

@BjaMaple 

converting a "vector" to a "set" is a high-risk activity for several reasons. The two most obvious are

  1. vectors can have repeated entries - so <1,1,1> is absolutely fine as a vector. However it makes no mathematical sense for a set to have repeated entries: so if  you convert the vector <1,1,1> to a set, you will get {1} - which may (or may not) be what you want
  2. vectors have the concept of "order" - so the vectors <1,2,3> and <3,2,1> are two different entities. However sets have no concept of "order". As a matter of convenience Maple stores/displays sets in "lexicographic" order. So converting either of these examples to a set will result in {1,2,3}, and any distinction between them will be lost.

Possibly you might want to convert vectors to lists? This will at least preserve repeated entries and "order"

@ActiveUser 

The keys/indexes for a table have to be unique,

The entries/values associated with a particular key/index can be pretty much anything (repeated as many times as you want)..

To progress any further with this issue, you are going to haev to provide a detailed example of your problem

@9009134 

how was this figure produced?

  1. It could be produced by plotting a function f(x,y,z)=0 for cartesian coordinates (x,y,z) - although in this case I would expect the function f() to be rather complicated
  2. It could be produced by plotting a function f(p, q, r)=0 for some "unspecified" coordinate system defined by the coordinate transformations p=P(x,y,z), q=Q(x,y,z), r=R(x,y,z). In this case the function f() might be "simpler", at the expense of having "more complicated" coordinate transformation functions P(), Q(), R()
  3. You really will have to provide some clue - because I'm not psychic!!

I notice that in an earlier post, you used the following definition

The toroidal coordinate has this equation.

 

X := (R+r*cos(theta))*sin(phi); Y := (R+r*cos(theta))*cos(phi); Z := sin(theta)

Now you can use this definition to plot various, more-or-less, toroidal shapes. In the attached I show

  1. the simplest possible torus, obtained by setting both the 'major' radius (ie R) and 'minor' radius (ie r) in the above definition to be constants
  2. a couple of more complicated toroidal structures obtained by setting the 'major' radius (ie R) to be a constant, and defining the 'minor' radius (ie r) to be more-or less random functions of the variables theta and phi
  3. I have also avoided the use of variable names 'x', 'y' , 'z', 'r', 'R', 'theta', phi' because these are just names, Just because you use the variable name 'theta' (say) this does not, in general, by some magical process imply a polar angle. This means that you will have to read the attached carefully, in order to figure out what coordinatesystem is being used where

#
# The simplest torus possible
#
  restart;
  Z:=4:
  s:=(p, q)->1:
  plot3d( [ ( Z+s(a,b)*cos(a))*sin(b),
            ( Z+s(a,b)*cos(a))*cos(b),
              s(a,b)*sin(a)
          ],
          a=-Pi..Pi,
          b=-Pi..Pi,
          scaling=constrained,
          axes=none
        );

 

#
# A "torus" with a variable minor radius
#
  s:=(p, q)-> 1+cos(4*p)*sin(4*q)/2:
  plot3d( [ ( Z+s(a,b)*cos(a))*sin(b),
            ( Z+s(a,b)*cos(a))*cos(b),
              s(a,b)*sin(a)
          ],
          a=-Pi..Pi,
          b=-Pi..Pi,
          scaling=constrained,
          axes=none
        );

 

#
# Another "torus" with a variable minor radius
#
  s:=(p, q)-> 2+cos(8*q)/2-cos(2*p)/2:
  plot3d( [ ( Z+s(a,b)*cos(a))*sin(b),
            ( Z+s(a,b)*cos(a))*cos(b),
              s(a,b)*sin(a)
          ],
          a=-Pi..Pi,
          b=-Pi..Pi,
          scaling=constrained,
          axes=none
        );

 

 

Download moreTori.mw

 

My previous comment about floating point equality verification in the procedure 'check' was incorrect. The arguments are always integers, so the equality checking is fine. I have edited my previous post to point out my stupidity

There *seems*  to be a typo in the procedure 'check'

check := proc(a, b, c, d, index)
              local flag, i, j;
              flag := 1;
              for i to index do
                  if   a = c[i]
                  then flag := 0;
                       break;
                  else for j to index do
                           if   b = d[i]
                           then flag := 0;
                                break;
                           end if;
                       end do;
                  end if; 
              end do;
              return flag;
         end proc

Nothing in the highlighted loop depends on the associated  loop index 'j': so the loop performs an absolutely identical calculation 'j' times which is completely pointless. I'm going to guess that what you meant was the following

check := proc(a, b, c, d, index)
              local flag, i, j;
              flag := 1;
              for i to index do
                  if   a = c[i]
                  then flag := 0;
                       break;
                  else for j to index do
                           if   b = d[j]
                           then flag := 0;
                                break;
                           end ifOK;
                       end do;
                  end if;
              end do;
              return flag;
         end proc

Note the change of index variable on d[..]

Even with this change there is still a fundamental issue about when/whether the following loop within the procedure 'f1' will ever exit

while(d>0.1 or flag=1) do
     rx:=rand(1..n):
     x[i]:=rx():
     ry:=rand(x[i]..n):
     y[i]:=ry():
     flag:=check(x[i],y[i],xdumb,ydumb,dumbindex):
     if   flag=1
     then dumbindex:=dumbindex+1:
          xdumb[dumbindex]:=x[i]:
          ydumb[dumbindex]:=y[i]:
     fi:
     p,d:=f(x[i],y[i])
end

Monitoring the values of the parameters 'd' and 'flag' using the debugger, all I can say is that the value of 'd' is rarely (if ever) <=0.1, and the value of 'flag' chenges pretty randomly. Thus the probability that

d<=0.1 and flag=0

for the loop to exit seems pretty remote - how many iterations do you expect????

I tried restricting the number of iterations in the above 'while' loop to 100, so with the enclosing 'for' loop running 50 times, a total of 50,000 calculations, but Igot bored and terminated this after about 30minutes execution time.

I only have three suggestions at this point

  1. Change the 'logic' of the code so that this 'while' loop does not run for an indeterminate (possibility infinite) time
  2. Rewrite the same "logic" in an efficient way - the current implementation is very poor, and I'd guess(!?) that could generate a speed-up of around 10X. Although you should be aware that reducing an infinite execution time by a factor of 10 still results in an infinite execution time!
  3. Leave it running - you never know - it might(?) finish eventually!

 

Edited: see commented in red text

which I forgot in my previous reply

Ignore the following - I misread the output of the Maple debugger: the 'check' procedure only receives integer arguments, so the equality checking is fine. Mea Culpa etc

The procedure 'check' returns an integer value (1 or 0) based on checking the equality of two floating point numbers. This is really bad programming practice: don't ever do it. If a and b are floats, then don't test a=b, test abs(a-b)<tol where 'tol' represents an appropriate "tolerance" for an equality check - say 10-6 (or 10-9 if you are feeling brave), depending on the setting of Digits.

@cencen_cj 

I have changed the parentheses as described in my earlier response. If this is incorrect please explain.

Your next problem is the use of the same names for indexed and unindexed quantities - this is never a good idea. In the procedure 'f1' you have the lines

vx := rand(-n .. n);
vx[i] := vx();
vy := rand(-n .. n);
vy[i] := vy();

where the indexed names are the same as the procedure names for the random number generators. If you change this to

rvx := rand(-n .. n);
vx[i] := rvx();
rvy := rand(-n .. n);
vy[i] := rvy();

then the recursion error goes away.

However the news is not all good - the code exeutes but never terminates!. The problem seems(?)to be in the "while" loop

while 0.1 < d or flag = 1 do
    rx := rand(1 .. n);
    x[i] := rx();
    ry := rand(x[i] .. n);
    y[i] := ry();
    flag := check(x[i], y[i], xdumb, ydumb, dumbindex);
    if flag = 1 then
        dumbindex := dumbindex + 1;
        xdumb[dumbindex] := x[i];
        ydumb[dumbindex] := y[i];
    end if;
    p, d := f(x[i], y[i]);
end do

where the two conditions (d<=0.1, and flag=0) never seem to occur simultaneously - they do occur separately.

It's past my bedtime now. If I get the time, I may be able to investigate further tomorrow

Since you do not supply the files "a1.txt", "a2.txt", "a3.txt", "a1bar.txt", "a2bar.txt", "a3bar.txt", "servrate.txt", your code is non-executable.

It is pretty much impossible to debug your code without being able to execute it!

Simply "reading" your code (a lousy way to attempt debug) there are a couple of places where I suspect problems may arise. In the procedure 'f1', you have the lines

vx[i] := w*vx[i] + 2*rp*[xl[i] - x[i]] + 2*rd*[xg - x[i]];
x[i] := x[i] + vx[i];
vy[i] := w*vy[i] + 2*rp*[yl[i] - y[i]] + 2*rd*[yg - y[i]];

The use of square brackets (ie '[]' ) for the term groupings  [xl[i] - x[i]] , [xg - x[i]], [yl[i] - y[i]], [yg - y[i]] looks very doubtful. This will result in the addition of a scalar to a (single-element?) list. Did you mean

vx[i] := w*vx[i] + 2*rp*(xl[i] - x[i]) + 2*rd*(xg - x[i]);
x[i] := x[i] + vx[i];
vy[i] := w*vy[i] + 2*rp*(yl[i] - y[i) + 2*rd*(yg - y[i]);

 

 

 

 

 

1 2 3 4 5 6 7 Last Page 1 of 143