acer

18951 Reputation

29 Badges

14 years, 275 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

When you load a package using the with command Maple rebinds the global name, such that a newly parsed expression will use whatever is the current binding of the name.

Several people throw around the term "overload" for this rebinding process, but that is a mistake. In Maple the term overload has quite a different technical meaning, and it's incorrect and wrong to mix up the two.

There is a procedure VectorCalculus:-`.` which is one of the package exports of the VectorCalculus package. When the VectorCalculus package is loaded, say by executing with(VectorCalculus) the name `.` is rebound to VectorCalculus:-`.` and that is what then gets called in your second example.

Your second example called with(VectorCalculus); but you omitted to show  the output displayed by that call. It contains a list of all the name that have been rebound to the package exports. One of those names printed is `.`, and by looking at that output and by reading the Help page for the with command you should be able to figure out that there is a VectorCalculus:-`.` export.

You can even utilize that export directly, without the whole package being loaded (as a prefix call rather that a infix call, however):

restart;
VC:=VectorCalculus:
f:=y*sin(x)+y^2:
VC:-`.`( VC:-Nabla, VC:-Nabla(f,[x,y]) );

                    2 - y sin(x)

The act of rebinding the name is done via the (undocumented) bind command, which the with command utilizes. Similarly the unwith command utilizes the (undocumented) unbind command. These work like a stack, such that calling unwith will restore the previous binding. I don't recommend anyone using bind and unbind directly; utilize with and unwith instead.

In this stack model, the last rebinding is what gets used. If you load the plots package then Maple rebinds the name arrow to plots:-arrow. If you subsequently load the plottools package then Maple futher rebinds the name arrow to plottools:-arrow, replacing the previous binding. If you then unload the plottools package then the binding of arrow reverts to the previous binding to plots:-arrow.

The name management functionality offered by with/unwith is there to assist people and make common operations easier. It cannot offer perfect management of the namespace because such a thing is impossible: there are cases where two different people want two different behaviors under otherwise identical circumstances.

So there will always be cases where qualifying commands names (using the colon-minus syntax to distinguish between different exports from different packages, and the global names) is necessary or highly preferable.

The mechanisms for "overloading" a name so that a procedure is able to dispatch according to the type of the arguments are different and distinct from rebinding. But the rebinding mechanism and the overloading mechanisms can be constructively utilized together, and together they make for even more power and flexibility of managing the command namespace.

In current Maple there is an overload option for procedures which can allow for a dynamic dispatch functionality according to argument type, say. In other words the code's branches off, according to the type of the arguments. But even before that option was introduced a procedure could be written to examine its arguments and dispatch dynamically.

Several packages have exports which examine the type of the arguments and, according to what is discovered, either do a computation or return a call to the global instance of the same or similar name. For example, Units:-Standard:-`+` will eventually return an expression in terms of the usual global :-`+` after dealing with units or if no units were detected. But some other commands (whose name happens to match a global version) will return as an unevaluated call to themselves rather than to as a call to the global version. There are no hard and fast rules -- the functionality needs, implementations, strengths and weaknesses vary.

Another related piece of functionality for dynamic dispatch involves the static exports of an object.

There is a great deal that could be written about these topics, and I'm not about to write a book here. I suggest reading the Programming Guide and relevant Help pages, and then asking specific questions about specific examples as wanted.

Hopefully I have not misunderstood your intention.

Once I'd set up cartesian coordinates under VectorCalculus I used command-completion in 2D Input.

I typed Del and pressed the Escape key (for Maple GUI command-completion on my Linux), and picked the entry "Del (divergence)" which shows the Nabla-dor markup.

Then I typed Gradient and pressed the Escape key, and picked gthe entry "PhysicsGriadient" because, well, because it seemed to show the right 2D markup.

Here's the result. I used 2D Input only for the line in question. I also include what the input converted to as 1D Notation.

restart;

with(VectorCalculus):

f := y*sin(x)+y^2:

SetCoordinates(cartesian[x,y]):

VectorCalculus[Nabla].Nabla(f)

2-y*sin(x)

(VectorCalculus:-Nabla) . (Nabla(f));

2-y*sin(x)

``

Download VC.mw

