tomleslie

5600 Reputation

17 Badges

10 years, 9 days

MaplePrimes Activity


These are answers submitted by tomleslie

A cou00le are shown in the attached worksheet

3dLine.mw

 

If you change the final command in your worksheet to specify the indeterminate function ( ie follow the instruction given in the error message), as in

pdsolve(sol, w(x,y));

then Maple will (correctly!) return

w(x, y) = 0.

Note that this satisfies all of your boundary conditions, and your PDE - so is definitely a solution for the problem

Somehow I have the feeling that this is not what you want!

In this case you are being warned that you should only use solve() to solve for names or functions, ie not expressions such as

a[1]*a[2]

presumably becuase the requested operation is not guaranteed to work correctly!

If you really, really want to, then you can turn off warnings (I do not recommend doing this). But if you precede your code with

interface(warnlevel=0);

then all warnings will be turned off. I would reset this to the default immendiately after your code, with

interface(warnlevel=3);

I suggest that you also check some of the terms in your definition of system3d: in particular

a[1](a[1]),

a[5](a[5]),

a[9](a[9])

Are these really what you meant???

 

This wold be OK if 'x' were guaranteed to be of type 'algebraic', Now what would you expect to happen if 'x' is not of type algebraic??

If you want to know what type/algebraic actually means, the check out the help page at ?type/algebraic

There are many Maple types which are not 'algebraic'!

L[1..5] returns an array whoses indices are 1..5 and whose values are L[1], L[2], L[3], L[4], L[5]

L[3..6] returns an array whose indices are 3 ..6 and whose values are L[3], L[4], L[5], L[6].

You seem to expect that L[3..6] would return an array whose indices are 1..4 and whose values are L[3], L[4], L[5], L[6]

So the obvious question you have to ask yourself is: when you ask for L[3..6] - what are the indices of the returned array?? Are they 3..6 or 1..4???

Now I could (if I had to) come up with argument in favour of either: but at some point, somewhere in history, somebody decided that L[3..6] would return an array whose indices are 3..6 (ie not 1..4), and then decide on the best way to represent it

Without an example it is dfficult to be specific, but you should investigate the colorscheme["valuesplit"....] plot option

Consider the following "toy" example

#
# Generate points for an example function
#
    pts:=[seq([k, sin(k)], k=-evalf(Pi)..evalf(Pi),0.25)]:
#
# Get the y-values of the plots and sort them in
# ascending order
#
   lims:= sort([seq( j[2], j in pts)]):
#
# valuesplit the colorscheme so that the maximum
# value (given by lims[-1]) is one color, and all
# others (range is given by lims[1]..lims[-2] )
# are a differentt color
#
   plots:-pointplot( pts,
                           colorscheme=[ "valuesplit",
                                                [ lims[1]..lims[-2]="Red",
                                                  lims[-2]..lims[-1]="Blue"
                                                ]
                                              ],
                           symbolsize=20
                       );



a:=[[x+y+z],[2*x+y+z],[x-y+z],[2*x+y+z]];
[{seq([ListTools:-SearchAll(j,a)], j in a)}[]];

seems to work pretty well, with the proviso that you are not too bothered aboout the order of the returned entries - which I am pretty sure I could fix with a sort command

  1. Can confirm that I get pretty much the same message as you when I try to multithread, os maybe I have to keep looking at this
  2. On the other hand if I make a small adjustemnt to your code, I (probably) wouldn't bother to multithread - the overhead of setup/result collation etc would (probably) outweigh actual calculation time in a single thread!

Modified code is as follows

  restart;
  with(LinearAlgebra):
  with(Student[MultivariateCalculus]):
  with(plots):
  with(Statistics):
  with(Threads):
  with(RealDomain): #I know, I don't need all of that here
   t1:=time(); #start timer
   cdf := x-> int(PDF(GammaDistribution(4.5/100,2.5),xi-0.068),xi=0.068..x);
   N:= 20:
#
# Change solve() to fsolve() - after all OP deo want a simple numeric
# answer
#
   sample := (i) -> fsolve( (1/(2*(N+1))+1/(N+1)*i)=cdf(x),x); #this is an equidistant sampling
