Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 169 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

workbooks...

tomleslie 13876
May 02 2020
2 1

A .maple file extension represents a Maple "workbook"

A Maple workbook is a collection of files which are necessary for a given problem. It will "usually" comprise a Maple worksheet (ie a .mw file), maybe some input data files (eg .csv or .xlsx),  possibly some Maple procedures (written as .mpl files) and indeed any other files which might be required.

The capability to generate Maple workbooks was introduced in Maple 2016. I can think of no way that such .maple files could be read/executed/whatever in any Maple version earlier than Maple 2016 :-(

As the name suggests, Maple 2016 was released in 2016, so the .maple file format has been around/usable for five years. You state that you are using Maple 15,which was released in 2011: so you are using ten year-old software and you expect to be able to read/execute file formats which were developed five years after your was released?? Really?

Worth trying (since it will take about 10 secs) but pretty much guaranteed to fail.

  1. If you have the file "filename.maple", change the file extension to ".mw" and try opening it it the usual way
  2. change the file extension to ".mpl", then from within Maple use the command  read("filename.mpl")

 

Enter full path...

tomleslie 13876
May 02 2020
2 1

 to the Geogebra executable.

The attached shows success/failure for launhing my (more or less default) editor

> 

  restart:
> 

#
# This doesn't work
#
  system[launch]("notepad++.exe");

Error, could not execute command, notepad++.exe; The system cannot find the file specified.

  
> 

#
# This does
#
  system[launch]("C:/Program Files (x86)/Notepad++/notepad++.exe");

1248

 (1)
> 

 

Download launch.mw

Since you do not provide the procedure C...

tomleslie 13876
May 01 2020
2 1

I made up a simple "toy" example. This seems to work without any errors, but

  1. when using the equality constraint, add(d[i], i = 1 .. 7)  =220 the returned answer violates the constraints
  2. when using "equivalent" inequality constraints, [ add(d[i], i = 1 .. 7)  > 219, add(d[i], i = 1 .. 7)  < 221  ], then the returned answer obeys the constraints

See the attached. To make further progress, I think you will have to supply the procedure CHWP()

> 

  restart;
  CHWP:= proc(v::Vector)
              add( v^~3);
         end proc:
> 

#
# With inequality constraints
#
  result1 := DirectSearch:-Search( CHWP,
                                   [ add(d[i], i = 1 .. 7)  > 219,
                                     add(d[i], i = 1 .. 7)  < 221
                                   ],
                                   variables = [d = Vector(7)],
                                   assume =nonnegint,
                                   strategy = globalsearch,
                                   method = cdos,
                                   evaluationlimit = 10000,
                                   objectivetarget = 10^(-10)
                                 );
  add(result1[2]);

Warning, initial point [d[1] = .900000000000000, d[2] = .900000000000000, d[3] = .900000000000000, d[4] = .900000000000000, d[5] = .900000000000000, d[6] = .900000000000000, d[7] = .900000000000000] does not satisfy the inequality constraints; trying to find a feasible initial point

  

Warning, the new feasible initial point is [d[1] = 24, d[2] = 12, d[3] = 13, d[4] = 101, d[5] = 50, d[6] = 3, d[7] = 17]

  

[2.17468*10^5, Vector(7, {(1) = 31, (2) = 32, (3) = 31, (4) = 31, (5) = 31, (6) = 32, (7) = 32}), 2644]

  

220

 (1)
> 

#
# With equality constraints
#
  result2 := DirectSearch:-Search( CHWP,
                                   [ add(d[i], i = 1 .. 7)  =220],
                                   variables = [d = Vector(7)],
                                   assume =nonnegint,
                                   strategy = globalsearch,
                                   method = cdos,
                                   evaluationlimit = 10000,
                                   objectivetarget = 10^(-10)
                                 );
  add(result2[2]);

[1.29353*10^5, Vector(7, {(1) = 26, (2) = 27, (3) = 26, (4) = 27, (5) = 26, (6) = 26, (7) = 27}), 1428]

  

185

 (2)
> 

 


 

Download DS.mw

 

Something very weird...

tomleslie 13876
April 30 2020
0 3

When I downloaded/executed/read OP's worksheet the output on a few commands is just downright weird. The indexes/subscripts change from integers (normal) to floating point numbers!!! I have no idea why/how this is happening and since it is now past my bedtime, I won't be able to figure it out tonight.

Anyhow check the attached where I have suppressed output on all commands which don't exhibit the weird behaviour. Anyone know why the variable indexes/subscripts appear to be "floats" !??

I couldn't figure out what was casuing this weird subscript behaviour - so I decide to avoid it by changing the indexed variables, to equivalents with "inert" subscripts - eg a[1] becomes a__1, etc. I also "tidied" the OP's  original to (mainly to save typing), and the second attachment now achieves a solution for the problem (with several of the variables being arbitrary). In an overdetermined system I don't fiind this too surprising, and it is possible I screwed up somewhere - but like I said it's past my bedtime

> 

restart:
> 

omega := 485.0135000:
> 

odeA1:=0.0001995000000*a[1]*omega^2*cos(omega*t) + 0.0001995000000*a[3]*omega^2*cos(3*omega*t) + 0.004987500000*a[5]*omega^2*cos(5*omega*t) - 93.86*a[1]*cos(omega*t) - 10.42888889*a[3]*cos(3*omega*t) - 93.86*a[5]*cos(5*omega*t) + 1.633844046*10^18*(A[1]*cos(omega*t) + A[3]*cos(3*omega*t)/9 + A[5]*cos(5*omega*t))^2*(a[1]*cos(omega*t) + a[3]*cos(3*omega*t)/9 + a[5]*cos(5*omega*t)) + 0.02250237053*a[1]*omega*sin(omega*t) + 0.007500790177*a[3]*omega*sin(3*omega*t) + 0.1125118526*a[5]*omega*sin(5*omega*t):
> 

odeA2:=A[1]*omega*sin(omega*t) + A[3]*omega*sin(3*omega*t)/3 + 5*A[5]*omega*sin(5*omega*t) - 1.499250375*10^(-7)*(A[1]*cos(omega*t) + A[3]*cos(3*omega*t)/9 + A[5]*cos(5*omega*t))*(2.110712380*10^10 - 1.633844050*10^18*(a[1]*cos(omega*t) + a[3]*cos(3*omega*t)/9 + a[5]*cos(5*omega*t))^2) + 1.499250375*10^(-7):
> 

expandOdeA1:=expand(odeA1):
> 

combineOdeA1:=combine(odeA1):
> 

expandOdeA2:=expand(odeA2):
> 

combineOdeA2:=combine(odeA2):
> 

Eq1:=coeff(combineOdeA1, cos(omega*t)):
> 

Eq3:=coeff(combineOdeA2, cos(omega*t)):
> 

Eq5:=coeff(combineOdeA1, cos(3*omega*t)):
> 

Eq7:=coeff(combineOdeA2, cos(3*omega*t)):
> 

Eq9:=coeff(combineOdeA1, cos(5*omega*t)):
> 

Eq11:=coeff(combineOdeA2, cos(5*omega*t)):
> 

Eq13:=coeff(combineOdeA1, cos(7*omega*t)):
> 

Eq15:=coeff(combineOdeA2, cos(7*omega*t)):
> 

Eq17:=coeff(combineOdeA1, cos(9*omega*t)):
> 

Eq19:=coeff(combineOdeA2, cos(9*omega*t)):
> 

Eq111:=coeff(combineOdeA1, cos(11*omega*t)):
> 

Eq112:=coeff(combineOdeA2, cos(11*omega*t));

756031267.0*a[3]^2*A[5]+0.6123853260e11*a[5]^2*A[1]+0.1224770652e12*a[1]*a[5]*A[5]+1512062534.*a[3]*a[5]*A[3]

 (1)
> 

Eq131:=coeff(combineOdeA1, cos(13*omega*t)):
> 

Eq132:=coeff(combineOdeA2, cos(13*omega*t));

0.1360856280e11*a[3]*a[5]*A[5]+6804281401.*a[5]^2*A[3]

 (2)
> 

Eq151:=coeff(combineOdeA1, cos(15*omega*t)):
> 

Eq152:=coeff(combineOdeA2, cos(15*omega*t));

0.6123853260e11*a[5]^2*A[5]

 (3)

NULL

> 

#sys := { Eq1, Eq3, Eq5, Eq7, Eq9, Eq11, Eq13, Eq15, Eq17, Eq19, Eq111, Eq112, Eq131, Eq132, Eq151, Eq152}:
> 

#fsolve( sys, { a[1], a[3], a[5], A[1], A[3], A[5]});
> 

#resultx:=evalf(eval(x(t), %));
> 

 

``

Download weird.mw

NULL

> 

  restart:
> 

  odeA1:=   0.0001995000000*a__1*omega^2*cos(omega*t)
          + 0.0001995000000*a__3*omega^2*cos(3*omega*t)
          + 0.004987500000*a__5*omega^2*cos(5*omega*t)
          - 93.86*a__1*cos(omega*t)
          - 10.42888889*a__3*cos(3*omega*t)
          - 93.86*a__5*cos(5*omega*t)
          + 1.633844046*10^18*(A__1*cos(omega*t)
          + A__3*cos(3*omega*t)/9
          + A__5*cos(5*omega*t))^2*(a__1*cos(omega*t)
          + a__3*cos(3*omega*t)/9 + a__5*cos(5*omega*t))
          + 0.02250237053*a__1*omega*sin(omega*t)
          + 0.007500790177*a__3*omega*sin(3*omega*t)
          + 0.1125118526*a__5*omega*sin(5*omega*t);
  odeA2:=   A__1*omega*sin(omega*t)
          + A__3*omega*sin(3*omega*t)/3
          + 5*A__5*omega*sin(5*omega*t)
          - 1.499250375*10^(-7)*(A__1*cos(omega*t)
          + A__3*cos(3*omega*t)/9
          + A__5*cos(5*omega*t))*(2.110712380*10^10
          - 1.633844050*10^18*(a__1*cos(omega*t)
          + a__3*cos(3*omega*t)/9
          + a__5*cos(5*omega*t))^2)
          + 1.499250375*10^(-7);

0.1995000000e-3*a__1*omega^2*cos(omega*t)+0.1995000000e-3*a__3*omega^2*cos(3*omega*t)+0.4987500000e-2*a__5*omega^2*cos(5*omega*t)-93.86*a__1*cos(omega*t)-10.42888889*a__3*cos(3*omega*t)-93.86*a__5*cos(5*omega*t)+0.1633844046e19*(A__1*cos(omega*t)+(1/9)*A__3*cos(3*omega*t)+A__5*cos(5*omega*t))^2*(a__1*cos(omega*t)+(1/9)*a__3*cos(3*omega*t)+a__5*cos(5*omega*t))+0.2250237053e-1*a__1*omega*sin(omega*t)+0.7500790177e-2*a__3*omega*sin(3*omega*t)+.1125118526*a__5*omega*sin(5*omega*t)

  

A__1*omega*sin(omega*t)+(1/3)*A__3*omega*sin(3*omega*t)+5*A__5*omega*sin(5*omega*t)-0.1499250375e-6*(A__1*cos(omega*t)+(1/9)*A__3*cos(3*omega*t)+A__5*cos(5*omega*t))*(0.2110712380e11-0.1633844050e19*(a__1*cos(omega*t)+(1/9)*a__3*cos(3*omega*t)+a__5*cos(5*omega*t))^2)+0.1499250375e-6

 (1)
> 

  sol:= eval
        ( solve
          ( [ seq
              ( coeff~
                ( expand~
                  ( [odeA1,odeA2]
                  ),
                  cos(j*omega*t)
                )[],
                j=1..15,2
              )
            ],
            [ a__1, a__3, a__5, A__1, A__3, A__5]
          )[],
          omega=485.0135000
       );

[a__1 = -2.333333333*a__3+115.0000000*a__5, a__3 = a__3, a__5 = a__5, A__1 = A__1, A__3 = 2.999999999*A__1+15.*A__5, A__5 = A__5]

 (2)
> 

  seq( coeff~(eval([odeA1, odeA2], [sol[],omega=485.0135000] ), cos(j*omega*t)), j=1..15,2);

[0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0]

 (3)

``

Download oddEqs.mw

Observation...

tomleslie 13876
April 30 2020
1 0

Whilst the "international characters in file/path name" is a possibility, this is trivial to check (and you should do so!)

IMHO a more likely cause is that Windows10 is "losing" file associations - and making them almost impossible to renew/restore.

Entering something like "Problems with Windows 10 file associations" in your favourite search engine will provide a depressingly long list of entries. These seem to have come in two separate waves, associated with either Windows version 1803 and 1903 (or sometimes versions 1809 and 1909). The latest outbreak of the problem (associated with the 1903/1909 versions) seems to be on-going, with no definitive solution.

the "overbar"...

tomleslie 13876
April 28 2020
1 1

in Maple represents the complex conjugate, and there is a fundamental difference between

coeff(expr, x); and

coeff(expr, conjugate(x) )

In your worksheet, there are no terms in sin(theta[4]). There are however terms in sin( conjugate(theta[4])), and the coefficients of the latter can be collected.

See the attached

> 

  restart:
  temp1b := d[2]*conjugate(sin(theta[4])*a[5]*cos(theta[5])) + d[2]*conjugate(a[4]*sin(theta[4])) + d[3]*conjugate(sin(theta[4])*a[5]*cos(theta[5])) + d[3]*conjugate(a[4]*sin(theta[4]));
  temp1b := expand(temp1b);
#
# Try to collect terms in sin(theta[4]) - but
# sin(theta[4]) doesn't exist in the
# expression - so no "collection" will be
# performed
#
  temp1b := collect(temp1b, sin(theta[4]));
#
# Collect terms in sin(conjugate(theta[4])),
# This works
#
  temp1b := collect(temp1b, sin(conjugate(theta[4])));
#
# Try to get the coefficient of sin(theta[4]) -
# but this term does not exist in the expression,
# so its coefficient will be zero
#
  temp2b := coeff(temp1b, sin(theta[4]));
#
# Get the coefficient of sin(conjugate(theta[4])) -
# which does exist in the expression
#
  temp2b := coeff(temp1b, sin(conjugate(theta[4])));

d[2]*conjugate(sin(theta[4])*a[5]*cos(theta[5]))+d[2]*conjugate(a[4]*sin(theta[4]))+d[3]*conjugate(sin(theta[4])*a[5]*cos(theta[5]))+d[3]*conjugate(a[4]*sin(theta[4]))

  

d[2]*sin(conjugate(theta[4]))*conjugate(a[5])*cos(conjugate(theta[5]))+d[2]*conjugate(a[4])*sin(conjugate(theta[4]))+d[3]*sin(conjugate(theta[4]))*conjugate(a[5])*cos(conjugate(theta[5]))+d[3]*conjugate(a[4])*sin(conjugate(theta[4]))

  

d[2]*sin(conjugate(theta[4]))*conjugate(a[5])*cos(conjugate(theta[5]))+d[2]*conjugate(a[4])*sin(conjugate(theta[4]))+d[3]*sin(conjugate(theta[4]))*conjugate(a[5])*cos(conjugate(theta[5]))+d[3]*conjugate(a[4])*sin(conjugate(theta[4]))

  

(conjugate(a[5])*cos(conjugate(theta[5]))*d[2]+conjugate(a[5])*cos(conjugate(theta[5]))*d[3]+conjugate(a[4])*d[2]+conjugate(a[4])*d[3])*sin(conjugate(theta[4]))

  

0

  

conjugate(a[5])*cos(conjugate(theta[5]))*d[2]+conjugate(a[5])*cos(conjugate(theta[5]))*d[3]+conjugate(a[4])*d[2]+conjugate(a[4])*d[3]

 (1)

``

Download collConj.mw

 

The other half of the answer...

tomleslie 13876
April 26 2020
1 2

As far as I can tell you, are not going top get an "analytic" solution for your ODE problem.

On the other hand if I make up a couple of simple initial conditons, eg A(0)=0, D(A)(0)=1, then a numeric solution is possible - see the attached.

You can try the same approach with other ICs - assuming you know what these should be!

> 

  restart:
  q1:=0.4: q:=0.4: s1:=0: mu:=2: mu1:=0.7:
  mu3:=0.2: s:=0.5: M:=1.9109: mu2:=1+mu3-mu1:

  V:=unapply( M^(2)-M*(M^(2)+2*A)^(1/(2))+mu1*(M^(2)*mu+s1)-mu1*((M^(2)*mu)/(2))^(1/(2))*(((M^(2)*mu+3*s1-2*A)+sqrt((M^(2)*mu+3*s1-2*A)^(2)-12*M^(2)*mu*s1))^(1/(2))+(4*M^(2)*mu*s1)*((M^(2)*mu+3*s1-2*A)+sqrt((M^(2)*mu+3*s1-2*A)^(2)-12*M^(2)*mu*s1))^(-3/(2)))-(2*mu2)/((3*q-1))*((1-(q-1)*A)^((3*q-1)/(2*q-2))-1)-(2*mu3)/(s*(3*q1-1))*((1+s*(q1-1)*A)^((3*q1-1)/(2*q1-2))-1),A);

  ode:=diff(A(x),x,x)+V(A(x))=0;
  odeSol:=dsolve([ode, A(0)=0, D(A)(0)=1], numeric);
  plots:-odeplot( odeSol, [x, A(x)], x=0..3);

 

proc (A) options operator, arrow; 17.76369314-1.9109*(2*A+3.65153881)^(1/2)-.9458472435*2^(1/2)*(7.30307762-2*A+((-2*A+7.30307762)^2)^(1/2))^(1/2)-5.000000000/(1+.6*A)^.1666666667-4.000000000/(1-.30*A)^.1666666667 end proc

  

diff(diff(A(x), x), x)+17.76369314-1.9109*(2*A(x)+3.65153881)^(1/2)-.9458472435*2^(1/2)*(7.30307762-2*A(x)+((-2*A(x)+7.30307762)^2)^(1/2))^(1/2)-5.000000000/(1+.6*A(x))^.1666666667-4.000000000/(1-.30*A(x))^.1666666667 = 0

  

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 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 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) = 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}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.5047658755841546e-2, (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) = 1.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) = 1.0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 1.0, (2) = 0.28187194800466386e-8}, 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] = A(x), Y[2] = diff(A(x),x)]`; if -2*Y[1]+evalf(((-2*Y[1]+7.30307762)^2)^(1/2)) < -7.30307762 then YP[1] := undefined; return 0 end if; if 2*Y[1] < -3.65153881 then YP[1] := undefined; return 0 end if; if .6*Y[1] < -1 then YP[1] := undefined; return 0 end if; if -.30*Y[1] < -1 then YP[1] := undefined; return 0 end if; if (-2*Y[1]+7.30307762)^2 < 0 then YP[1] := undefined; return 0 end if; YP[2] := -17.76369314+1.9109*evalf((2*Y[1]+3.65153881)^(1/2))+1.33762999969091*evalf((7.30307762-2*Y[1]+evalf(((-2*Y[1]+7.30307762)^2)^(1/2)))^(1/2))+5.000000000/(1+.6*Y[1])^.1666666667+4.000000000/(1-.30*Y[1])^.1666666667; YP[1] := Y[2]; 0 end proc, -1, 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] = A(x), Y[2] = diff(A(x),x)]`; if -2*Y[1]+evalf(((-2*Y[1]+7.30307762)^2)^(1/2)) < -7.30307762 then YP[1] := undefined; return 0 end if; if 2*Y[1] < -3.65153881 then YP[1] := undefined; return 0 end if; if .6*Y[1] < -1 then YP[1] := undefined; return 0 end if; if -.30*Y[1] < -1 then YP[1] := undefined; return 0 end if; if (-2*Y[1]+7.30307762)^2 < 0 then YP[1] := undefined; return 0 end if; YP[2] := -17.76369314+1.9109*evalf((2*Y[1]+3.65153881)^(1/2))+1.33762999969091*evalf((7.30307762-2*Y[1]+evalf(((-2*Y[1]+7.30307762)^2)^(1/2)))^(1/2))+5.000000000/(1+.6*Y[1])^.1666666667+4.000000000/(1-.30*Y[1])^.1666666667; YP[1] := Y[2]; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..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 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) = [x, A(x), diff(A(x), x)], (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 odeSol.mw

 

Half an answer...

tomleslie 13876
April 26 2020
2 0

Consider the following diffrential equation:

 

(dA(x)/dx )^2+V(A)=0,
1- How to plot dA/dx vs A?

Can't do much with this unless the function V() is defined :-(

However for your second point, ie

2- How to find the coordinate of intersection point of f(x)=e-b*exp(x)-c*(1-2*x/c)^(1/2) with df/dx? e,b,c are constants.

see the attached - and in fact there seem to be two intersection points

> 

  restart;
  f:=x->e-b*exp(x)-c*(1-2*x/c)^(1/2);
#
# x-Value(s) of the intersection(s)
#
  sols:=[solve(f(x)=diff(f(x),x), x)];
#
# f(x) value(s) of the intersection(s)
#
  simplify(f(sols[1]));
  simplify(f(sols[2]));

proc (x) options operator, arrow; e-b*exp(x)-c*(1-2*x/c)^(1/2) end proc

  

[-(1/4)*e*(e+(e^2-4*c)^(1/2))/c+(1/2)*c+1/2, (1/4)*e*(-e+(e^2-4*c)^(1/2))/c+(1/2)*c+1/2]

  

e-b*exp(-(1/4)*(e*(e^2-4*c)^(1/2)-2*c^2+e^2-2*c)/c)-(1/2)*c*2^(1/2)*((e*(e^2-4*c)^(1/2)+e^2-2*c)/c^2)^(1/2)

  

e-b*exp((1/4)*(e*(e^2-4*c)^(1/2)+2*c^2-e^2+2*c)/c)-(1/2)*c*((-2*e*(e^2-4*c)^(1/2)+2*e^2-4*c)/c^2)^(1/2)

 (1)
> 

 


 

Download solEq.mw

Since you don't supply the worksheet...

tomleslie 13876
April 26 2020
0 0

In the attached have used a "toy" example, and determied the root using "regula falsi"

> 

  restart:
  with(Student[Precalculus]):
#
# define a yoy function whose root is desired
# and two -xValues which "bracket" the root
#
  f:=x->x*exp(x)-3:
  X__1:=0:X__2:=2:
#
# Specify the accuracy required
#
  tol:=1.0e-06:
#
# Iterate until tolerance is achieved
#
  for j from 1 do
    #
    # Get the x-intercept of the line between
    # the two points  given by {X__1, f(X__1)]
    # and [X__2, f(X__2)]
    #
      p:=evalf( [ Line
                  ( [X__1, f(X__1)],
                    [X__2, f(X__2)]
                  )
                ][4]
              );
    #
    # Exit if required tolerance has be achieved
    #
      if   abs(f(p))<tol
      then break
    #
    # Otherwise update the range end-points
    #
      elif f(p)>0
      then X__2:=p:
      else X__1:=p;
      fi;
  od:
  root=p, residual=f(p);

root = 1.049908802, residual = -0.544e-6

 (1)
> 

 

Download regulaFalsi.mw

Not accessing global variables?!...

tomleslie 13876
April 24 2020
1 1

I suggest you check exactly what is being passed to the procedure simplify2(), as shown in the last line of the attached.

Is that what you expected/wanted to pass as an argument?

> 

restart:
with(LinearAlgebra): with(ArrayTools): with(codegen): with(CodeGeneration): with(StringTools):
simplify2 := proc (Term)
  # Simplify a term. Additionally check if the simplification is advantageous.
  # In some cases the term gets longer after the simplification. Then discard.
  # Give the local procedure variable a long name. They must not be identical in Term.
  local simplify_tmpvar_c1, simplify_tmpvar_c1sum, simplify_tmpvar_c2, \
        simplify_tmpvar_c2sum, simplify_tmpvar_nms, simplify_tmpvar_k, Term2:
  # Assume all contents of the Term to be local variables to this procedure
  # Otherwise, variables like "l1" can be overwritten by occurences in the calling worksheet.
  simplify_tmpvar_nms:=convert(indets(Term,name),list): # get list of symbols
  for simplify_tmpvar_k from 1 to ColumnDimension(simplify_tmpvar_nms) do
    if simplify_tmpvar_nms[simplify_tmpvar_k] = Pi or simplify_tmpvar_nms[simplify_tmpvar_k] = 0 then
      next: # 0 and Pi can not be variables
    end if:
    eval(sprintf("local %s;", String(simplify_tmpvar_nms[simplify_tmpvar_k]))); # assume local
   end do:

  # Get computational effort before simplification
  simplify_tmpvar_c1:=add(cost~(Term)):
  simplify_tmpvar_c1sum := diff(simplify_tmpvar_c1,functions)+ \
                           diff(simplify_tmpvar_c1,additions)+ \
                           diff(simplify_tmpvar_c1,multiplications)+\
                           diff(simplify_tmpvar_c1,divisions):
  # Perform simplification. Attention: tries to use global variables to simplify the expression (see above)
  Term2 := simplify(Term):
  # Get effort after simplification
  simplify_tmpvar_c2:=add(cost~(Term2)):
  simplify_tmpvar_c2sum := diff(simplify_tmpvar_c2,functions)+\
                           diff(simplify_tmpvar_c2,additions)+\
                           diff(simplify_tmpvar_c2,multiplications)+\
                           diff(simplify_tmpvar_c2,divisions):
  if simplify_tmpvar_c2sum > simplify_tmpvar_c1sum then
    # The simplification was not successful. Take old version.
    Term2 := Term:
  end if:
  return Term2:
end proc:
T := m2*(l1+l2)*(qD2)^2; # dummy-expression for kinetic energy
l1 := length(T);
T_simpl := 'simplify2'(T);

m2*(l1+l2)*qD2^2

  

26

  

simplify2(m2*(26+l2)*qD2^2)

 (1)
> 

 


Download locGlob.mw

 

 

The attached...

tomleslie 13876
April 21 2020
1 2

shows how to import the data and "organise" it in a way whihc makes it reasonaly easy to access by year and category. Depending on exactly what you want to do with this data, this choice of organisation may not be optimal - but it is a start

> 

  restart;
> 

#
# Read the data - NB will have ti change the path
# name and file name to something appropriate for
# his/her installation
#
  someData:=ExcelTools:-Import("C:/Users/TomLeslie/Desktop/testDAta.xls"):
#
# There are many ways to "organise" this data into
# ways which are convenient to access - the following
# is just one possibility
#
  data:= table
         ( [ seq
             ( 2013+j=table
                      ( [ cases=someData[4..,6+4*j],
                          labCfd=someData[4..,7+4*j],
                          deaths=someData[4..,8+4*j],
                          CFR=someData[4..,9+4*j]
                        ]
                      ),
               j=0..4
             )
           ]
         ):
#
# Access data by year and category, where the former is
# one of 2013, 2014,....2017 and the latter is one of
# cases, labCfd, deaths, CFR
#
  data[2013][deaths];
  data[2017][cases];

_rtable[18446744074355802590]

  

_rtable[18446744074355804270]

 (1)
> 

 

> 

 

Download dataOrg.mw

 

 

You can use...

tomleslie 13876
April 21 2020
2 3

either pdsolve() or dsolve(), since the system doesn't really comprise partial differential equations. There are no deivatives with respect to the variabl 'T' and this variable doesn't appear explicitly anywhere.

Two possible solutions are shown in the attached

Download odeSol.mw

with plot3d...

tomleslie 13876
April 21 2020
1 2

it is generally better to "improve" a plot by changing th grid=[....] ooption whihc will increase the number of points at which function evaluation is performed - the default is 49*49

The attached code, produces the attached plot. I only use style=surface, becuase I (usually) prefer it to the default (wireframe)

imprPlot.mw

 

I agree with acer...

tomleslie 13876
April 20 2020
1 0

see the attached if you need evidence.

As well as defining 'e' in terms of 'e', you also use 'e' and 'e(t)' as well as 'i' , 'i[t]' and 'i(t)' in your equations - so ask yourself very carefully:- what is the relationship (if any) between these quantities

And the capability to inline files on this site is *STILL* broken

Download badNames.mw

patmatch() bug??...

tomleslie 13876
April 19 2020
3 0

Im guessing that patterns are "matched" using the patmatch() command.

The following 'toy' examples show a surprising (for me) result

patmatch( s^3, s^(y::realcons));
patmatch( 1/s^3, 1/s^(y::realcons));
patmatch( s^(-3), s^(y::realcons));

which returns

true
false ????????
true

 

First 77 78 79 80 81 82 83 Last Page 79 of 207