It also worked for me if (in 2D Input mode) I used the "Del" entry from the "Common Symbols" palette, along with either a typed period or the center-dot from the same pallete. But I used just f, having already set the coordinates.

Escape the doublequotes, eg.

-c 'path:=\"/home/eduardo/examples\";'

@janhardo Your procedures are named InflectionF which is capitalized, but you accidentally tried to call them as uncapitalized inflectionF.

You should generally not have your procedures rely on having loaded a package at the top-level by the with command. It's better (IMO) to utilize the uses facility within a procedure body, or to utilzie the fully qualified names for commands. I show one way to do that below.

Your declaration of G as local is enough. Your attempt to declare G[1..6,1] as local is invalid. I corrected that.

I leave it to you to add the interval as an additional parameter for your procedures.

 

 

------------------------------------------------------------------------------------------------

Seems that task 6 is about drawing tangentlines in graph for 6 points (x0..x5) and dividing the interval a..b in 5 intervals
So the procedure input  must be
procedure ( function,interval a..b)
Did not add the interval a..b in the interval in the procedure and  the tw procedures below are not working yet.

================================================================

Example 2.3 code ..

restart;

with(Student[Calculus1]):

 

f := x -> x^sin(x):
df := D(f):
ddf := D(D(f)):
G := Array(1..6, 1..4):
P := Array(1..6, 1..2):
G[1..6, 1] := <Student:-Calculus1:-Roots(ddf(x), x=0..16, numeric)>:
for i from 1 to 6 do
  G[i, 2] := f(G[i, 1]):
  G[i, 3] := df(G[i, 1]):
  G[i, 4] := G[i,2]+G[i,3]*(x-G[i,1]):
  P[i, 1] := plot(G[i,4], x=G[i,1]-1..G[i,1]+1,
                  colour=red, adaptive=false, numpoints=2):
  P[i, 2] := plots:-pointplot([G[i,1],G[i,2]], symbol=solidcircle,
                              symbolsize=15, color=blue):
end do:
plots:-display(plot(f, 0..16, color=black),
               seq(seq(P[i,j], i=1..6), j=1..2),
               size=[500,300]);

 

 

Using this code for making a procedure.
Function input in procedure : as expression or as function ?

Procedurename : InflectionF (function, interval,two points in interval  )

Note: in code above the interval has a fixed value, so a variable intvx=a..b is needed

 

restart;

InflectionF := proc(f)
  local df, ddf, G, P, i, j ;
  uses Student:-Calculus1, plots;   
  df := D(f):
  ddf := D(D(f)):
  G := Array(1..6, 1..4):
  P := Array(1..6, 1..2):
  G[1..6, 1] := <Student:-Calculus1:-Roots(ddf(x), x=0..16, numeric)>:
  for i from 1 to 6 do
    G[i, 2] := f(G[i, 1]):
    G[i, 3] := df(G[i, 1]):
    G[i, 4] := G[i,2]+G[i,3]*(x-G[i,1]):
    P[i, 1] := plot(G[i,4], x=G[i,1]-1..G[i,1]+1,
                    colour=red, adaptive=false, numpoints=2):
    P[i, 2] := pointplot([G[i,1],G[i,2]], symbol=solidcircle,
                         symbolsize=15, color=blue):
  end do:
  display(plot(f, 0..16, color=black),
          seq(seq(P[i,j], i=1..6), j=1..2),
         size=[500,300]);
end proc:

InflectionF(x -> x^sin(x));  

 

 

-----------------------------------------------------------------------------------------------------------

Other solution of task .. procedure ?

restart;

with(Student:-Calculus1):
with(plots):

f := x -> x^sin(x);

proc (x) options operator, arrow; x^sin(x) end proc

df := D(f);

proc (x) options operator, arrow; x^sin(x)*(cos(x)*ln(x)+sin(x)/x) end proc

ddf := D(D(f));

proc (x) options operator, arrow; x^sin(x)*(cos(x)*ln(x)+sin(x)/x)^2+x^sin(x)*(-sin(x)*ln(x)+2*cos(x)/x-sin(x)/x^2) end proc

Inflpts := [ fsolve(ddf(x), x=0..16, maxsols=6) ];

