dharr

Dr. David Harrington

6115 Reputation

21 Badges

19 years, 329 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are replies submitted by dharr

@srmaple Actually, I was comparing with the beta__c=0 solution, not Mathematica's, and suggested c2 = c2*beta__c as a simple function that went to zero. I even suggested c2*beta__c^2 was OK. But perhaps this is a better justification?

restart;

with(DEtools):

ode1 := diff(c(T),T)=-2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2);
ode2 := diff(p__c(T),T)=2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T);
sys := [ode1,ode2]:

diff(c(T), T) = -2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2)

diff(p__c(T), T) = 2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T)

First solve the system with beta__c = 0.
We want the solution for non-zero beta__c to reduce to this in the limit of beta__c -> 0.

sys0:=eval(sys,beta__c=0);
dsolve(sys0);

[diff(c(T), T) = -2*c(T), diff(p__c(T), T) = 2*p__c(T)]

{c(T) = c__2*exp(-2*T), p__c(T) = c__1*exp(2*T)}

sol2:=simplify([eval(dsolve(sys,[c(T),p__c(T)],explicit),c__2=f(beta__c))]) assuming beta__c>0;

[{c(T) = -2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = -2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = -2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = -2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}]

So p__c(T) already tends to the right solution for f(beta__c) -> 0.

What other condition on f(beta) is required for c(T)? Take the first solution for example

cT1:=eval(c(T),sol2[1]);

-2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2)

Take a Taylor series for f(beta__c) with f(0)=0

ff:=eval(series(f(beta__c),beta__c),f(0)=0);

 

series((D(f))(0)*beta__c+((1/2)*((D@@2)(f))(0))*beta__c^2+((1/6)*((D@@3)(f))(0))*beta__c^3+((1/24)*((D@@4)(f))(0))*beta__c^4+((1/120)*((D@@5)(f))(0))*beta__c^5+O(beta__c^6),beta__c,6)

series(eval(cT1,f(beta__c)=ff),beta__c,2) assuming beta__c::positive;
eval(convert(%,polynom),beta__c=0); #first term

series(-2*2^(1/4)*((D(f))(0)/(c__1*exp(8*T))^(1/2))^(1/2)+((1/2)*2^(1/4)*((D(f))(0)/(c__1*exp(8*T))^(1/2))^(1/2)*(-((D@@2)(f))(0)*c__1*exp(8*T)+8*(D(f))(0)^2)/(c__1*exp(8*T)*(D(f))(0)))*beta__c+O(beta__c^2),beta__c,2)

-2*2^(1/4)*((D(f))(0)/(c__1*exp(8*T))^(1/2))^(1/2)

Above agrees with the beta__c=0 solution for any function of beta__c that goes to zero and whose first derivative is non-zero at beta__c=0, e.g., const*beta__c

NULL

Download dsolve.mw

@Ronan A couple of comments

1. I wasn't using makehelp, but see here for how to do subfolders with that method

Creating Help - An Exploration Of Makehelp The browser entries have 3 items. But notice that he did not succeed in altering the order.

2. I was supposing that you might just replace the TOC with a new one, by making the appropriate Record structure. In my case I wanted my Orbitals package under Physics, so that the the top level Record is for Physics, and it has a child Record of Orbitals, which has child Records for the individual help pages. Here's what the final TOC looks like, from lprint(eval(TOC)):

