tomleslie

12175 Reputation

19 Badges

12 years, 176 days

MaplePrimes Activity


These are answers submitted by tomleslie

Unfortunately your system of equations is inconsistent, which is demonstrated in a couple of ways in the attached

restart

Eq[0, 0] := 1.33120000000000000000000000000*lambda[0, 1]+1.33120000000000000000000000000*lambda[0, 2]+1.33120000000000000000000000000*lambda[0, 3]+1.33120000000000000000000000000*lambda[0, 4] = .916487142969312002551492271668

Eq[0, 1] := 1.32901933598375616624661615670*lambda[0, 1]+1.32901933598375616624661615670*lambda[0, 2]+1.32901933598375616624661615670*lambda[0, 3]+1.32901933598375616624661615670*lambda[0, 4] = 1.09232395220587507357427365904

Eq[0, 2] := 1.37120000000000000000000000000*lambda[0, 1]+1.37120000000000000000000000000*lambda[0, 2]+1.37120000000000000000000000000*lambda[0, 3]+1.37120000000000000000000000000*lambda[0, 4] = 1.25415129307905065856083635281

Eq[0, 3] := .966980664016243833753383843299*lambda[0, 1]+.966980664016243833753383843299*lambda[0, 2]+.966980664016243833753383843299*lambda[0, 3]+.966980664016243833753383843299*lambda[0, 4] = 1.37114174964252179339832329224

Eq[0, 0]/(1.33120000000000000000000000000); Eq[0, 1]/(1.32901933598375616624661615670); Eq[0, 2]/(1.37120000000000000000000000000); Eq[0, 3]/(.966980664016243833753383843299)

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = .6884669043

 

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = .8219022270

 

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = .9146377574

 

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = 1.417961911

(1)

with(LinearAlgebra); A, b := GenerateMatrix([Eq[0, 0], Eq[0, 1], Eq[0, 2], Eq[0, 3]], indets(Eq[0, 0])); LinearSolve(A, b)

Matrix(%id = 36893488148293117220), Vector[column](%id = 36893488148293117100)

 

Error, (in LinearAlgebra:-BackwardSubstitute) inconsistent system

 

``

Download inconSys.mw

precisely what you want (or why!), Rather obviously, if

xi=xi-1+b

then the answer will be

{x__0, x__0 + b, x__0 + 2*b, x__0 + 3*b, x__0 + 4*b, x__0 + 5*b, x__0 + 6*b, x__0 + 7*b, x__0 + 8*b, x__0 + 9*b, x__0 + 10*b}

While you can generate this set recursively (as shown in the attached) it really seems like doing it the *hard way*!

  restart:
  n:=10:
  rProc:= proc(N, xv:=x__0)
               if   N=0
               then return xv;
               else return xv, rProc(N-1, xv+b)
               fi;
          end proc:

{rProc(n)};

{x__0, x__0+b, x__0+2*b, x__0+3*b, x__0+4*b, x__0+5*b, x__0+6*b, x__0+7*b, x__0+8*b, x__0+9*b, x__0+10*b}

(1)

 

Download recur.mw

in what context you want to do this - it just *seems* a bit weird - but the attached will do it

  restart:

  nArgs:= expr-> numelems([op([1..-1], expr)]);

  nArgs('f'(x1, x2, x3));
  nArgs('g'(x1, x2, x3, x4, x5, x6, x7));

proc (expr) options operator, arrow; numelems([op([1 .. -1], expr)]) end proc

 

3

 

7

(1)

 

Download args.mw

what you are trying to achieve here.- so the attached just contains some hints and tips which you might find useful - or maybe not!

There are many, ways in Maple to  perform multiple assignments - the attached shows a couple. In general I don't think it is a good idea to ask how to replicate the behaviour of a command in some other language. Much better to specify exactly what you want Maple to do.

The *simplest* way of getting "user input" is probably via the readline() command. Provided that you realise that anything entered in response to a  readline() command will be a *string*, and therefore you may have to "convert" appropriately, this approach can be surprisingly flexible. The attached shows how to do a simple variable assignment based on user input and also a simple guessing game, which I have shamelessly replicated from the  "help" for the readline() command.

So if these do not satisfy your requirement(s), then ithink yu need to specify your requirements more carefully without reference to other software packages.

  restart:
