tomleslie

Some suggestions...

Answered: tomleslie 13876

March 26 2019

2 2

Maple 7 was superseded seventeen years ago, so I doubt if anyone here still has a copy: I certainly don't! I only have Maple releases going back five years, so I can't test anything for your release - just make some suggestions which you can try

If the [$65..95] construct is casing a problem - don't use it. The construct [seq(i, i=65..95)] *ought* to do the same thiing
Although the convert(...., `bytes`) command appears to work, I notice from the AMple update history that the StringTools package was introduced in Maple 7 and this contains specific commands for converting integers to letters and vice versa. Thus the construct [seq(StringTools[Char](i), i=65..95)] *ought to do exactly the same as convert( [$65..95], 'bytes' ). And possibly without the "Error" massage?!
The warnng messages about the redefnition of "changecoords" and "arrow" are not cause by anything in this worksheet. I can only assume that these arise because of code you have in your Maple initialisation file (maple.ini). If you want to get rid of these warning messages you will need to locate/edit this initialisation file. Depending on your installation, the relevant file could be in a variety of locations. Try typing "?worksheet,reference,installation" or "?maple.ini" (without the quotes) at the Maple prompt and you ought to get some idea where your maple.ini file is located

Parallelization...

Answered: tomleslie 13876

March 25 2019

0 8

