tomleslie

5348 Reputation

15 Badges

9 years, 294 days

MaplePrimes Activity


These are replies submitted by tomleslie

@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]);

 

 

 

 

 

Forgot the attachment - see below

U[1, 6]:=(x, theta)->0.03215257166*(sin(-2.350000000 + 9.400000000*x) - 0.1369508410*sinh(-2.350000000 + 9.400000000*x))*cos(6*theta);

proc (x, theta) options operator, arrow; 0.3215257166e-1*(sin(-2.350000000+9.400000000*x)-.1369508410*sinh(-2.350000000+9.400000000*x))*cos(6*theta) end proc

(1)

with(plots):

#
# OP's original plot. Note the command 'cylinderplot()'
# has been deprecated since Maple 8 which was superseded
# in 2003 and is now replaced by plot3d(...coords=cylindrical)
#
  cylinderplot( U[1, 6](x, theta) - 0.1,
                theta = 0..2*Pi,
                x = 0..0.5
              );
  plot3d( U[1, 6](x, theta) - 0.1,
          theta = 0..2*Pi,
          x = 0..0.5,
          coords = cylindrical
        )

 

 

  cylinderplot( 1,
                theta = 0..2*Pi,
                x = 0..0.5
              );
  plot3d( 1,
          theta = 0..2*Pi,
          x = 0..0.5,
          coords = cylindrical
        )

 

 

  cylinderplot( 0.5*x,
                theta = 0..2*Pi,
                x = 0..0.5
              );
  plot3d( 0.5*x,
          theta = 0..2*Pi,
          x = 0..0.5,
          coords = cylindrical
        )

 

 

 

Download xomeplots.mw

@9009134 

You have already obtained a plot of U[1,6] in your original worksheet using the deprecated command cylinderplot().

This should (really) be replaced by plot3d(...., coords = cylindrical) unless you are using Maple 8 (or an earlier version), which was supereded in 2003. The two commands give the same result - see the attached.

You have already obtained these plots. What is wrong with them???

The most common mistake people make when using these commands is to assume that they are producing a plot of h(r, theta) in cylindrical coordinates (r, theta, z), when in fact they produce a plot of r(theta, z)

So, for example

cylinderplot( 1,
                theta = 0..2*Pi,
                x = 0..0.5
              );

will plot a cylinder of radius 1, where the 'height' coordinate varies from 0 to 0.5. Similarly

cylinderplot( 0.5*x,
                theta = 0..2*Pi,
                x = 0..0.5
              );

will produce a "cone", since the radius is computed radius is proportional to the height.

These are illustrated in the attached

the attached, in which I have modified the function so that it returns "FAIL" if t<1.

restart;
f:= t-> `if`(t>=1, combine(int(1/x, x = 1 .. t)), FAIL);
f(10);
f(10.0);
f(1);
f(1.0);
f(0.1);

proc (t) options operator, arrow; `if`(1 <= t, combine(int(1/x, x = 1 .. t)), FAIL) end proc

 

ln(10)

 

2.302585093

 

0

 

0.

 

FAIL

(1)

 


 

Download doInt.mw

Throughout your expression, you use the '^' character to denote exponentiation, except in the last line where you use '**'. Is the latter usage deliberate or a typo?

That the results "look" better if one uses "thicker" lines: suggest changing the doContour() procedure in my previous post to the following

doContour:= proc( nsteps::posint)
                    local j;
                    uses plots, ColorTools:
                    return display
                           ( [ seq
                               ( implicitplot
                                 ( sin(x*y)=-1+2*(j-1)/nsteps,
                                   x = -2 .. 2,
                                   y = -2 .. 2,
                                   axes = framed,
                                   color = WavelengthToColor
                                           ( 450+(j-1)*(625-450)/nsteps,
                                             method="linear"
                                           ),
                                   if   type(j, odd)
                                   then legend=(-1+2*(j-1)/nsteps)
                                   else NULL
                                   fi,
                                   thickness=4,
                                   scaling=constrained,
                                   size=[800,800]
                                 ),
                                 j = 1..nsteps+1
                               )
                             ]
                           );
              end proc:

 

@ogunmiloro 

I have fixed many of these in the attached - but I doubt if I caught them all.

Please read the comments in the attached very very carefully


 

restart:

#
# Statement not terminated in the following - fixed it
#
  B[0]:= 1.7e8:  C[0]:=0:  P[0]:=400: E[0]:=20:
#
# Why does this (commented) Digits command exist?
#
#Digits:=10:

#
# OP never uses 'gamma' in this worksheet - why the
# hell unprotect an inbuilt variable for no purpose???
#
unprotect('gamma'):


 A__h := 10:
 beta__1 := 0.34e-1:
 psi__o := 0.25e-1:
 mu := 0.4e-3:
 theta := .7902:
 alpha := .11:
 k := 0.136e-3:
 z := 0.5e-1:
 delta := .7:
 eta := .134:
 xi := .12:
 beta := 0.1e-1:
 Q1 := 1:
 Q2 := 1:
 h := .1:
#
# The variables 'tau' and 'kappa' are used in the following,
# but neither of these is defined. Gave them a "random"
# values here just to assist debug
#
  tau:=1.0;
  kappa:=2.0
;

1.0

 

2.0

(1)

#Step 1
n:=100:
#
# Fixed the names of all of the following for some kind
# of consistency
#
lambda1[n] := 0: lambda2[n] := 0: lambda3[n] := 0: lambda4[n] := 0: u0[0] := 0:

for i from 0 to n-1 do

     B[i+1] := h*(theta*E[i]+A__h+B[i])/(1+ psi__o*h*P[i])/(P[i]+xi)+mu+u0[i]+beta*C[i];
#
# Following statement is a recursive assignment.
# It was a recursive assignment the last time the
# OP posted some code, and it is still a recursive
# assignment now!!!!
#
# Let's see if I can make this explanation simple enough
# for the brain-dead.
#
# C[i+1] is undefined. Therefore, assigning C[i+1] in
# terms of C[i+1] is doomed to failure.
#
# Changing C[i+1] to C[i] on the rhs is *may* be incorrect
# but at least the statement can now be evaluated!!
#
#     C[i+1] := h*B[i+1](beta*C[i]-psi[o]*P[i])/(P[i]+xi)-(alpha+mu+k)*C[i+1];
#
# I'm also assumin that psi[o] ought to be psi__o, so
# I made this change
#
     C[i+1] := h*B[i+1](beta*C[i]-psi__o*P[i])/(P[i]+xi)-(alpha+mu+k)*C[i];
     P[i+1] := h*z*C[i+1]/(1+delta+eta);
     E[i+1] := h*(alpha*C[i+1]+B[i+1]*u0[i]+E[i])/(1+theta+mu);

#
# The following statement has one more oopening parenthesis (ie '(')
# than closing parenthesis - again basic syntax
#
# I randomly added a closing parenthesis, just to make it syntactically
# correct, but the odds that I got this closing parenthesis in the
# desired location are almost nil!!!!!
#
#   lambda1[n-i-1] := -h*(Q1+(beta*C[i]+(psi__o*P[i]/(P[i]+xi)))*lambda2[n-i]-lambda1[n-i]-mu+u0[i]/(1+beta*C        [i]+P[i]/(P[i]+xi));
#
  lambda1[n-i-1] := -h*(Q1+(beta*C[i]+(psi__o*P[i]/(P[i]+xi)))*lambda2[n-i]-lambda1[n-i]-mu+u0[i])/(1+beta*C[i]+P[i]/(P[i]+xi));
#
# The variable 'kappa' in the following is undefined
#
  lambda2[n-i-1] := -h*(-beta*B[i+1]*lambda1[n-i-1]+Q2-alpha-mu-kappa+z*lambda3[n-i]+alpha*lambda4[n-i]-lambda2[n-i])/(-beta*B[i+1]+1);

  lambda3[n-i-1] := h*((psi__o*B[i+1]*P[i+1]/(P[i+1]+xi)^2+theta)*lambda1[n-i-1]-lambda2[n-i-1]*psi__o*B[i+1]*P[i+1]/(P[i+1]+xi)^2-lambda3[n-i])/(1+delta+eta);

  lambda4[n-i-1] := h*lambda4[n-i]/(1-theta-mu);
#
# The variable 'tau' was undefined, so (obviously) the
# following cannot be fully evaluated
#

     R1:= (lambda1[n-i-1]-lambda4[n-i-1])*B[i+1]/tau;