#smpl := [Seq(sample(i),i=1..(N-1))]; #I try to run a multithread seq -> got it from the documentation
   smpl:=seq([i, sample(i)],i=1..N-1);
   t2:=time()-t1;

Note the use of fsolve() rather than solve() in the definition of the sample() function.

I also include the index value in the output of the above smpl() command, because when I tried the single-threaded approach with the your original sample/solve() command, some index values did return a value at all - I wanted to check this issue. But I did not resolve this. All I establishe was that for some index values solve() returned no solution

Using fsolve(), the included timer commands indicate that the single-threaded version provides the required answers in 3.3secs (on my machine). Doubt if running multi-threaded will do much better - setup and return overheads versus actual calculation time etc, etc

Of course your original calculation should run multithreaded. At the moment I can't figure out why it doesn't, but my Maple installation seems to be partially broken at the moment ( :-( ) ) so I have to get ot the bootom of that first

 

p1 := pds:-plot(x = .16, t = 0 .. 10, numpoints = 10000, color = red);
p2 := pds:-plot(x = .21, t = 0 .. 10, numpoints = 10000, color = blue);
plots:-display([p1, p2]);

Obviusly you may want a different range for t, or more/fewer point for "smoothness"

It seems that you have to define the foci option using 'points' rather than a simple list of values, as in

restart;
with(geometry):
point(A, 0,1);
point(B, 4,1);
ellipse(e1,['foci'=[A,B], 'MajorAxis' = 8],[x,y]);
detail(e1);

Manual does say (my emphasis)

foi is a list of two points which are the foci of the ellipse

 

  1. It is perfectly possible for a procedure to update 'global' values - although whether you regard this as 'acceptable' programmng practice will get you into deep "philosophical" discussion. There are people who regard the use of 'global' parameters as completely unacceptable - makes code difficult to 'comparementalise', 're-use' and debug, etc, etc, etc. Personally I prefer to avold the use of global variables, but I don't get too "fundamentatlist" over it. I do recommend using a 'global' statement within a procedure, whenever global values in the calling code are necessary. This is never strictly required but serves as a reminder, that the correct execution of the procedure may depend on the correct definition of variables outwith the procedure
  2. When returning values from a procedure, the default is to return the output of the last statement executed. However if you have multiple conditional statements (or other similar structures) it can be difficult to determine where the return points actually are (especially a week or so after you have written the code!). For this reason, I tend to use the explicit 'return' command. So far as I am aware this is never absolutely necessary. On the other hand if I look at a procedure with half a dozen 'return' statements then I know that there are exactly half a dozen way for this procedure to exit. This aids debug/reuse (at least for me)
  3. The argument(s) for a return statement can be pretty much anything - eg a single value (which may be a list of multiple values) or a sequence of multiple values. In the case of the latter, the only thing you have to ensure is that the calling statement contains enough arguments for the outputs to be assigned. For example consider the "toy" example

  powers:= proc(x)
                       return x, x^2, x^3:
                end proc;

          along with the calling sequence

  a,b,c:=powers(2);

          Then 'a' will be assigned to 2, 'b' to 4 and 'c' to 8.

          An alternative would be to write the procedure to return a single list, as in

          powers:= proc(x)
                               return [x, x^2, x^3]:
                        end proc;

          With the calling sequence

          a:=powers(2);

          In which case 'a' is the list [2, 4, 8] and you can subsequently access values as a[1], a[2] etc.

Determining the 'best' way to return multiple values from a procedure is never an 'absolute' - it depends on subjective criteria such as readability, debuggability, reusability etc

I think it is easier to see what is goiing on if you plot, B[1](t), B[2](t), and C(T) separately as functions of t. You can see that B[1]() asymtotes to about 0.29, B2() asymptotes to about 0.17 and C(T) asymptotes to about 100. This makes the output of your DEPlot3D command easier (for me at least to interpret).

In one of your equations you had B[1], rather than B[1](t) - I have corrected this. Check out the following