Record('entry' = Physics,'topic' = (NULL),'priority' = (NULL),'language' = "en"
,'product' = "Maple",'category' = "Help Page",'children' = [Record('entry' = 
Orbitals,'topic' = (NULL),'priority' = (NULL),'language' = "en",'product' = 
"Maple",'category' = "Help Page",'children' = [Record('entry' = "Overview",'
topic' = "Orbitals,Overview",'priority' = 100,'language' = "en",'product' = 
"Maple",'category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(atom4e, atomJ, atomK) - atomic exchange, coulomb and four-electron integrals"
,'topic' = "Orbitals,atom4e",'priority' = 70,'language' = "en",'product' = 
"Maple",'category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(cartesian, fullcartesian) - convert to cartesian coordinates",'topic' = 
"Orbitals,cartesian",'priority' = 65,'language' = "en",'product' = "Maple",'
category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(HYD) - hydrogenic orbital",'topic' = "Orbitals,HYD",'priority' = 60,'language
' = "en",'product' = "Maple",'category' = "Help Page",'children' = (NULL)), 
Record('entry' = "(HYDr) - radial part of hydrogenic orbital",'topic' = 
"Orbitals,HYDr",'priority' = 55,'language' = "en",'product' = "Maple",'category
' = "Help Page",'children' = (NULL)), Record('entry' = 
"(overlap) - overlap integral",'topic' = "Orbitals,overlap",'priority' = 50,'
language' = "en",'product' = "Maple",'category' = "Help Page",'children' = (
NULL)), Record('entry' = "Plotting orbitals examples",'topic' = 
"Orbitals,plots",'priority' = 5,'language' = "en",'product' = "Maple",'category
' = "Help Page",'children' = (NULL)), Record('entry' = 
"(realY) - real combinations of spherical harmonics",'topic' = "Orbitals,realY"
,'priority' = 45,'language' = "en",'product' = "Maple",'category' = "Help Page"
,'children' = (NULL)), Record('entry' = "Orbital package references",'topic' =
"Orbitals,references",'priority' = 2,'language' = "en",'product' = "Maple",'
category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(STO) - Slater type orbital",'topic' = "Orbitals,STO",'priority' = 35,'
language' = "en",'product' = "Maple",'category' = "Help Page",'children' = (
NULL)), Record('entry' = "(STOr) - radial part of Slater type orbital",'topic'
= "Orbitals,STOr",'priority' = 30,'language' = "en",'product' = "Maple",'
category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(Y) - spherical harmonic",'topic' = "Orbitals,Y",'priority' = 25,'language' =
"en",'product' = "Maple",'category' = "Help Page",'children' = (NULL)), Record(
'entry' = "Orbitals package examples",'topic' = "Orbitals,examples",'priority'
= 10,'language' = "en",'product' = "Maple",'category' = "Help Page",'children'
= (NULL))])])

If you think it's useful, I could send the help files themselves, so you could try converting them to the help database.

@MaPal93 Your treatment neglects the fact that the relationship between your new Sigma_v and sigma__v1 and sigma__v2 is lost if you only proceed with p3 and don't add the new auxilliary equation to the analysis. I fixed this and the scale factor up.

For two equations convert each to polynomials by removing the radicals. Then each polynomial is normalized so the first terms is 1 and then all the remaining terms for the two polynomials and the auxilliary equation(s)  are collected (even if there is a lot of redundancy, the matrix procedure will take care of it.)

Case_with_radicals.mw

@MaPal93 

1. It doesn't matter because const*p(x) has the same roots as p(x). If simplify didn't simplify to that form, then the old variables wouldn't be removed. I suppose checking for removal of the old variables could be automated, but I was just looking at it to check all the old ones were gone and did that manually to get a form where they are all gone.

2. Hubert's theory only works for polynomials where the monomials have integer exponents, or rational functions with the same condition. You could make new variables for the quantities with square roots, and then apply the theory to the system of equations, e.g. defining s= sqrt(a^2+b^2) will add terms a^2/s^2 and b^2/s^2 and then use the original eqn with s in. 

@sursumCorda I usually don't like to use undocumented commands, but that does work nicely for both rectangular and sparse matrices (by sparse I mean storage=sparse vs storage=rectangular).

Is the list l always the sequence of positive integers? Surely it must go to 10 if there are 10 twins and 10 is in your final list? If both these are true, I get 7 in the final list, so I'm not sure I understand. See attached.

twins.mw

@MaPal93 Since plot will evaluate the expressions numerically anyway, you don't want to simplify. If you plot [allvalues(Lambda)] it does plot, but it is hard to see which surface is which. You can plot them 1 at a time (I updated the above to Maple 2023 with a plot of the first one). If they are mostly complex, there won't be much to see.

The warning is because of the 5x51 matrix K output, but it is fast to rerun the worksheet so you can answer yes or no to "save these results", or just put a colon to not display the K matrix.

I'm guessing that since p3_complex has the same (by eye) number of real variables as p3, the new Lambda probably depends also on two nondim parameters. For p3_mostcomplex there is one more real variable, so maybe one more non-dim, but you might be lucky.

