## 5162 Reputation

9 years, 243 days

## Don't understand the question...

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

## Try...

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);
 (1)
 >

## Observation...

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?

## Just realised...

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:
```

## You still have lots of basic syntax erro...

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

 > 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)
 > #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:=<|>: #points:= {seq(B[i],i=0..100)};
 > SIRsys; a(0);g(0);h(0); dsolve([SIRsys[], a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]);
 > p1:=pointplot(L, color=blue, symbol=diamond): display([p1,aplot]);
 >
 >

## Well...

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

## In which case...

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);
 (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);
 (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);
 (3)
 >

## Observations...

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]));
 (1)
 >

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

## More of an observation...

You use the initial condition

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

and the assumption

assuming t>0

These are contradictory - is t=0 within the domain of the solution or not?

If t=0 is within the domain of the solution, then the assumption is wrong. If t=0 is not in the domain of the solution, then the initial condition is (at best) superfluous. So I'm confused! (And maybe Maple is as well?).

So for simple logical consistency, once should have something like

```restart;
#infolevel[pdsolve]:=10;
unassign('z,t,r,u');
lap:=diff(u(r,z,t),r\$2)+ 1/r*diff(u(r,z,t),r)+diff(u(r,z,t),z\$2);
bc:=u(r,0,t)=0,u(r,1,t)=0, u(1,z,t)=0;
ic:=u(r,z,0) = f(r,z);
pdsolve([diff(u(r,z,t),t) = lap,bc,ic],u(r,z,t)) assuming t>=0
```

This still produces the same error :-(

However, the above correction required some thought about how "meaningful" the initial condition actually is. After all if one has a function u(r,z,t), then the function u(r,z,0) is (formally) a function of the two arguments 'r' and 'z' - so writing u(r, z, 0)=f(r, z) is really rather pointless. It contains no useful information. Interestingly, if one removes this meaingless(?) initial condition from your code, then it executes with no error. So the code

```restart;
unassign('z,t,r,u');
lap:=diff(u(r,z,t),r\$2)+ 1/r*diff(u(r,z,t),r)+diff(u(r,z,t),z\$2);
bc:=u(r,0,t)=0,u(r,1,t)=0, u(1,z,t)=0;
ic:=u(r,z,0) = f(r,z);
pdsolve([diff(u(r,z,t),t) = lap,bc],u(r,z,t)) assuming t>=0
```

will return u(r, z, t) = 0 !

This site has a wonderful invention which allows users to upload executable code, rather than some mixture of (non-executable) code+comments

This invention is the big green up-arrow in the Mapleprimes toolbar

I suggest you use it

## Genuinely seeking enlightenment...

Read your post with interest and tried to figure out why your solution was "numerically stable" and the one I originally posted was not.

I note that you made the same (possibly unjustified??) assumptions about boundary conditions which I did. As far as I can tell, the only "significant" difference between our solutions is that you removed the 'spacestep=1" option from the pdsolve() command.

On reflection, looking at the Curvefit() of the function A(x), this seems like a very sensible strategy and I can't think why I didn't do it. Mea Culpa etc.

Can you confirm that the significant difference between our solutions is removal of the "spacestep" option? (Removing the "timestep=0.005" command doesn't seem to make much difference either way).

Reposting my original showing the impact of keeping/removing the "spacestep" option.

 >
 >
 >
 >
 >
 >
 >
 >
 >
 (1)
 >
 >
 >
 >

## Specific to Maple 2015...

I checked the same worksheet in Maple 18, Maple 2015, Maple 2016, Maple 2017, Maple 2018 and Maple 2019

Only Maple 2015 showed the problem

## My debug process...

which may (or may not !) create some enlightenment

Stage1

I could reproduce an "anomaly" on your worksheet dodgy0.mw. The commands

sol1 := dsolve(sys, numeric, method=rkf45);
sol2 := dsolve(sys, numeric, method=ck45);
sol3 := dsolve(sys, numeric, method=rosenbrock);

*ought* to produce the outputs

sol1:=proc(x_rkf45) ... end proc
sol2:=proc(x_ck45) ... end proc
sol3:=proc(x_rosenbrock) ... end proc

but in fact they produce nothing. This is weird! Note that I am running Maple2015.2 on 64-bit Windows 7. So the anomaly is not OS-related

Stage2

Somehow you have some weird (non-printing?) character in your worksheet. So I started with a clean workesheet and retyped the "minimal" example shown in the attached. So this contains no possibility of random non-printing characters.

Even on this minimal woksheet - no output from the command

sol1:=dsolve(sys, numeric);

 > restart;
 > interface(version);
 (1)
 > sys:={ diff(x(t), t)=(2+3*I)*x(t)-(1-I)*t, x(0)=1 };
 (2)
 > sol1:=dsolve(sys, numeric);
 >

Stage3

Quick check on the *same* worksheet loaded in the immediately preceding and succeeding Maple versions, ie Maple 18 and Maple 2016

 > restart;
 > interface(version);
 (1)
 > sys:={ diff(x(t), t)=(2+3*I)*x(t)-(1-I)*t, x(0)=1 };
 (2)
 > sol1:=dsolve(sys, numeric);
 (3)
 >

 > restart;
 > interface(version);
 (1)
 > sys:={ diff(x(t), t)=(2+3*I)*x(t)-(1-I)*t, x(0)=1 };
 (2)
 > sol1:=dsolve(sys, numeric);
 (3)
 >
 >

So at this point we have established that the anomaly is not OS-related and it it is specific to Maple 2015.

At this point further diagnosis gets a bit "random", so don't ask why I'm doing any of this, it is just random stuff I tried

Stage4

Remove the assignment on the dsolve() command and the anomaly disappears! I mean WTF!. And at this stage this site will not display the contents of the following, but take my word for it, just replacing

sol1:=dsolve(sys, numeric);

with

dsolve(sys, numeric);

and as if by magic

proc(x_rkf45) ... end proc

appears

Stage5

I also concur with your observation that replacing 'I' with pretty much anything else in the definition of the ODE also makes this anomaly go away

Conclusions

1. The anomaly is not OS-related
2. The anomaly is specific to Maple 2015
3. The anomaly is somehow related to the occurrence of 'I' in the ODE definition
4. The anomaly is somehow related to whether or not the output of the dsolve() command is assigned to anything
5. The anomaly is much too subtle for me to work out WTF is going on - sorry:-(

## Can't reproduce...

Error, (in dsolve/numeric) 'events' can only be used with real valued IVP/DAE problems

I am doing this on 64-bit Windows 7.

A pretty random suggestion - you had some kind of save error previously. Depending on how you subsequently managed to achieve the save and what you did with filenames afterwards, you may have left a "dodgy" .bak file in your working directory. (These .bak files are produced by Maple's autosave operation). If your current file (events.mw) has the same name as a possibly-corrupted .bak file, this *may* be the source of your problem

If I am honest I always thought that Logic:-Satisfy() was only really suitable for "small" problems - such as the sort of examples which are given on its help page.

I only began to grasp the possibilities of this command a couple of months ago when I came across this solution to the "worls's hardest Sudoku problem"

https://www.maplesoft.com/applications/view.aspx?SID=154483

I had written my own "heuristic-based" Sudoku solver in Maple, but it failed on this example. So, like your code, this solution impressed the hell out of me.

Until I saw the Sudoku example (and now yours), I would have discounted the viability of programmatically generating hundreds (thousands?, more?) of logical conditions and expecting Satisfy() to come with a solution.

I still have no real notion of how "big" one could go but like I said - it is "underadvertised"

 4 5 6 7 8 9 10 Last Page 6 of 147
﻿