[1.395288666, 2.916095859, 7.258616747, 8.576145755, 13.57205647, 14.75675948]

Q := map(p->[p,f(p)], Inflpts);

[[1.395288666, 1.388167079], [2.916095859, 1.270355627], [7.258616747, 5.161057836], [8.576145755, 5.015577540], [13.57205647, 9.048000408], [14.75675948, 8.947326153]]

# The tangent lines, using point-slope form
T := seq(plot(df(p)*(x-p)+f(p), x=p-1..p+1), p=Inflpts):

display(FunctionPlot(f(x), x=0.0..16.0, sign=[], slope=[],
                     caption="", pointoptions=[symbolsize=1],
                     concavity=[color(cyan,magenta),filled(coral,wheat)]),
        T,
        pointplot(Q, symbolsize=10, symbol=solidcircle,
                  color=blue, legend="inflection points"),
        axis[1]=[tickmarks=Inflpts], size=[800,400]);

 

-----------------------------------------------------------------------

restart;

InflectionF := proc(f)
  local df, ddf, G, P, Inflpts, Q, T;
  uses Student:-Calculus1, plots;
  df := D(f);
  ddf := D(D(f));
  Inflpts := [ fsolve(ddf(x), x=0..16, maxsols=6) ];
  Q := map(p->[p,f(p)], Inflpts);
  T := seq(plot(df(p)*(x-p)+f(p), x=p-1..p+1), p=Inflpts):
  display(FunctionPlot(f(x), x=0.0..16.0, sign=[], slope=[],
                       caption="", pointoptions=[symbolsize=1],
                       concavity=[color(cyan,magenta),filled(coral,wheat)]),
          T,
          pointplot(Q, symbolsize=10, symbol=solidcircle,
                    color=blue, legend="inflection points"),
          axis[1]=[tickmarks=Inflpts], size=[800,400]);
end proc:

InflectionF(x -> x^sin(x));

 

 

Download betounes_ex_set_2_opg_6_ac.mw

I suggest that you also add the name x as a local to your procedure, if you intend to continue to utilize it as a dummy name therein. Otherwise your procedure will not work if the name x has been assigned, say, a numeric value at the higher level.

You can use the maplemint command to get an idea of which names in your procedure should be guarded against colliding with assigned values at the higher level. I include this attachment, for some ideas on making your procedures robust against collision with assigned global names.
  betounes_ex_set_2_opg_6_maplemint.mw

Your procedure foo returns an unevaluated name if it is passed a name and is not passed an actual table. And not just any name. Your procedure foo returns the very same name that it is passed (when passed a name).

Even though the formal parameter of foo is specified to have be of type table, it also accepts any name which evaluates to a table. It allows that because tables have last-name-eval, so that they can be operated on in-place as a usefully behaving mutable structure.

The returned name happens to evaluate to a table, unless you somehow unassign it -- which can be done in a few ways.

In your first example B does not get assigned a table per se, but rather B gets assigned the returned, unevaluated global name A which was passed as argument in the call to foo. Your first example unhelpfully used a colon to suppress the output. I change that to a semicolon, which shows that a call to procedure foo returns the name that was passed. Since that A evaluates to a table then a subsequent, full evaluation of B does also. But notice that a 1-level evaluation of B returns just A the name.

restart;

foo:=proc(A::table)
  A["c"]:=3; #add new field
  return A;
end proc:

A:=table():
A["a"]:=1:
A["b"]:=2:

B:=foo(A);
                      B := A

eval(B);
          TABLE(["a" = 1, "c" = 3, "b" = 2])

eval(B,1);
                         A