@Ronan My survey of the various filetypes used (.off, .ply, .stl etc) is that they are all for 3D objects and do not contain other information, with the possible exception of polygon face "textures" and "materials", which are for gaming applications and 3D printing respectively. So I think this is not a Maple issue. The face coloring I am less certain about, and the ability to support color is different for different files. But Maple doesn't seem to attempt to export color for any of them, which may be just that most 3D objects have a single color and not separate colors for individual faces. That can be fixed by messing with the PLOT3D structure, but that is significant effort.  

@dharr @acer's routine is brilliant for a single sum, but does not take advantage of efficiencies if you need a sequence of sums. @sursumCorda's approach can be adapted to give an efficient way to generate a sequence.

Download ijknew.mw

[Edit: faster routine now added:]

restart;

Modification of @sursumCorda's routine that works out a term for each permutation of a partition.

This version  (i) accumulates the sums and partitions as we go (`combstruct/allstructs/Partition/allpartk` uses option remember), and

                      (ii) does not work out all the permutations, just takes one term for each partition and multiplies by the number of permutations

sums := proc(m::posint)
  local n, s, parts, part, ss;
  ss:=Vector(m);        
  s := 0:
  for n to m do
    parts := `combstruct/allstructs/Partition/allpartk`(n,3); # partitions of n with 3 parts
    s := s + add(combinat:-numbperm(part)/mul(part), part in parts);
    ss[n] := s;
  end do;
  ss
end proc:

L:=CodeTools:-Usage(sums(200),iterations=4):

memory used=357.30MiB, alloc change=72.00MiB, cpu time=406.25ms, real time=3.39s, gc time=23.44ms

Further efficiencies are possible by generating information for each partition iteratively only once..
The partitions themselves do not have to be generated.
If p(n, k) are the partitions of p with k parts, then these may be found by

(i) appending 1 to each of the p(n-1, k-1)

    The product of the parts is unchanged.
    The number of 1's increases by 1
    The number of permutations is multiplied by (np+1)/(m1+1), where m1 is the number of 1's in the old partition,
       and np is the number of parts in the old partition.
     For k = 3 quantityp2 = a*b+a*c+b*c can be calculated as prod+sum for [a, b]

and

(ii) adding 1 to each part in each of the p(n-k, k)

    The new product is prod+p2+sum+1, with p2 = 0 for k = 2 and p2 = a*b+a*c+b*c for [a, b, c] k = 3
     The new p2 is p2+2*sum+3
    The number of 1's is now zero.

    The number of permutations is unchanged.