#
# Since R1 has not been fully evaluated, u0[i+1] cannot
# be fully evaluated either - and so the chaos continues!
#
    
      u0[i+1] := min(1, max(R1, 0));    
    end do:

#
# These two command are assigned to nothing. They are
# are never used for anything - why do they exist??
#
seq(i,i=0..30):
seq(B[i],i=0..30):

with(plots):

with(DEtools):

#
# The quantity 'h', was earlier defined to be 0.1
# Thus h(t)=0.1, h(0)=0.1 and so on. Probably not
# what the OP desired. Even funnier is when the
# OP uses initial conditions such as h(0)=400,
# which of course will evaluate to 0.1=400. This
# is a whole new level of misunderstanding
#
# Obviously this mean that the definition of these
# ODEs is completely incorrect - so it is a complete\
# waste of time trying to solve them!
#

DE1 := diff(a(t),t) = A__h-beta__1*a(t)*g(t)-psi__o*a(t)*h(t)/(h(t)+xi)-mu*a(t)+theta*m(t);
DE2 := diff(g(t),t) = beta__1*a(t)*g(t)-psi__o*a(t)*h(t)/(h(t)+xi)-(alpha+mu+k)*g(t);
DE3 := diff(h(t),t) = z*g(t)-(delta+eta)*h(t);
DE4 := diff(m(t),t) = alpha*g(t)-(theta+mu)*m(t);

SIRsys := [DE1, DE2, DE3, DE4];
#
# The commands DE1plot(), DE2plot(), DE3plot() DE4plot()
# do not exist in Maple so the following will never
# evaluate. OP (probably??) means DEplot(). However with
# the incorrect use of the name 'h' described above,
# even a syntactically correct command will not
# produce anything meaningful. Correct command with
# crap arguments will still fail!!
#

aplot := DE1plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, a = 0.167e9 .. 0.17e9, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, a(t)], thickness = 2, linecolor = red, stepsize = .1);

gplot := DE2plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, g = 0 .. 400, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, g(t)], thickness = 2, linecolor = red, stepsize = .1):

hplot := DE3plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, h = 0 .. 0.1e4, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, h(t)], thickness = 2, linecolor = red, stepsize = .1):

mplot := DE4plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, m = 0 .. 0.1e4, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, m(t)], thickness = 2, linecolor = red, stepsize = .1):

#
# Did no syntax checking beyond this point becuase I got
# really, REALLY BORED
#

L:=<<seq(i,i=0..100)>|<seq(B[i],i=0..100)>>:
#points:= {seq(B[i],i=0..100)};
 

diff(a(t), t) = 10-0.34e-1*a(t)*g(t)-0.1176363636e-1*a(t)+.7902*m(t)

 

diff(g(t), t) = 0.34e-1*a(t)*g(t)-0.1136363636e-1*a(t)-.110536*g(t)

 

0 = 0.5e-1*g(t)-0.834e-1

 

diff(m(t), t) = .11*g(t)-.7906*m(t)

 

[diff(a(t), t) = 10-0.34e-1*a(t)*g(t)-0.1176363636e-1*a(t)+.7902*m(t), diff(g(t), t) = 0.34e-1*a(t)*g(t)-0.1136363636e-1*a(t)-.110536*g(t), 0 = 0.5e-1*g(t)-0.834e-1, diff(m(t), t) = .11*g(t)-.7906*m(t)]

 

Error, (in DEtools/DEplot) vars must be declared as a list, e.g. [x(t),y(t),...]

 

SIRsys; a(0);g(0);h(0);

dsolve([SIRsys[], a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]);

[diff(a(t), t) = 10-0.34e-1*a(t)*g(t)-0.1176363636e-1*a(t)+.7902*m(t), diff(g(t), t) = 0.34e-1*a(t)*g(t)-0.1136363636e-1*a(t)-.110536*g(t), 0 = 0.5e-1*g(t)-0.834e-1, diff(m(t), t) = .11*g(t)-.7906*m(t)]

 

a(0)

 

g(0)

 

.1

 

Error, (in dsolve) found the following equations not depending on the unknowns of the input system: {1/10 = 400}

 

p1:=pointplot(L, color=blue, symbol=diamond):
display([p1,aplot]);