In your second example you are assigning the unevaluated name A (both passed in the call to foo and also returned by that call to foo) to A itself. You are not explicitly (re)assigning any table to A, but rather you are assigning a name that would evaluate to a table (unless it were reassigned something else -- which is precicely what's happening here!). Recall, foo returns a name, not an actual table.

In the second example your call,

   A:=foo(A);

is assigning the unevaluated name A to A itself. You have probably seen that kind of thing before. It is a common way to unassign a name:

   A := 'A';

In other words, your second example happens to assign to A its very own name, unevaluated. That unassigns it.

That should be enough to explain your second example. It assigns the unevaluated name A to itself. And so A no longer evaluates to anything else.

Of course, if you passed foo an actual table (for example, the value to which A evaluated after its original assignment) then it would also return an actual table.

A:=table():
A["a"]:=1:
A["b"]:=2:

eval(A);
                   TABLE(["a" = 1, "b" = 2])

A:=foo(eval(A)):
eval(A);

               TABLE(["a" = 1, "c" = 3, "b" = 2])

As I mentioned at the start, your procedure foo returns the same name as it gets passed. I think that your examples would be a bit less confusing if it used another name for the formal parameter of the procedure foo. The examples would behave the same if foo were defined as follows:

foo:=proc(H::table)
  H["c"]:=3; #add new field
  return H;
end proc:

Lastly, there are other things with last-name-eval which would also be unassigned by your second kind of example. For instance, procedures also have last-name-eval, so your procedure foo also unassigns f below.

restart;

foo:=proc(K)
  K["c"]:=3; #add new field
  return K;
end proc:

f := x -> sin(x):
eval(x);

x

f := foo(f);

f

eval(f);

f

 

And here are your original two examples

 

A:=table():
A["a"]:=1:
A["b"]:=2:

B:=foo(A):

eval(B,1);

A

eval(B);

table( [( "b" ) = 2, ( "a" ) = 1, ( "c" ) = 3 ] )

A:=table():
A["a"]:=1:
A["b"]:=2:

A:=foo(A):

eval(A);

A

 

Download LNE_example.mw

 

Is this adequate?

restart;

f := proc(Z::[identical(a,b,c,d),identical(a,b,c,d)])
  print("ok");
end proc:

f([a,a]), f([a,b]), f([d,c]);

"ok"

"ok"

"ok"

f([a]);
f([a,b,c]);
f([a,x]);

Error, invalid input: f expects its 1st argument, Z, to be of type [identical(a, b, c, d), identical(a, b, c, d)], but received [a]

Error, invalid input: f expects its 1st argument, Z, to be of type [identical(a, b, c, d), identical(a, b, c, d)], but received [a, b, c]

Error, invalid input: f expects its 1st argument, Z, to be of type [identical(a, b, c, d), identical(a, b, c, d)], but received [a, x]

 

Download type_list.mw

Was the "extra" information to ban repeats there when I first answered this? I don't recall seeing it. But perhaps I didn't notice. If I'd noticed then my answer would include a check for nops of the list. The key thing is the identical(a,b,c,d) I suppose.

The implicit declaration of `a` as local happens when the procedure is created, as opposed to when it is run. In your second example both `a` and `b` are implicitly declared local variables.

So, after the procedure is created, it is the same as if it were defined using,

   local a,b;

So the line

   b:=a+10;

will access the local instance of `a`.

It's unclear to me why you would expect 20 for the second example's returned value. (You didn't explain your reasoning.)

The following behaves just as you have originally described. The originally posted question was this:

   How can I create an even function in Maple—i.e. g(−x) = g(x)—so
   that g(−3) returns g(3), g(10) returns g(−10), and so on?

I'm, not sure that's really what you after, though. Did you want it to compute anything more concretely?

g:=proc(ee)                                        
   return 'procname'(-ee);                          
end proc:

g(-3);
                g(3)

g(10);
                g(-10)

g(-x);
                g(x)

g(x); 
                g(-x)

g(0);
                g(0)

 

Here are three ways to deal with the issue. They each require providing a numeric value for a, but they proceed in different ways. See the comments.

(I merged the plots only to keep the sheet short. You can do as you prefer, of course.)

restart

 

1) You can use the parameters option of the dsolve(...,numeric) command to
    specify that variable a does not yet have a numeric value.


This is more efficient if you are going to re-run the computations with

several choices of value for a. You do not need to call dsolve multiple

times in order to do this. Instead, you simply would call,

    S1(parameters=[somevalue]

 

I have renamed the previously assigned variable parameters to params.

 

inf := 50

lambda := L(t)+k*A(t)

equ1 := diff(S(t), t) = Lambda-(`&beta;__1`+`&beta;__2`)*S(t)*lambda-mu*S(t)+r*R(t)

equ2 := diff(E(t), t) = `&beta;__1`*S(t)*lambda-(alpha+mu+d__1+(1-p)*epsilon+p*eta)*E(t)

equ3 := diff(L(t), t) = p*eta*E(t)+e*A(t)-(`&sigma;__2`+mu+d__3+y__2)*L(t)

equ4 := diff(A(t), t) = (1-p)*epsilon*E(t)-(`&sigma;__1`+mu+d__2+y__1+e)*A(t)

equ5 := diff(R(t), t) = y__2*L(t)+y__1*A(t)-(mu+r)*R(t)

equ6 := diff(V(t), t) = `&beta;__2`*S(t)*lambda+alpha*E(t)+`&sigma;__1`*A(t)+`&sigma;__2`*L(t)-a*V(t)

S__0 := 1000000; E__0 := 1; L__0 := 0; A__0 := 0; R__0 := 0; V__0 := 0

ICs := S(0) = S__0, E(0) = E__0, L(0) = L__0, A(0) = A__0, R(0) = R__0, V(0) = V__0

params := [Lambda = .22, mu = .22, `&beta;__1` = 0.3e-1, `&beta;__2` = 0.3e-1, r = .5, alpha = .6, p = .1, k = 0.4e-1, y__1 = .8, y__2 = .8, d__1 = 0.2e-1, d__2 = 0.6e-1, d__3 = .2, `&sigma;__1` = 0.2e-1, `&sigma;__2` = 0.2e-1, epsilon = 0.2e-1, eta = 0.2e-1, e = 0.1e-1]

S1 := dsolve(eval([equ1, equ2, equ3, equ4, equ5, equ6, ICs], params), {A(t), E(t), L(t), R(t), S(t), V(t)}, range = 0 .. inf, numeric, 'parameters' = [a])

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := [a = a]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 6, (2) = 6, (3) = 0, (4) = 0, (5) = 1, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 1, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .0, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..7, {(1) = 0., (2) = 1., (3) = 0., (4) = 0., (5) = 1000000., (6) = 0., (7) = Float(undefined)})), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..6, {(1) = .1, (2) = .1, (3) = .1, (4) = .1, (5) = .1, (6) = .1}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0}, datatype = integer[8]), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = A(t), Y[2] = E(t), Y[3] = L(t), Y[4] = R(t), Y[5] = S(t), Y[6] = V(t)]`; YP[1] := 0.18e-1*Y[2]-1.11*Y[1]; YP[2] := 0.3e-1*Y[5]*(Y[3]+0.4e-1*Y[1])-.860*Y[2]; YP[3] := 0.2e-2*Y[2]+0.1e-1*Y[1]-1.24*Y[3]; YP[4] := .8*Y[3]+.8*Y[1]-.72*Y[4]; YP[5] := .22-0.6e-1*Y[5]*(Y[3]+0.4e-1*Y[1])-.22*Y[5]+.5*Y[4]; YP[6] := 0.3e-1*Y[5]*(Y[3]+0.4e-1*Y[1])+.6*Y[2]+0.2e-1*Y[1]+0.2e-1*Y[3]-Y[7]*Y[6]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = A(t), Y[2] = E(t), Y[3] = L(t), Y[4] = R(t), Y[5] = S(t), Y[6] = V(t)]`; YP[1] := 0.18e-1*Y[2]-1.11*Y[1]; YP[2] := 0.3e-1*Y[5]*(Y[3]+0.4e-1*Y[1])-.860*Y[2]; YP[3] := 0.2e-2*Y[2]+0.1e-1*Y[1]-1.24*Y[3]; YP[4] := .8*Y[3]+.8*Y[1]-.72*Y[4]; YP[5] := .22-0.6e-1*Y[5]*(Y[3]+0.4e-1*Y[1])-.22*Y[5]+.5*Y[4]; YP[6] := 0.3e-1*Y[5]*(Y[3]+0.4e-1*Y[1])+.6*Y[2]+0.2e-1*Y[1]+0.2e-1*Y[3]-Y[7]*Y[6]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..7, {(1) = 0., (2) = 0., (3) = 1., (4) = 0., (5) = 0., (6) = 1000000., (7) = 0.}); _vmap := array( 1 .. 6, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4), ( 5 ) = (5), ( 6 ) = (6)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, A(t), E(t), L(t), R(t), S(t), V(t)], (4) = [a = a]}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc


Of course we must still set a value for parameter a, before a numeric

solution can be obtained.

S1(parameters = [1])

[a = 1.]

P1 := plots:-odeplot(S1, [[t, S(t)]], 0 .. inf, linestyle = 1, color = blue, labels = ["time", "Population"], labeldirections = ["horizontal", "vertical"], legend = ["S(t) - Susceptible", "R(t) - Recovery"], axes = "boxed")

P2 := plots:-odeplot(S1, [[t, E(t)]], 0 .. inf, linestyle = 2, color = brown, labels = ["time", "Population"], labeldirections = ["horizontal", "vertical"], legend = ["S(t) - Susceptible", "R(t) - Recovery"], axes = "boxed")

P3 := plots:-odeplot(S1, [[t, L(t)]], 0 .. inf, linestyle = 3, color = red, labels = ["time", "Population"], labeldirections = ["horizontal", "vertical"], legend = ["S(t) - Susceptible", "R(t) - Recovery"], axes = "boxed")

P4 := plots:-odeplot(S1, [[t, A(t)]], 0 .. inf, linestyle = 4, color = green, labels = ["time", "Population"], labeldirections = ["horizontal", "vertical"], legend = ["S(t) - Susceptible", "R(t) - Recovery"], axes = "boxed")

P5 := plots:-odeplot(S1, [[t, R(t)]], 0 .. inf, linestyle = 5, color = purple, labels = ["time", "Population"], labeldirections = ["horizontal", "vertical"], legend = ["S(t) - Susceptible", "R(t) - Recovery"], axes = "boxed")