sums2:=proc(m::posint)
        description "outputs length m Vector with V[n] = sum(1/(i*j*k), 1 <= (i,j,k) <= n)";
        local prods, m1, nperms, p2, n, s, ss, key, nn, nm1, nm2, nm3;
        if m = 1 then return Vector([0]) end if;
        if m = 2 then return Vector([0, 0]) end if;
        if m = 3 then return Vector([0, 0, 1]) end if;
        # initialize Array
        # prods = products of partition parts, e.g. (2,2) for 2 in 2 parts = [1,1]; prod =1
        prods := Array(1..4, 2..3, {
                (1,2) = Vector(),    (1,3) = Vector(),
                (2,2) = Vector([1]), (2,3) = Vector(),
                (3,2) = Vector([2]), (3,3) = Vector([1])});
        # multiplicity of 1's, e.g., (2,2) has 2 1's
        m1 := Array(1..4, 2..3, {
                (1,2) = Vector(),    (1,3) = Vector(),
                (2,2) = Vector([2]), (2,3) = Vector(),
                (3,2) = Vector([1]), (3,3) = Vector([3])});
        # number of permutations, e.g. 2 ways for (3,2) = [1,2] or [2,1]
        nperms := Array(1..4, 2..3, {
                (1,2) = Vector(),    (1,3) = Vector(),
                (2,2) = Vector([1]), (2,3) = Vector(),
                (3,2) = Vector([2]), (3,3) = Vector([1])});
        # for partitions with 3 parts only, p2 is sum of (products of two parts)
        # for [a,b,c], p2 = a*b + b*c +a*c        
        p2 := Array(1..4, {(1) =  Vector(), (2) =  Vector(), (3) =Vector([3])});
        s := 1;
        ss := Vector(m, {1 = 0, 2 = 0, 3 = s});
        # key has pointers to Array rows; reindexing allows reuse of rows in rotationg fashion
        key := [1, 2, 3, 4];        
        for n from 4 to m do
                nm3, nm2, nm1, nn := key[]; # mn3 = n-3 etc
                # update Array for this row
                # first part of Vector is for appending 1 to p(n-1,k-1)
                # second part of Vector is for adding 1's to p(n-k,k)
                # k = 2;
                        # (n-1,1) = [n-1] -> [n-1,1]; adding 1's -> no 1's
                        m1[nn, 2] := Vector(numelems(m1[nm2, 2])+1, {1=1});         # [1, 0, 0, ...]
                        # [n-1,1] has 2 perms; adding 1's leaves nperms unchanged             
                        nperms[nn, 2] := Vector([2, nperms[nm2, 2]]);
                        # [n-1,1] has prod n-1; adding 1's gives prodnew = (prod + sum + 1) of p(n-2,2)
                        prods[nn, 2] := Vector([n-1, prods[nm2, 2] +~ (n-1)]);                        
                # k = 3;
                        # appending 1 increases m1 by 1; adding 1's -> no 1's
                        m1[nn, 3] := Vector([ m1[nm1, 2] +~ 1, Vector(numelems(m1[nm3, 3])) ]);        
                        # appending 1: new nperms = (old nperms)*k/(m1+1); adding 1's leaves nperms unchanged                          
                        nperms[nn, 3] := Vector([ zip((x,y) -> x*3/(y+1), nperms[nm1, 2], m1[nm1, 2]), nperms[nm3, 3]]);
                        # [a,b] -> [a,b,1], p2 = (a*b) + (a+b); [a+1,b+1,c+1] has p2 = (a*b + b*c + a*c) + 2*(a+b+c) + 3
                        p2[nn] := Vector([ prods[nm1, 2] +~ (n-1), p2[nm3] +~ (2*n-3) ]);
                        # [a,b,1] has same prod as [a,b]; adding 1's gives prodnew = prod + p2 + sum + 1
                        prods[nn, 3] := Vector([ prods[nm1, 2], prods[nm3, 3] +~ p2[nm3] +~ (n-2) ]);
                # calculate and store sum
                s := s + add(nperms[nn, 3] /~ prods[nn, 3]);
                ss[n] := s;
                # reindex Array rows for next iteration
                key := ListTools:-Rotate(key, 1);
        end do;
        ss
end proc:

L2:=CodeTools:-Usage(sums2(200),iterations=4):

memory used=85.29MiB, alloc change=12.27MiB, cpu time=78.25ms, real time=452.50ms, gc time=3.91ms

EqualEntries(L,L2);

true

NULL

NULL

Download ijknew2.mw

@raj2018 There are two pieces becouse there are two solutions for each kc and deltah, e.g., for kc=2.3, deltah = 0.8, M can be 0.0007832833567 or 1.134283416. If you don't want to see the bottom piece on the plot, just adjust the M range accordingly.

To get data out, it is best just to solve the equation for the values you want. An example is given in the attached 

implicitplot3d.mw

@sursumCorda If you know the answer, you can cheat. Both convert(z*hypergeom([1/4, 1/2], [5/4], z**4), FormalPOwerSeries) and convert(EllipticF(z, I), FormalPowerSeries) give the same power series, verifying they are the same.

They are converted to apparently different expressions with integrals with convert(..., Int), so interconverting might be possible by rearranging, but it wasn't obvious to me.

@acer Nice solution! I'm a bit surprised the Iterator method is so much slower than @sursumCorda's more-or-less equivalent method with combinat routines.

Just to add to @mmcdara's comments. I especially agree with his last comment. Please clarify that you want to find each of  the 12 variables n1,n,2,n3,n4,4,a1,a2,a3,a4,s1,s2,s3,s4 in terms of the 4 parameters b,q,r,theta. So as noted you need 12 equations for your 12 variables. 