#
# Assigning multiple matrixes of the same size
#
  with(LinearAlgebra):
  randomize():
  n:=3:
  m:=5:
  A, B, C:= seq( RandomMatrix(n, m), j=1..3):
  A;
  B;
  C;

Matrix(3, 5, {(1, 1) = 2, (1, 2) = -22, (1, 3) = 47, (1, 4) = -44, (1, 5) = 31, (2, 1) = -37, (2, 2) = -1, (2, 3) = -96, (2, 4) = 21, (2, 5) = -90, (3, 1) = 6, (3, 2) = 90, (3, 3) = 27, (3, 4) = 45, (3, 5) = -46})

 

Matrix(3, 5, {(1, 1) = 91, (1, 2) = 74, (1, 3) = -55, (1, 4) = -91, (1, 5) = -52, (2, 1) = 65, (2, 2) = -46, (2, 3) = 28, (2, 4) = -2, (2, 5) = -3, (3, 1) = -92, (3, 2) = -17, (3, 3) = -70, (3, 4) = 96, (3, 5) = 68})

 

Matrix(%id = 36893488147886702828)

(1)

  restart:
#
# Initialize one matrix, then define several of the
# same size
#
  with(LinearAlgebra):
  M:=Matrix(3, 5);
  A, B, C:= seq( Matrix( op(1, M)), j=1..3):
  A;
  B;
  C;
  

Matrix(3, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0})

 

Matrix(3, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0})

 

Matrix(3, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0})

 

Matrix(%id = 36893488147952599156)

(2)

  restart:
#
# Really simple user input
#
 printf( "enter a value for the variable x\n"):
 x:= parse(readline()):
 x;

enter a value for the variable x

 

3.13

(3)

  restart;
#
# Using readline() to get user input
#
  GetNumber:= proc()
                   local number;
                   randomize();
                   number := rand( 1 .. 5 )();
                   printf( "Guess a number between 1 and 5:\n>> " );
                   while( parse( readline() ) <> number ) do
                          printf( "Try again!\n>> " );
                   end do:
                   printf( "Correct! The number was %d.\n", number );
              end proc:
   GetNumber();

Guess a number between 1 and 5:
>>

 

Try again!
>>

 

Try again!
>>

 

Correct! The number was 4.

 

 

 

Download tips.mw

 

 

without using the geometry() package.

  restart;
  with(plots):
  a := 11:
  b := 5:
  ell := (x*cos(theta) + y*sin(theta))^2/a^2 + (x*sin(theta) - y*cos(theta))^2/b^2 = 1:

  pltAll:=  proc(ang)
                 local p1,d;
                 uses plots, plottools:
                 p1:= implicitplot( eval(ell, theta=ang),
                                    x = -12 .. 12,
                                    y = -12 .. 12,
                                    thickness = 3,
                                    color = blue
                                  ):
                 d:= getdata(p1)[2]:
                 display( [ p1,
                            plot( [ [op(1, d[1]), op(1, d[2])],
                                    [op(1, d[1]), op(2, d[2])],
                                    [op(2, d[1]), op(2, d[2])],
                                    [op(2, d[1]), op(1, d[2])],
                                    [op(1, d[1]), op(1, d[2])]
                                  ],
                                  color=red
                                )
                          ]
                        );
            end proc:
  animate( display,
           ['pltAll'(theta)],
           theta=0..2*Pi
         );

 

 

Download sameAnim.mw

since this poster's questions have a habit of disappearing after they have been answered, I can come up with two possible animations which might be desired. Both are shown in the attached,