restart;
Model := [diff(B[1](t), t) = k[a1]*C(t)*(R-B[1](t)-B[2](t))-k[d1]*B[1](t),
          diff(B[2](t), t) = k[a2]*C(t)*(R-B[1](t)-B[2](t))-k[d2]*B[2](t),
          diff(C(t), t) = (-(k[a1]+k[a2])*C(t)*(R-B[1](t)-B[2](t))+k[d1]*B[1](t)+k[d2]*B[2](t)+k[m]*((I)(t)-C(t)))/h
         ];
DissMod := subs((I)(t) = 0, Model);
AssMod := subs((I)(t) = C[T], Model);


Pars := [k[a1] = 6*10^(-4), k[d1] = 7*10^(-3), k[a2] = 5*10^(-4), k[d2] = 10^(-2), R = .5, k[m] = 10^(-4), C[T] = 100, h = 10^(-6)];

DETools:-DEplot3d(subs(Pars, AssMod), [B[1](t), B[2](t), C(t)], t = 0 .. 1000, number = 3, B[1] = 0 .. .5, B[2] = 0 .. .5, [[B[1](0) = 0, B[2](0) = 0, C(0) = 0]], scene = [B[1](t), B[2](t), C(t)], maxstep = .1, maxfun = 0);
desys:=[subs(Pars, AssMod)[],B[1](0) = 0, B[2](0) = 0, C(0) = 0];
sols:=dsolve( desys, numeric, maxfun=0)
plots:-odeplot(sols,[t, B[1](t)], t=1..1000);
plots:-odeplot(sols,[t, B[2](t)], t=1..1000);
plots:-odeplot(sols,[t, C(t)], t=1..1000);

Try this

restart:
Pars := [k[a1] = 6*10^(-4), k[d1] = 7*10^(-3), k[a2] = 5*10^(-4), k[d2] = 10^(-2), R = .5, k[m] = 10^(-4), C[T] = 100, h = 10^(-6)];
eqns:=[C = -(k[d2]*B[2]+C[T]*k[m]+k[d1]*B[1])/((B[1]+B[2]-R)*k[a1]+(B[1]+B[2]-R)*k[a2]-k[m]),

       C = k[d1]*B[1]/(k[a1]*(R-B[1]-B[2]))
      ];
eqns2:=eval(eqns, Pars);
p1:=plots:-implicitplot3d(eqns2[1], B[1]=0..0.5, B[2]=0..0.5, C=0..100, style=patchnogrid, color=red, grid=[50,50,50], transparency=0.5):
p2:=plots:-implicitplot3d(eqns2[2], B[1]=0..0.5, B[2]=0..0.5, C=0..100, style=patchnogrid, color=blue, grid=[50,50,50], transparency=0.5):
p3:=plots:-intersectplot(eqns2[1],eqns2[2], B[1]=0..0.5, B[2]=0..0.5, C=0..100, thickness=6,color=black):
plots:-display([p1,p2,p3])

The following code provides a partial answer

restart:
Pars := [k[a1] = 6*10^(-4), k[d1] = 7*10^(-3), k[a2] = 5*10^(-4), k[d2] = 10^(-2), R = .5, k[m] = 10^(-4), C[T] = 100, h = 10^(-6)];
eqns:=[C = -(k[d2]*B[2]+I*k[m]+k[d1]*B[1])/((B[1]+B[2]-R)*k[a1]+(B[1]+B[2]-R)*k[a2]-k[m]),

       C = k[d1]*B[1]/(k[a1]*(R-B[1]-B[2]))
      ];
eqns2:=eval(eqns, Pars);
plots:-implicitplot3d(eqns2[1], B[1]=0..0.5, B[2]=0..0.5, C=0..100);
plots:-implicitplot3d(eqns2[2], B[1]=0..0.5, B[2]=0..0.5, C=0..100);

The reason that nothing is plotted for eqns2[1] is that your first equation contains 'I' which will be interpreted as the square root of -1. Is that intentional?

then try

T:=1..18;
TB:={1,2,3,4,7,8,15,18};
DD:={5,7,11,12,14,15,17};
e:=1..9;
K:=1..17;
add( add( add( `if`( j=i, 0, x[i,j,k]+x[j,i,k] ), i in TB), j in TB), k in DD);
add( add( add( `if`( i=2*j-1, 0, x[2*j-1,i,k] ), i in T), j in e), k in K);

First 98 99 100 101 102 103 104 Last Page 100 of 123