If I take the fastest code I had in my earlier answer and increase the size of the problem by 10X to nelems=10000000, then this sequential code takes about 55secs on my machine (compared with ~6secs for nelems=1000000, so execution time looks pretty linear with size

If I parallelize the code into 4 Tasks, then the same size problem (ie nelems=10000000), then the execution time drops to ~23secs on my machine. (Parallelized code is a bit rough, I was in a hurry - sorry)

So a further 2.5X speed-up

Check the attached - be aware that total execution times is about 85 secs (on my machine)

>

#
# Using the fastest sequential code I cam up with
# with nelems=10000000. Takes about 55 secs on my
# machine
#
  restart;
  nelems := 10000000:
  n := 374894756873546859847556:
  op(n):
  A := Array(1 .. 4, 1 .. nelems):
  length(n):
  st := time[real]():
  A[1,1..nelems]:= Array([seq( i^10*n, i=1..nelems)]):
  A[2,1..nelems]:= Array([seq( length(A[1, i])-1, i=1..nelems)]):
  A[3,1..nelems]:= Array([seq(iquo(A[1,i], 10^(A[2, i]-2)),i=1..nelems)]):
  A[4,1..nelems]:= Array([seq(irem(A[1,i],1000),i=1..nelems)]):
  time[real]() - st;
  A:
  A[1, -2];
  A[2, -2];
  A[3, -2];
  A[4, -2];

(1)

>

#
# Parallelize the above across four tasks, Takes about
# 23 secs on my machine, so a speed up of 2.5X on the
# above sequential code
#
  restart;
  nelems := 10000000:
  n := 374894756873546859847556:
  op(n):
  length(n):
  st := time():
  doTask:=proc(sv, nv)
               local rng:=nv-sv+1,
                     B:=Array(1..4,1..nv-sv+1),
                     i :
               B[1,1..rng]:= Array([seq( (sv-1+i)^10*n, i=1..rng)]):
               B[2,1..rng]:= Array([seq( length(B[1, i])-1, i=1..rng)]):
               B[3,1..rng]:= Array([seq(iquo(B[1,i], 10^(B[2, i]-2)),i=1..rng)]):
               B[4,1..rng]:= Array([seq(irem(B[1,i],1000),i=1..nv-sv+1)]):
               return B;
         end proc:
  time() - st;
  contTask:=proc( V,W,X,Y )
                  return copy(ArrayTools:-Concatenate(2, V, W, X, Y)):
           end proc:

  setupTask:=proc(n)
                  Threads:-Task:-Continue
                  ( contTask,
                    Task=[doTask, 1, n/4 ],
                    Task=[doTask, n/4+1, n/2],
                    Task=[doTask, n/2+1, 3*n/4],
                    Task=[doTask, 3*n/4+1, n]
                  );
            end proc:

  st := time[real]():
  A:=Threads:-Task:-Start( setupTask, nelems):
  time[real]() - st;
  A[1, -2];
  A[2, -2];
  A[3, -2];
  A[4, -2];

>

(2)

>	interface(rtablesize=10);

(3)

>

Download speedUP2.mw

I think I'd call thi a side-effect...

Answered: tomleslie 13876

March 25 2019

0 0

I'm pretty sure that any Matrix initialiser function which takes two arguments can be used.

For example, the following

#
# Consider an initialiser function
# which takes two arguments
#
f:=(x,y)->x+y;
Matrix(3, 3, f);

will "work".

When you use eg `+`, this is just the "prefix form" of the + operator, and will accept two (or more) arguments. So for example

`+`(1,2,3,4,5,6);

will correctly reurn 21

Can confirm...

Answered: tomleslie 13876

March 25 2019

2 1

that this is a version issue, because I see the same problem in my Maple 2018.2 installation - see the attached.

It seems to be a "bug" in odetest(), which got "fixed" with the release of Maple 2019

Using Maple 2018.2, one can check validity without using odetest(). Just substitute the proposed solutions into the original ODEs/boundary conditions. This confirms that the solutions are OK, so the problem must be with odetest()..

See the attached.

>	# # Initialize and show version: possibly a version # Maple issue on this problem # restart: kernelopts(version); Physics:-Version();

$"C:\Users\TomLeslie\maple\toolbox\2018\Physics Updates\lib\Physics Updates.maple", `2019, March 14, 23:22 hours, version in the MapleCloud: 341, version installed in this computer: 326.`$

(1)

>

#
# Using lists rather than sets tests two solutions
# and four BCs, so odetest() should return 6 zeros
# if everything is good - but it doesn't.
#
# Last two entries in the output list are non-zero,
# suggesting that the problem arises with the final
# two boundary conditions
#
  laplace := -phiL(z) + diff(phiL(z),z$2)=0,
             -phiM(z) + diff(phiM(z),z$2)=0:

  BCs:= phiL(d1)=0,
         phiM(-d1)=0,
         phiL(0)=phiM(0),
         D(phiL)(0)-D(phiM)(0)=-n:

  sol:= [ phiM(z) = n/2/coth(d1)*(cosh(z)+coth(d1)*sinh(z)),
          phiL(z) = n/2/coth(d1)*(cosh(z)-coth(d1)*sinh(z))
        ]:

  odetest(sol, [laplace, BCs]);

(2)

>

#################################################
# Test the ode system with the obtained solutions.
#
# Should give 0=0, 0=0, which it does!
#
  eval( laplace, sol);
#################################################
# Test the first BC, phiL(d1)=0 with the relevant
# solution.
#
# Should result in 0, which it does!
#
  simplify
  ( eval
    ( eval
      ( phiL(z), sol),
      z=d1
    ),
    trig
  );
##################################################
# Test the second BC phiM(-d1)=0 with the relevant
# solution.
#
# Should result in 0, which it does!
#
  simplify
  ( eval
    ( eval
      ( phiM(z), sol),
      z=-d1
    ),
    trig
  );
##################################################
# Test the third BC phiL(0)=phiM(0) with the
# relevant solutions.
#
# Should result in 0, which it does!
#
  simplify
  ( eval
    ( eval
      ( phiL(z), sol),
      z=0
    ),
    trig
  )-
  simplify
  ( eval
    ( eval
      ( phiM(z), sol),
      z=0
    ),
    trig
  );
###############################################
# Test the fourth BC D(phiL)(0)-D(phiM)(0)=-n
# with the relevant solutions.
#
# Should result in -n, which it does!
#
  eval
  ( diff
    ( eval
      ( phiL(z), sol),
      z
    ),
    z=0
  )-
  eval
  ( diff
    ( eval
      ( phiM(z), sol),
      z
    ),
    z=0
  );

(3)

>

Download odeTest2.mw

Two issues...

Answered: tomleslie 13876

March 24 2019

1 3

Simple problem - you use sets rather than lists for both the ode system and the solution functions. Sets do not maintain the 'order' in which you enter elements. In fact elements will be sorted (more-or-less) lexicographically. This will cause arguments to the odetest() command to be p[resented in a slightly different order - so the output may(?) be in a (correspondingly? different order.Since you have two solutions and four boundary conditions, there ought to be 6 outputs (ideally all zero). If you use sets as arguments to odetest() then these six outputs will (ideally) collapse to a set with the single element 0. since sets cannot contain duplicate elements. Since you are getting three outputs, this implies that you are getting four zeros and two non-zeros: the set container collapses the four zeros to a single zero, so you get an output set of three eleemnts, I recommend that you use list arguments rather than sets, as in the attached: this will maintain all six (hopefully) zeros and if anything is non-zero you can work out which solution/bc "failed"
I'm restricted in identifying which two solutions/BCs are causing issues, because as the attached shows, using either sets or lists, I get either {0} or [0,0,0,0,0,0], indicating that odetest() thinks that everything is fine. This may be a version issue. What version are you running - check by using the same commands I used at the top of the attached.

For some reason this worksheet won't display here - no idea why: you will have to download to view

Download odeTest.mw

Simple(?) 5x speed-up...

Answered: tomleslie 13876

March 24 2019

1 1

The attached performs your calculation three different ways: your original code, then two "improved" versions for nelems=1000000, which is long enough to make timings more-or-less reliable, but short enough that I don't get bored waiting. Before considering parallelizing I think it is important to make the "simple" code as fast as possible

The overall speed-up in the attached is about 5X. I may get a chance to look at improvements by parallelization later

>

#
# OP's original - takes ~35secs for nelems=1000000
# on my machine
#
  restart;
  nelems := 1000000:
  n := 374894756873546859847556:
  op(n):
  A := Array(1 .. 4, 1 .. nelems):
  length(n):
  st := time():
  for i to nelems do
      A[1, i] := i^10*n;
      A[2, i] := length(A[1, i]) - 1;
      b := convert(A[1, i], string);
      A[3, i] := parse(b[1 .. 3]);
      A[4, i] := parse(b[-3 .. -1]);
  end do:
  time() - st;
  A:
  A[1, -2];
  A[2, -2];
  A[3, -2];
  A[4, -2];

(1)

>

#
# Simple code rewrite. Gives about 4X speed-up.
# Takes about 9.5secs for nelems=1000000 on my
# machine
#
  restart;
  nelems := 1000000: # ~ takes about 50mins when n=100,000,000
  n := 374894756873546859847556:
  op(n):
  A := Array(1 .. 4, 1 .. nelems):
  length(n):
  st := time():
  for i to nelems do
      A[1, i]:= i^10*n;
      A[2, i]:= length(A[1, i]) - 1;
      A[3, i]:= iquo(A[1,i], 10^(A[2, i]-2));
      A[4, i]:= irem(A[1,i],1000);
  end do:
  time() - st;
  A:
  A[1, -2];
  A[2, -2];
  A[3, -2];
  A[4, -2];

(2)

>

#
# Use seq() commands rather than for-loop. Gives
# about 30% further speed-up. Takes about 6.5secs
# for nelems=1000000 on my machine
#
  restart;
  nelems := 1000000:
  n := 374894756873546859847556:
  op(n):
  A := Array(1 .. 4, 1 .. nelems):
  length(n):
  st := time():
  A[1,1..nelems]:= Array([seq( i^10*n, i=1..nelems)]):
  A[2,1..nelems]:= Array([seq( length(A[1, i])-1, i=1..nelems)]):
  A[3,1..nelems]:= Array([seq(iquo(A[1,i], 10^(A[2, i]-2)),i=1..nelems)]):
  A[4,1..nelems]:= Array([seq(irem(A[1,i],1000),i=1..nelems)]):
  time() - st;
  A:
  A[1, -2];
  A[2, -2];
  A[3, -2];
  A[4, -2];

(3)

>	interface(rtablesize=10);

(4)

>

Download speedUP.mw

In your original post...

Answered: tomleslie 13876

March 23 2019

0 0

that is before you edited/removed the pictures from your textbook, the first example had the initial conditions tau=0, prey(0)=40, pred(0)=16, although variable names 'pred' and 'prey' were given as 'x' and 'y'.

The attached reproduces the first few graphs from your textbook piccies with these conditions.

Observation

In a predator/prey model, the numbers of both quantities should be integers, and you should be using difference equations rather than differential equations: non-integer numbers of either predator or prey do not make sense - unless you can define something like (say) 0.1 of a lion!!

I supppose that if the numbers of both predator and prey are "sufficiently large", then using continuou variables may be an acceptable(?) approximation.

However in your revised post you start with pred(0)=1 and prey(0)=1. This fails the "sufficiently large numbers" test mentioned above. It is also (logically) nonsense because

with prey(0)=1, no reproduction is possible. Even without a predator, the number of prey would stay at 1.Furthermore, given a moderately peckish predator, the number of prey will decline from 1 to 0 quite quickly. Some time later the predator will also die of starvation!

Attached worksheet won't display inline here: it contains the "rtablesize" workaround so the problem is probably one of the plot options I have used ( 'size=' maybe?). You will have to download to view :-(

Download predPrey.mw

Yet another...

Answered: tomleslie 13876

March 22 2019

2 0

slightly different(?) way

>

  restart;
  with(plots):
  nFrames:=100:
  step:=evalf(2*Pi/nFrames):
  display( [ pointplot
             ( [ seq( [j*Pi, 0], j=0..6 )],
               symbol=solidcircle,
               symbolsize=10,
               color=red
             ),
             display
             ( [ seq( plot( [ sin(x+phi),
                              sin(x-phi),
                              sin(x+phi)+sin(x-phi)
                            ],
                            x=0..6*Pi,
                            color=[blue, green, black]
                          ),
                      phi=0..step*(nFrames-1), step
                    )
               ],
               insequence=true
             )
           ],
           axes=boxed,
           gridlines=false
         );

>	interface(rtablesize=10);

(1)

>

Download stWave.mw

Observation...

Answered: tomleslie 13876

March 21 2019

0 1

Your "desired" result contains the variable name G_r but this variable does not occur in the problem statement. It is therefore impossible to obtain the former from the latter!

It is however relatively(?) simple to get the term containing B_r - see the attached. And if G_r were defined anywhere it might be possible to do somehting similar

>

>	restart: with(DETools): a1:=m^2((1-N)/(2-N))(r/2)Dp: a2:=((1-N)m^2/(2-N))(((Nb-Nt)/64)(r^5/6-h^2r/2)-(B[r]/4)(r^3/4-h^2r/2)(Nt/Nb)): a3:=a1-a2: a4:=r^2a3: DE1:=r^2diff(v(r),r,r)+rdiff(v(r),r)-(m^2r^2+1)*v(r)=a4: b1:=dsolve(DE1,v(r));

v(r) = BesselI(1, m*r)*_C2+BesselK(1, m*r)*_C1-(1/128)*(((-(1/3)*r^4+h^2)*Nb^2+(((1/3)*r^4-h^2)*Nt+64*Dp)*Nb-16*Nt*(h^2-(1/2)*r^2)*B[r])*m^4+(-8*Nb^2*r^2+8*Nb*Nt*r^2+64*Nt*B[r])*m^2-64*Nb*(Nb-Nt))*r*(-1+N)/(Nb*m^4*(-2+N))

(1)

>	# # Get the term in B[r] # simplify(collect(coeff(rhs(collect(expand(b1), [B[r], Nt])), B[r]),r))*B[r];

(1/16)*(2*h^2*m^2-m^2*r^2-8)*r*Nt*(-1+N)*B[r]/(Nb*m^2*(-2+N))

(2)

>	interface(rtablesize=10):

>

Download sinpProb.mw

Noy completely...

Answered: tomleslie 13876

March 20 2019

0 0

clear (to me) exactly what you want to be included in the animation. Possibly the one shown in the attachment

>

with(plots):
local O:
R:= 5:
alpha:= (1/9)*Pi:
beta:= (1/3)*Pi:
n:= 100:
dt:= 2*Pi/n:
C1:= plot([R*cos(t), R*sin(t), t = 0 .. 2*Pi], color = blue):
O:= [0, 0]:
M:= [R*cos(beta), R*sin(beta)]:
A:= [R*cos(alpha), R*sin(alpha)]:
B:= [R*cos(alpha+Pi), R*sin(alpha+Pi)]:
AB:= plot([A, B]):
OM:= plot([O, M]):
P:= [R*cos(t0*dt), R*sin(t0*dt)]:
Q:= [R*cos(dt*t0+Pi), R*sin(dt*t0+Pi)]:
diamob:= seq
           ( display
             ( [ plot( [P, Q], color = green),
                 plot([M, P]),
                 plot([M, Q])
               ]
             ),
             t0 = 0 .. n
           ):
tp:= textplot( [ [A[1]+.3, A[2], "A"],
                   [B[1]-.3, B[2], "B"],
                   [M[1]+.3, M[2]+.3, "M"]
                 ]
               ):
display( [ AB, C1, OM, tp,
            display
            ( diamob,
              insequence = true
            )
          ],
          scaling = constrained
        );

>	interface(rtablesize=10):

>

Download anim.mw

Assuming your circle is 2D...

Answered: tomleslie 13876

March 18 2019

0 0

then the attached is one alternative. (it would be relatively trivial to change for 3D, or handle 2D and 3D, but I don't feel inclined to cover all the bases unless it is really necessary)

Ignore the interface(rtablesize=10) command at the end: this is just a workaround to ensure that worksheets created in Mapale 2019 display properly on this site

>

  restart;
  getCirc:= proc( k, l, m)
                  uses geometry, plots;
                  local ci;
                #
                # Check collinearity
                #
                  point(A, k[1], k[2]);
                  point(B, l[1], l[2]);
                  point(C, m[1], m[2]);
                  if   AreCollinear(A, B, C)
                  then error "Points are collinear"
                  else circle
                       ( c,
                         [A, B, C],
                         [x,y]
                       );
                       return Equation(c),
                              [ pointplot
                                ( [k, l, m],
                                  symbol=solidcircle,
                                  symbolsize=20,
                                  color=blue
                                ),
                                draw
                                ( c,
                                  axes=none,
                                  color=red
                                )
                              ]
                  fi;
            end proc:
#
# If we hae points which are not collinear
#
  expr, plt:= getCirc( [2.5, 1], [0, 5], [-1,2] ):
  plots:-display( plt,
                  size=[1000, 500],
                  title=typeset( "A plot of the circle\n",
                                  expr,
                                 "\nwith defining points"
                               ),
                  titlefont=[times, bold, 20]
                );
#
# Or if we have points which are collinear
#
  expr, plt:= getCirc( [0, 0], [2, 0], [-3, 0]);

Error, (in getCirc) Points are collinear

>	interface(rtablesize=10):

>

Download drawCirc.mw

Use setVertexPositions()...

Answered: tomleslie 13876

March 18 2019

2 1

If I undertsand correctly, you can set the positions of the vertexes in your graph with the SetVertexPositions() command,

See the attached - ignore the interface(rablesize=10) command at teh end of the attached - this is just needed to make things display prooperly on this site

>

  restart;
  with(GraphTheory):
#
# Define adjacency Matrix
#
  A:= Matrix(9, 9, [ [0, 0, 1, 0, 1, 1, 1, 0, 0],
                     [1, 0, 0, 0, 0, 0, 0, 0, 0],
                     [0, 1, 0, 0, 0, 0, 0, 0, 0],
                     [1, 0, 0, 0, 0, 0, 0, 0, 0],
                     [1, 0, 0, 0, 0, 0, 0, 0, 0],
                     [1, 0, 0, 0, 0, 0, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0, 1],
                     [0, 0, 0, 1, 0, 0, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 1, 0]
                   ]
          ):
#
# Define graph based on adjacency matrix and then
# draw it
#
  G1:=Graph(A):
  DrawGraph( G1 );
#
# Define and set vertext positions for the graph
# then draw it
#
  vp:=[ [40, 40], [22, 22], [36, 26], [21, 45], [45, 35],
        [55, 20], [55, 45], [26, 59], [55, 65]
      ]:
  SetVertexPositions(G1, vp):
  DrawGraph(G1);

>	interface(rtablesize=10):

>

Download drGraph.mw

With minimal comment...

Answered: tomleslie 13876

March 17 2019

0 0

@vv

The IntegrationTools:-Parts() command will provide the "desired"(??) solution - which looks "suspicious" to me in the absence of limits

Download parts.mw

Caught out by using 2D input - again!...

Answered: tomleslie 13876

March 15 2019

0 0

Your basic problem is that (in certain execution groups) you have attempted to use 2D-input with manual code indentation (presumably for readability). Just don't!!! The 2-D parser will place all sorts of "interesting" interpretations on the extra "whitespace" which arises from this indentation. Obviously this never happens if you use 1-D input because, because in this format "whitespace" is just "whitespace"

If you must use 2D-input, you *might* get away with using manual indentation provided you put commands and their associated opening bracket on the same line. So write

seq(

but never

seq
(

Apply the same principle for all other commands - eg Transpose etc. Randomly typing "shift C/R" and "whitespace" (for indentation purposes) whilst using 2D-input mode is just asking for trouble!

In the attached, I have simply removed the "indentation" from your code, and everything now works

>

>	# # Load packages "properly" and check Maple # version # with(DirectSearch): with(LinearAlgebra): kernelopts(version)

(1)

>

(2)

>

(3)

>

(4)

>

(5)

>

(6)

>

$for n to 20 do e1 := t[0] = (-p[0, fail]^7+1)/((1-p[0, fail])*((-p[0, fail]^7+1)/(1-p[0, fail])+(17/2+(17/2)*p[0, fail]+(33/2)*p[0, fail]^2+(33/2)*p[0, fail]^3+(65/2)*p[0, fail]^4+(65/2)*p[0, fail]^5+(65/2)*p[0, fail]^6)/p[0, idle]+(1-r[0])*(p[0, fail]^6+(-p[0, fail]^6+1)*p[0, succ]/(1-p[0, fail]))/r[0])); e2 := t[1] = (-p[1, fail]^7+1)/((1-p[1, fail])*((-p[1, fail]^7+1)/(1-p[1, fail])+(9/2+(9/2)*p[1, fail]+(17/2)*p[1, fail]^2+(17/2)*p[1, fail]^3+(33/2)*p[1, fail]^4+(33/2)*p[1, fail]^5+(33/2)*p[1, fail]^6)/p[1, idle]+(1-r[1])*(p[1, fail]^6+(-p[1, fail]^6+1)*p[1, succ]/(1-p[1, fail]))/r[1])); e3 := t[2] = (-p[2, fail]^7+1)/((1-p[2, fail])*((-p[2, fail]^7+1)/(1-p[2, fail])+(5/2+(5/2)*p[2, fail]+(9/2)*p[2, fail]^2+(9/2)*p[2, fail]^3+(17/2)*p[2, fail]^4+(17/2)*p[2, fail]^5+(17/2)*p[2, fail]^6)/p[2, idle]+(1-r[2])*(p[2, fail]^6+(-p[2, fail]^6+1)*p[2, succ]/(1-p[2, fail]))/r[2])); e4 := t[3] = (-p[3, fail]^7+1)/((1-p[3, fail])*((-p[3, fail]^7+1)/(1-p[3, fail])+(3/2+(3/2)*p[3, fail]+(5/2)*p[3, fail]^2+(5/2)*p[3, fail]^3+(9/2)*p[3, fail]^4+(9/2)*p[3, fail]^5+(9/2)*p[3, fail]^6)/p[3, idle]+(1-r[3])*(p[3, fail]^6+(-p[3, fail]^6+1)*p[3, succ]/(1-p[3, fail]))/r[3])); e5 := t[4] = (-p[4, fail]^7+1)/((1-p[4, fail])*((-p[4, fail]^7+1)/(1-p[4, fail])+(1+p[4, fail]+(3/2)*p[4, fail]^2+(3/2)*p[4, fail]^3+(5/2)*p[4, fail]^4+(5/2)*p[4, fail]^5+(5/2)*p[4, fail]^6)/p[4, idle]+(1-r[4])*(p[4, fail]^6+(-p[4, fail]^6+1)*p[4, succ]/(1-p[4, fail]))/r[4])); e6 := p[0, fail] = n*t[0]*(1-(1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n)+0.1590000000e-1*n*t[0]*(1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e7 := p[1, fail] = n*t[1]*(1-(1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n)+0.1590000000e-1*n*t[1]*(1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e8 := p[2, fail] = n*t[2]*(1-(1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n*(1-t[3])^n*(1-t[4])^n)+0.1590000000e-1*n*t[2]*(1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e9 := p[3, fail] = (1/2)*n*t[3]*(1-(1-t[3])^(n-1)*(1-t[4])^n)+(1/2)*n*t[3]*(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n)+(0.1590000000e-1*((1/2)*n*t[3]*(1-t[3])^(n-1)*(1-t[4])^n+(1/2)*n*t[3]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e10 := p[4, fail] = (1/2)*n*t[4]*(1-(1-t[4])^(n-1)*(1-t[3])^n)+(1/2)*n*t[4]*(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n)+(0.1590000000e-1*((1/2)*n*t[4]*(1-t[4])^(n-1)*(1-t[3])^n+(1/2)*n*t[4]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e11 := p[0, succ] = .9841000000*n*t[0]*(1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e12 := p[1, succ] = .9841000000*n*t[1]*(1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e13 := p[2, succ] = .9841000000*n*t[2]*(1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e14 := p[3, succ] = (.9841000000*((1/2)*n*t[3]*(1-t[3])^(n-1)*(1-t[4])^n+(1/2)*n*t[3]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e15 := p[4, succ] = (.9841000000*((1/2)*n*t[4]*(1-t[4])^(n-1)*(1-t[3])^n+(1/2)*n*t[4]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e16 := p[0, idle] = (1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n; e17 := p[1, idle] = (1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n; e18 := p[2, idle] = (1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n; e19 := p[3, idle] = (1/2)*(1-t[3])^(n-1)*(1-t[4])^n+(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n; e20 := p[4, idle] = (1/2)*(1-t[4])^(n-1)*(1-t[3])^n+(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n; CSTR := `~`[`=`]({p[0, fail], p[0, idle], p[0, succ], p[1, fail], p[1, idle], p[1, succ], p[2, fail], p[2, idle], p[2, succ], p[3, fail], p[3, idle], p[3, succ], p[4, fail], p[4, idle], p[4, succ], t[0], t[1], t[2], t[3], t[4]}, 0 .. 1); residuals := `~`[lhs-rhs]([e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e20]); ansLS[n] := Optimization:-LSSolve(residuals, op(CSTR), iterationlimit = 10000); ansDS[n] := DirectSearch:-SolveEquations(residuals, CSTR, initialpoint = ansLS[n][2]) end do$

>

Matrix(20, 20, [seq(ansDS[j][3], j = 1 .. 10)]); Matrix(20, 20, [seq(ansLS[j][2], j = 1 .. 20)])

RTABLE(18446744074184966142, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 20, 1 .. 20)

RTABLE(18446744074184966294, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 20, 1 .. 20)

(7)

>

RTABLE(18446744074184966446, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 21, 1 .. 21)

RTABLE(18446744074184966598, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 21, 1 .. 21)

(8)

>

mi := max[index](DSMat-LSMat); (DSMat-LSMat)[mi]

(9)

>

Download dsls.mw

Confirming mmcdara' response...

Answered: tomleslie 13876

March 14 2019

0 0

You should never, ever use an indexed variable (eg c[1], c[2]), and a nonindexed variable (eg 'c') in the same calculation!

It si also not a good idea to use non-commutative multiplication ( ie the '.' symbol) rather than simple multiplication (ie the '*' symbol), unless your are really sure that you need the former

I changed refereence to 'c' to 'cc' and all '.' characters to '*' and came up with the attached. Be aware that the dsolve command takes ~3mins to completet on my machine - yours may be faster (or slower!)

>

restart;
with(plots):
r := 0.5e-2; k := 10000; a := .4; alpha := .25; epsilon := 0.2e-1;
mu := 0.4e-1; delta := 0.3e-2; Lambda := 0.2e-1;
beta[2] := .45; beta[1] := .2; c[1] := 2; c[2] := 5;
w[1] := 10; w[2] := 30; T := 3;
u[1] := min(max(0, z), 1); z := beta[2]*s(t)*i(t)*(lambda[2](t)-lambda[1](t))/(w[1]*(s(t)+i(t)+e(t)));
u[2] := min(max(0, cc), 1);
cc := (beta[1]*s(t)*p(t)*(lambda[2](t)-lambda[1](t)))/(w[2]*(a+p(t)))+(lambda[2](t)*e(t)+(i(t)+alpha.e(t))*lambda[3](t)-(gamma.i(t)+p(t)).lambda[4](t))/w[2];

sys := diff(s(t), t) = r*s(t)*(1-(s(t)+i(t)+e(t))/k)-beta[1]*s(t)*p(t)*(1-u[2])/(a+p(t))-beta[2]*s(t)*i(t)*(1-u[1])/(s(t)+i(t)+e(t)),
       diff(e(t), t) = beta[1]*s(t)*p(t)*(1-u[2])/(a+p(t))+beta[2]*s(t)*i(t)*(1-u[1])/(s(t)+i(t)+e(t))-(mu+alpha+u[2])*e(t),
       diff(i(t), t) = (alpha+u[2])*e(t)-(mu+epsilon+u[2])*i(t),
       diff(p(t), t) = Lambda+(epsilon+u[2])*i(t)-delta*p(t),
       diff(lambda[1](t), t) = -lambda[1](t)*(r*(1-(2*s(t)+i(t)+e(t))/k)-beta[1]*p(t)*(1-u[2])/(a+p(t))-beta[2]*i(t)*(1-u[1])/(s(t)+i(t)+e(t)))-lambda[2](t)*(beta[1]*p(t)*(1-u[2])/(a+p(t))-beta[2]*i(t)*(1-u[1])/(s(t)+i(t)+e(t))),
       diff(lambda[2](t), t) = -c[1]+lambda[1](t)*r*s(t)/k+lambda[2](t)*(mu+alpha+u[2])-(1-u[2])*alpha*lambda[3](t),
       diff(lambda[3](t), t) = -c[2]+lambda[1](t)*(r*s(t)/k+beta[2]*s(t)*(1-u[1])/(s(t)+i(t)+e(t)))-lambda[2](t)*beta[2]*s(t)*(1-u[1])/(s(t)+i(t)+e(t))+lambda[3](t)*(u[2]+mu+gamma)-lambda[4](t)*gamma*(1-u[2]),
       diff(lambda[4](t), t) = ((lambda[1](t)*beta[1])*s(t)*a*(1-u[2]))/(a+p(t))^2-((lambda[2](t)*beta[1])*s(t)*a*(1-u[2]))/(a+p(t))^2-lambda[4](t)*(delta+u[2]),
       s(0) = 1000, e(0) = 10, i(0) = 0, p(0) = 100, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0, lambda[4](T);
p1 := dsolve({sys}, numeric, abserr=1.0e-03, maxmesh=1024);

0.4500000000e-1*s(t)*i(t)*(lambda[2](t)-lambda[1](t))/(s(t)+i(t)+e(t))

0.6666666667e-2*s(t)*p(t)*(lambda[2](t)-lambda[1](t))/(.4+p(t))+(1/30)*lambda[2](t)*e(t)+(1/30)*(i(t)+.25*e(t))*lambda[3](t)-(1/30)*((gamma*i(t)+p(t)).lambda[4](t))

(1)

>

plots:-odeplot(p1, [t, e(t)],t=0..3);
plots:-odeplot(p1, [t, i(t)],t=0..3);
plots:-odeplot(p1, [t, p(t)],t=0..3);
plots:-odeplot(p1, [t, s(t)],t=0..3);
plots:-odeplot(p1, [t, lambda[1](t)],t=0..3);
plots:-odeplot(p1, [t, lambda[2](t)],t=0..3);
plots:-odeplot(p1, [t, lambda[3](t)],t=0..3);
plots:-odeplot(p1, [t, lambda[4](t)],t=0..3);