I'd guess the OP probably wants the first one (but my guesses are often wrong!)

  restart:
  with(plots):
  with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':

  point(OO, [0,0]):
  a := 11:
  b := 5:
  ellipse( ell, x^2/a^2+y^2/b^2=1):

  aProc:= proc(ang)
               local ellr,plt,d;
               uses geometry:
               rotation(ellr, ell, ang, counterclockwise, OO):
               plt:=draw(ellr):
               d:=plottools:-getdata(plt)[2];
               point(pt1,[op(1, d[1]), op(1, d[2])]):
               point(pt2,[op(1, d[1]), op(2, d[2])]):
               point(pt3,[op(2, d[1]), op(2, d[2])]):
               point(pt4,[op(2, d[1]), op(1, d[2])]):
               segment(s1, [pt1, pt2]):
               segment(s2, [pt2, pt3]):
               segment(s3, [pt3, pt4]):
               segment(s4, [pt4, pt1]):
               draw( [ ellr(color=blue, thickness=3),
                       s1(color=red),
                       s2(color=red),
                       s3(color=red),
                       s4(color=red)
                     ],
                      axes=none
                    ):    
          end proc:

  bProc:= proc(ang)
               uses geometry:
               rotation(ellr, ell, ang, counterclockwise, OO):
               point(p1, [ -a, -b]):
               point(p2, [ -a,  b]):
               point(p3, [  a,  b]):
               point(p4, [  a, -b]):
               segment( s1, [p1, p2]):
               segment( s2, [p2, p3]):
               segment( s3, [p3, p4]):
               segment( s4, [p4, p1]):
               rotation( s1r, s1, ang, counterclockwise, OO ):
               rotation( s2r, s2, ang, counterclockwise, OO ):
               rotation( s3r, s3, ang, counterclockwise, OO ):
               rotation( s4r, s4, ang, counterclockwise, OO ):
               draw( [ ellr(color=blue, thickness=3),
                       s1r(color=red),
                       s2r(color=red),
                       s3r(color=red),
                       s4r(color=red)
                     ],
                     axes=none
                   ):
         end proc:

  animate( display,
           ['aProc'(theta)],
           theta=0..2*Pi
         );
  animate( display,
           [  ['bProc'(theta)]],
           theta=0..2*Pi
         );

               

 

 

 

Download anims.mw

 

your pde is second order in time and second order in space so a total of four boundary/initial conditions would be needed for a complete solution. One could (in principle) get an analytic solution containing an arbitrary constant. But since Maple cannot find an analytic solution, then the only option is to solve numerically. For a numeric solution four boundary/initial condition is mandatory - so I made one up (higlighted in blue in the attached).

Now a numeric solution is possible - see the attached, which for some reason is refusing to display inline on this site!

Download pdSolve.mw

that there is a "feasible region" whihc satisfies all four inequalities - see the attached

  ineqs:=[  alpha <= 0.0002500000000*(-18000.*m^2 + 47271.*m + 39514. + sqrt(3.24000000*10^8*m^4 - 1.701756000*10^9*m^3 - 4.266980559*10^9*m^2 - 3.036299412*10^9*m - 6.95987804*10^8))/(9.*m^2 + 12.*m + 4.),
  0.00005000000000*(-90000.*m^2 + 237237.*m + 198158. + sqrt(8.100000000*10^9*m^4 - 4.270266000*10^10*m^3 - 1.069976858*10^11*m^2 - 7.612670111*10^10*m - 1.744924704*10^10))/(9.*m^2 + 12.*m + 4.) <= alpha,
  -0.6666666667 < m,
   m < -0.6665522013
  ];

[alpha <= 0.2500000000e-3*(-18000.*m^2+47271.*m+39514.+(324000000.0*m^4-1701756000.*m^3-4266980559.*m^2-3036299412.*m-695987804.0)^(1/2))/(9.*m^2+12.*m+4.), 0.5000000000e-4*(-90000.*m^2+237237.*m+198158.+(8100000000.*m^4-0.4270266000e11*m^3-0.1069976858e12*m^2-0.7612670111e11*m-0.1744924704e11)^(1/2))/(9.*m^2+12.*m+4.) <= alpha, -.6666666667 < m, m < -.6665522013]

(1)

  p1:=plots:-inequal( ineqs[1], m=lhs(ineqs[3])..rhs(ineqs[4]), alpha=15000..200000, color=red):
  p2:=plots:-inequal( ineqs[2], m=lhs(ineqs[3])..rhs(ineqs[4]), alpha=15000..200000, color=blue):
  plots:-display([p1,p2]);

 

 

Download ineqs.mw

but maybe somethig in the attached will help?

  restart:
  exact:=u(x,t)=sin(x)*cos(t);
  approx:=u(x,t)=sin(x)*sin(t);
  plot3d( [rhs(exact), rhs(approx)], x=0..2*Pi, t=0..3*Pi, color=[red, blue], style=surface);

u(x, t) = sin(x)*cos(t)

 

u(x, t) = sin(x)*sin(t)

 

 

 

Download plt3.mw

 