The solutions for b=0 have arbitrary values of n3,n4 and s4. Do you expect this to be true in the general case? A dimensional analysis might help you to reduce the number of variables or parameters.

The reason it works for b=0 is that it can be cast as a polynomial system for that case. That can be done also for the full system, which might make it feasible to solve. But it would be better to simplify first, and we definitely need number of eqns = number of variables. Are the ai, ni and si components of 4-vectors? If so, perhaps the original equation(s) in matrix vector form would be useful.

@dharr And here's how to do it with an event.

restart;

k12:=0.0452821;
k21:=0.0641682;
k1x:=0.00426118;
Dose:=484;
#Inj:=ifelse(t=1,Dose,0);

ode_sys:=diff(SA2(t),t)=SA1(t)*k12-SA2(t)*k21,
         diff(SA1(t),t)=-k1x*SA1(t);

ics:=SA2(0)=0,SA1(0)=0;

0.452821e-1

0.641682e-1

0.426118e-2

484

diff(SA2(t), t) = 0.452821e-1*SA1(t)-0.641682e-1*SA2(t), diff(SA1(t), t) = -0.426118e-2*SA1(t)

SA2(0) = 0, SA1(0) = 0

Solution1:=dsolve({ode_sys,ics},{SA1(t),SA2(t)},type=numeric,events=[[t=1,SA1(t)=SA1(t)+Dose]]);

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 ) = ([Array(1..2, 1..21, {(1, 1) = 1.0, (1, 2) = 2.0, (1, 3) = 2.0, (1, 4) = .0, (1, 5) = 1.0, (1, 6) = .0, (1, 7) = 1.0, (1, 8) = undefined, (1, 9) = undefined, (1, 10) = 1.0, (1, 11) = undefined, (1, 12) = undefined, (1, 13) = undefined, (1, 14) = undefined, (1, 15) = undefined, (1, 16) = undefined, (1, 17) = undefined, (1, 18) = undefined, (1, 19) = undefined, (1, 20) = undefined, (1, 21) = undefined, (2, 1) = 1.0, (2, 2) = .0, (2, 3) = 100.0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = undefined, (2, 9) = undefined, (2, 10) = 0.10e-6, (2, 11) = undefined, (2, 12) = .0, (2, 13) = undefined, (2, 14) = .0, (2, 15) = .0, (2, 16) = undefined, (2, 17) = undefined, (2, 18) = undefined, (2, 19) = undefined, (2, 20) = undefined, (2, 21) = undefined}, datatype = float[8], order = C_order), proc (t, Y, Ypre, n, EA) EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) if e = 1 then Y[1] := Y[1]+484 end if; return 0 end proc, Array(1..1, 1..2, {(1, 1) = undefined, (1, 2) = undefined}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..65, {(1) = 2, (2) = 2, (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) = 1, (17) = 0, (18) = 29, (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) = 0.10e-5, (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..2, {(1) = .0, (2) = .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..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 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}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .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..2, {(1) = 0, (2) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.0, (2) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = SA1(t), Y[2] = SA2(t)]`; YP[1] := -0.426118e-2*Y[1]; YP[2] := 0.452821e-1*Y[1]-0.641682e-1*Y[2]; 0 end proc, -1, 0, 0, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) if e = 1 then Y[1] := Y[1]+484 end if; return 0 end proc, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = SA1(t), Y[2] = SA2(t)]`; YP[1] := -0.426118e-2*Y[1]; YP[2] := 0.452821e-1*Y[1]-0.641682e-1*Y[2]; 0 end proc, -1, 0, 0, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) if e = 1 then Y[1] := Y[1]+484 end if; return 0 end proc, 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..2, {(1) = 0., (2) = 0.}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2)  ] ); _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, SA1(t), SA2(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

plots:-odeplot(Solution1,[[t,SA1(t),color="red",thickness=1],[t,SA2(t),color="blue",thickness=1]],t=0..50,labels=["Time t","Salicylic acid"],labeldirections=[horizontal,vertical]);

NULL

Download Event.mw

@vv Thanks, I just assumed that was the default order. As you undoubtedly know, the exact order isn't important, just that it is consistent for the left and right eigenvectors.

[Fixed and comment added.]

4 5 6 7 8 9 10 Last Page 6 of 63