Error, (in plots:-pointplot) points cannot be converted to floating-point values

 

Error, (in plots:-display) expecting plot structures but received: [p1]

 

J := `<|>`(`<,>`(seq(i, i = 0 .. 100)), `<,>`(seq(C[i], i = 0 .. 100))); p2 := pointplot(J, color = blue, symbol = soliddiamond); display([p2, gplot])

J := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

Error, (in plots:-pointplot) points cannot be converted to floating-point values

 

Error, (in plots:-display) expecting plot structures but received: [p2]

 

NULL


 

Download corrections2.mw

@saemi123 

the following code will get a lot closer to an answer

Rb:=13;
X:=x->(Rb+s[0](x))*cos(x)-Ds[0]*sin(x);

Y:=x->(Rb+s[0](x))*sin(x)+Ds[0]*cos(x);

t:=evalf(Pi);

seq(X(phi_0),phi_0=0..t,0.1);
seq(Y(phi_0),phi_0=0..t,0.1);

It wont produce complete answers, because the functions X() and Y() which you desire require calculation of the (indexed?) function s[0]() which is undefined, and the (indexed?) variable Ds[0], which is also undefined

 

 

 

@Rouben Rostamian  

The OP's solution is as given in the attached.

My comment about the boundary conditions arose mainly because of the first column in the output matrix, which corresponds to u(x=0..20, t=0), where u(0,0) is 1, and all other values are Float(undefined).

Since this first column corresponds to u(x,0), and there is a supplied boundary condition u(x,0)=0, I rather expected Maple's numerical solution to agree with the boundary condition along this boundary. Strangely (for me) a defined value (ie 1) is obtained at u(0,0). To me this seems odd, since the supplied boundary conditions mean that u(x,0) is "well-defined" when x>0 and "ill-defined" when x=0. So I'm stil a little nervous about this solution

  restart:
  interface(rtablesize=10);

[10, 10]

(1)

#
# Define pde, and bcs/ics, then solve it
#
  eq1:= diff(u(x,t),t)=-4 *diff(u(x,t),x):
  cond:= {u(0,t)=sin(.5*Pi), u(x,0)=0}: # conds are inconsistent
                                        # what is u(0,0)???
  sol:= pdsolve(eq1,cond, numeric, time=t, range=0..20);NULL

_m707119744

(2)

#
# Generate a 3D plot of the solution
#
  pp:=sol:-plot3d(t=0..20, x=0..20);

 

#
# Generate a 21x21 matrix of results.
# Columns correspond to t=0..20 in steps of 1
# Rows correspond to x=0..20 in steps of 1
# So, for example
#
#   opM[3,3]= u(2,2)
#
  opM:= Matrix
        ( 21,
          21,
          (i,j)-> rhs~(sol:-value(u(x,t))(i-1,j-1))[3]
        );
#
# Export results matrix to a file. OP will need to change
# file path/name to something appropriate for his/her
# installation/OS
#
  ExportMatrix("C:/Users/TomLeslie/Desktop/opRes.dat", opM);

_rtable[18446744074373019398]

 

8514

(3)

 

NULL

Download pdeSol2.mw

@nm 

The Laplacian of a time-dependent function in cylindrical polars (with no theta-dependence) can be simply obtained with no "substitutions" required, which is probably how I would prefer to do it.

See the attached

restart;
pde:=diff(u(r,z,t),t)= expand(VectorCalculus:-Laplacian(f(r, z, t), cylindrical[r, theta,z]));

diff(u(r, z, t), t) = (diff(f(r, z, t), r))/r+diff(diff(f(r, z, t), r), r)+diff(diff(f(r, z, t), z), z)

(1)

 

Download lapl.mw

The whole t>0, t>=0 thing.

When debugging Maple errors, I think it is good practice to remove any potential ambiguities or inconsistencies from the code generating the error. In this case, removing the inconsistency didn't affect the error, but as a general technique, I think this approach is still worthwhile.

 

 

I still suspect that your initial condition

ic:=u(r,z,0) = f(r,z)

where f(r,z) is undefined is going to cause all sorts of weird problems - because it isn't really a meaningful inital condition: it provides no useful information

As I stated in my first post, removing this initial condition makes the Error go away and a solution is obtained

5 6 7 8 9 10 11 Last Page 7 of 149