if I have to do it the hard way (why?), then it would probably look something like the attached

  restart:
  with(plottools):
  with(plots):
  _local(D):

  a:= 7:
  b:= a/2:
  r:= b/2:
  c:= 2*sqrt(10)*r:

  Cir1:= disk([c/2 - r, -c/2 + r], r, color = blue):
  Cir2:= disk([-c/2 + r, -c/2 + r], r, color = blue):
  Cir3:= disk([-c/2 + r, c/2 - r], r, color = blue):
  Cir4:= disk([c/2 - r, c/2 - r], r, color = blue):

  rot:= Matrix(2, 2, [[cos(theta), -sin(theta)], [sin(theta), cos(theta)]]):
  T:= rot.(Vector(2, [x, y])):
  Ell:= x^2/a^2 + y^2/b^2 = 1:
  doRot:= ang -> expand(eval(Ell, eval([x = T[1], y = T[2]], theta = ang))):
  e1:= doRot(Pi/4):
  e2:= doRot(-Pi/4):

  p1:= shadebetween( rhs(solve(e2,[y])[2][]), -c/2, x=-c/2..0, scaling=constrained):
  p2:= display([p1, reflect(p1, [[ 0, -c/2], [0, c/2]])]):
  p3:= display([p2, reflect(p2, [[ -c/2, 0], [c/2, 0]])]):
  p4:= display([p3, rotate(p3, Pi/2, [0,0])], color=red):

  Ell1:= implicitplot(e1, x = -6 .. 6, y = -6 .. 6, thickness = 3, color = blue):
  Ell2:= implicitplot(e2, x = -6 .. 6, y = -6 .. 6, thickness = 3, color = blue):

  display([p4, Ell1, Ell2, Cir1, Cir2, Cir3, Cir4], scaling = constrained, axes = none);

 

 

 

Download hardWay.mw

 

when you specify the range o the independent variable as 0..3*Pi, then the calculation is (probably?) done in steps which are a fraction of 'Pi' - whihc means it is very likely to exactly "land on" the value Pi, at which point you have reached the limit specified for the range of the independent variable, so the calculation stops.

Two workarounds

  1. Don't supply a range for the dependent variable - for many ODEs you may not be able work this out a priori in any case
  2. Or make sure that the specified range of the dependet variable is "larger" than any value it can possibly achieve

Both methods are shown in the attached

  restart;
  interface(version);
  ode := diff(y(t), t$2) + y(t)=0;
#
# No specified range for the dependent variable
#
  DEtools:-DEplot( ode,
                   y(t),
                   t=0 .. 3*Pi,
                   [[y(0)=1, D(y)(0)=0]],
                   linecolor=blue
                 );
#
# "Large" specified range for the dependent variable
#
  DEtools:-DEplot( ode,
                   y(t),
                   t=0 .. 3*Pi,
                   y=-1.1..1.1,
                   [[y(0)=1, D(y)(0)=0]],
                   linecolor=blue
                 );

`Standard Worksheet Interface, Maple 2022.1, Windows 7, May 26 2022 Build ID 1619613`

 

diff(diff(y(t), t), t)+y(t) = 0

 

 

 

 

Download deplot.mw

to use the Maple's geometry() package. One advantage of this approach is that its draw() command has the concept of a stacking order for geometric objects. So, for example, in the attached

  1. the red square is at the bottom of the stack
  2. then two filled (with color=white) ellipses are placed on top of the red square
  3. then two ellipse "outlines" are placed on top of 1+2 above
  4. then four filled circles are placed on top of 1+2+3 above.

See the attached.

(As usual, the figure renders rather better in a Maple worksheet than it does on this site- honest!)

  restart;
  with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':
  point(OO,0,0):

  a := 7:
  b := a/2:
  r := b/2:
  c := 2*sqrt(10)*r:
#
# The "background" square
#
  square( Sq,
          [ point(sq1, -c/2, -c/2),
            point(sq2, -c/2,  c/2),
            point(sq3,  c/2,  c/2),
            point(sq4,  c/2, -c/2)
          ]
        ):
#
# The four circles
#
  circle( cir1,
          [ point(c1, c/2 - r, -c/2 + r),
            r
          ]
        ):
  circle( cir2,
          [ point(c2, c/2 - r,  c/2 - r),
            r
          ]
        ):
  circle( cir3,
          [ point(c3, -c/2 + r, -c/2 + r),
            r
          ]
        ):
  circle( cir4,
          [ point(c4, -c/2 + r, c/2 - r),
            r
          ]
        ):
#
# The two ellipses
#
  ellipse(e1, x^2/a^2 + y^2/b^2 = 1):
  rotation(e2, e1, Pi/4, clockwise, OO):
  rotation(e3, e1, Pi/4, counterclockwise, OO):