plots:-display(P1, P2, P3, P4, P5, axis[2] = [mode = log])

restart

 

2) Set a numeric value for a before calling dsolve(...,numeric).

 

This is simplest, but not most efficient if the computation is to be re-run

for multiple choises of value for a. In order to recompute for different

values for a the dsolve command could be repeated after each

re-assignment to a.

 

inf := 50

lambda := L(t)+k*A(t)

equ1 := diff(S(t), t) = Lambda-(`&beta;__1`+`&beta;__2`)*S(t)*lambda-mu*S(t)+r*R(t)

equ2 := diff(E(t), t) = `&beta;__1`*S(t)*lambda-(alpha+mu+d__1+(1-p)*epsilon+p*eta)*E(t)

equ3 := diff(L(t), t) = p*eta*E(t)+e*A(t)-(`&sigma;__2`+mu+d__3+y__2)*L(t)

equ4 := diff(A(t), t) = (1-p)*epsilon*E(t)-(`&sigma;__1`+mu+d__2+y__1+e)*A(t)

equ5 := diff(R(t), t) = y__2*L(t)+y__1*A(t)-(mu+r)*R(t)

equ6 := diff(V(t), t) = `&beta;__2`*S(t)*lambda+alpha*E(t)+`&sigma;__1`*A(t)+`&sigma;__2`*L(t)-a*V(t)

S__0 := 1000000; E__0 := 1; L__0 := 0; A__0 := 0; R__0 := 0; V__0 := 0

ICs := S(0) = S__0, E(0) = E__0, L(0) = L__0, A(0) = A__0, R(0) = R__0, V(0) = V__0

parameters := [Lambda = .22, mu = .22, `&beta;__1` = 0.3e-1, `&beta;__2` = 0.3e-1, r = .5, alpha = .6, p = .1, k = 0.4e-1, y__1 = .8, y__2 = .8, d__1 = 0.2e-1, d__2 = 0.6e-1, d__3 = .2, `&sigma;__1` = 0.2e-1, `&sigma;__2` = 0.2e-1, epsilon = 0.2e-1, eta = 0.2e-1, e = 0.1e-1]

a := 1

S1 := dsolve(eval([equ1, equ2, equ3, equ4, equ5, equ6, ICs], parameters), {A(t), E(t), L(t), R(t), S(t), V(t)}, range = 0 .. inf, numeric)