#
# Draw everything - note that the "stacking"
# order is important
#
  draw( [ cir1( filled=true, color=blue),
          cir2( filled=true, color=blue),
          cir3( filled=true, color=blue),
          cir4( filled=true, color=blue),
          e2( color=blue, thickness = 3),
          e3( color=blue, thickness = 3),
          e2( filled=true, color=white),
          e3( filled=true, color=white),
          Sq( filled=true, color=red)
        ],
        axes=none,
        size=[800, 800],
        gridlines=false
      );
 

 

 

Download geomDraw.mw

 

So you wantt a solution for a system of odes which contains some totally arbitrary function h(t)?

Just for amusement in the attached, I have shown that Maple will provide solutions for a variety of possibilities of this  otherwise completely indeterminate function - I used h(t)=0, h(t)=cos(t) and h(t)=t^2. In all cases, a solution is found. But for h(t) entirely arbitrary - get real!

You might want to consider the attached where I have used some more-or-less  *random* guesses for h(t), as in   h(t)=0,  h(t)=cos(t), and h(t)=t^2. All of these generate "solutions". But for h(t) entirely arbitrary - no chance!

  restart:
  assume(lambda>0);
  interface(showassumed=0):
  odeSys:= diff(x(t),t)=y(t),
           diff( y(t), t) =z(t),
           diff(z(t),t)=-lambda*y(t)-h(t);
  dsolve(eval( [odeSys], h(t)=0));
  dsolve(eval( [odeSys], h(t)=cos(t)));
  dsolve(eval( [odeSys], h(t)=t^2));
  

diff(x(t), t) = y(t), diff(y(t), t) = z(t), diff(z(t), t) = -lambda*y(t)-h(t)

 

{x(t) = (_C1*lambda^(1/2)-cos(lambda^(1/2)*t)*_C2+sin(lambda^(1/2)*t)*_C3)/lambda^(1/2), y(t) = _C2*sin(lambda^(1/2)*t)+_C3*cos(lambda^(1/2)*t), z(t) = lambda^(1/2)*(cos(lambda^(1/2)*t)*_C2-sin(lambda^(1/2)*t)*_C3)}

 

{x(t) = (_C1*lambda^(3/2)-lambda*cos(lambda^(1/2)*t)*_C3+lambda*sin(lambda^(1/2)*t)*_C2-sin(t)*lambda^(1/2)-_C1*lambda^(1/2)+_C3*cos(lambda^(1/2)*t)-_C2*sin(lambda^(1/2)*t))/(lambda^(1/2)*(lambda-1)), y(t) = sin(lambda^(1/2)*t)*_C3+cos(lambda^(1/2)*t)*_C2-cos(t)/(lambda-1), z(t) = (lambda^(1/2)*cos(lambda^(1/2)*t)*_C3*(lambda-1)-lambda^(1/2)*sin(lambda^(1/2)*t)*_C2*(lambda-1)+sin(t))/(lambda-1)}

 

{x(t) = -_C3*cos(lambda^(1/2)*t)/lambda^(1/2)+_C2*sin(lambda^(1/2)*t)/lambda^(1/2)-(1/3)*t^3/lambda+_C1+2*t/lambda^2, y(t) = (sin(lambda^(1/2)*t)*_C3*lambda^2+cos(lambda^(1/2)*t)*_C2*lambda^2-lambda*t^2+2)/lambda^2, z(t) = (lambda^(3/2)*cos(lambda^(1/2)*t)*_C3-lambda^(3/2)*sin(lambda^(1/2)*t)*_C2-2*t)/lambda}

(1)

 

Download badProb.mw

 

but probably the simplest is to use the surfdata() command as shown in the attached. Many options are available with the surfdata() command, none of which I have bothered with

Two points to note

  1. You will have to change the file path in th ExcelTools:-Import() command to something appropriate for your machine
  2. It is not obvious from the data yu supply which direction you regard as 'x' and which is 'y', so I madee an arbitrary choice

   restart;
   with(plots):
   f := ExcelTools:-Import("C:/Users/TomLeslie/Desktop/PlotTest.xlsx"):
   L:= [ seq
         ( [ seq
             ( [ f[i+1,1], f[1,j+1], f[i+1, j+1] ],
                 i=1..9
             )
           ],
           j=1..9
         )
       ]:
   surfdata(L);

 

 

NULL

Download surfplot.mw

to get a numerical solution for your problem - although it appears that Maple cannot find an analytic solution.

See the attached.

After the point of achieving a solution, I have no idea what your subsequent code is trying to achieve, so I have ignored it!

  restart:
  with(plots):
  _local(I):
  Digits := 15:
  odeSys:= (1 - p)*(diff(S(t), t) + mu*S(t)) + p*(diff(S(t), t) + mu*S(t) + beta*S(t)*I(t) - rho*R(t) - varepsilon),
           (1 - p)*(diff(E(t), t) + (alpha1 + mu)*E(t)) + p*(diff(E(t), t) + (alpha1 + mu)*E(t) - beta*S(t)*I(t)),
           (1 - p)*(diff(I(t), t) + (alpha2 + delta + mu)*I(t)) + p*(diff(I(t), t) + (alpha2 + delta + mu)*I(t) - alpha1*E(t)),
           (1 - p)*(diff(R(t), t) + (mu + rho)*R(t)) + p*(diff(R(t), t) + (mu + rho)*R(t) - alpha2*I(t)):
  params:= [ mu=0.133*10^(-5),
             varepsilon=0.99879,
             delta=0.004554,
             beta=0.1009*10^(-6),
             alpha1=0.0008999,
             alpha2=0.1997,
             rho=0.00090021,
             p=1
           ]:
  ibvc := S(0) = 2304219, E(0) = 84929, I(0) = 299, R(0) = 71411:
#
# Won't find an analytic solution :-(
#
  solA:=dsolve(eval( [odeSys ,ibvc], params));
#
# Will find a numeric solution :-)
#
  solN:=dsolve(eval( [odeSys,ibvc], params), numeric);
  odeplot(solN, [t, S(t)], t=0..10, color=[red, green, blue, black]);
  odeplot(solN, [t, E(t)], t=0..10, color=[red, green, blue, black]);
  odeplot(solN, [t, I(t)], t=0..10, color=[red, green, blue, black]);
  odeplot(solN, [t, R(t)], t=0..10, color=[red, green, blue, black]);

I

 

Warning, The imaginary unit, I, has been renamed _I

 

 

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 := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 28, [( 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..65, {(1) = 4, (2) = 4, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (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, (64) = -1, (65) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .9832164330847996, (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..4, {(1) = 84929.0, (2) = 299.0, (3) = 71411.0, (4) = 2304219.0}, datatype = float[8], order = C_order)), ( 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..4, {(1) = .1, (2) = .1, (3) = .1, (4) = .1}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 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}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8]), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..4, {(1) = 84929.0, (2) = 299.0, (3) = 71411.0, (4) = 2304219.0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = -7.024349237099997, (2) = 15.35526343, (3) = -4.669572940000002, (4) = -7.2971383928999956}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = E(t), Y[2] = I(t), Y[3] = R(t), Y[4] = S(t)]`; YP[1] := -0.901230000000000e-3*Y[1]+0.100900000000000e-6*Y[4]*Y[2]; YP[2] := -.204255330000000*Y[2]+0.8999e-3*Y[1]; YP[3] := -0.901540000000000e-3*Y[3]+.1997*Y[2]; YP[4] := -0.133000000000000e-5*Y[4]-0.100900000000000e-6*Y[4]*Y[2]+0.90021e-3*Y[3]+.99879; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = E(t), Y[2] = I(t), Y[3] = R(t), Y[4] = S(t)]`; YP[1] := -0.901230000000000e-3*Y[1]+0.100900000000000e-6*Y[4]*Y[2]; YP[2] := -.204255330000000*Y[2]+0.8999e-3*Y[1]; YP[3] := -0.901540000000000e-3*Y[3]+.1997*Y[2]; YP[4] := -0.133000000000000e-5*Y[4]-0.100900000000000e-6*Y[4]*Y[2]+0.90021e-3*Y[3]+.99879; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 27 ) = (""), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0), ( 28 ) = (0)  ] ))  ] ); _y0 := Array(0..4, {(1) = 0., (2) = 84929., (3) = 299., (4) = 71411.}); _vmap := array( 1 .. 4, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4)  ] ); _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); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `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 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 elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) 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, E(t), I(t), R(t), S(t)], (4) = []}); _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

 

 

 

 

 

 

Download odeProb.mw

1 2 3 4 5 6 7 Last Page 1 of 187