Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 169 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

I agree with dharr...

tomleslie 13876
March 08 2020
1 0

Fixing some (possibly non-printing?) strange characters, then there is nothing syntactically wrong with this code (see the attached)

as to the very basic solution - see the attached.

On the other hand, I have no idea what this code snippet is intended to achieve - or even whether your intention is achieved??

Why do I say this?

Well

  1. Although Maple provides an 'assign()' statement, in about ten years of writing Maple code, I can count on one hand the number of times I have ever used it. So forgive me if I suggest that it is *probably* unnecessary for your application
  2. Backquotes (ie ``) basically means treat the enclosed literal as a name. Forward (ie unevaluation) quotes (ie '') mean do not evaluate whatever is enclosed. So when you combine these as in '`something`', it is not syntactically incorrect, but do you really really understand what this will do? Again this is a construct which I have very, very rarely used - and only in exceptional circumstances
  3. From (1) and (2) above, I can only conclude that whatever you are trying to achieve, the code you supply (probably) isn't the way to do it!
> 

  restart;
  assignVariables:= proc( Class )
                          description "oppdatere globale variabler
                                       i 2D Math notasjon";
                          global `f__m,k`,    `f__t,0,k`,  `f__c,0,k`,
                                 `f__c,90,k`, `f__v,k`,    `f__r,k`,
                                 `E__0,mean`, `E__0,05`,   `E__90,mean`,
                                 `E__90,0,05`,`G__mean`,   `G__0,05`,
                                 `G__r,mean`, `G__r,0,05`, `rho__k`,
                                 `rho__mean`;
                          assign
                          ( '`f__m,k`',
                             Property( Class, "f_m,k" )     
                          );
                          assign
                          ( '`f__t,0,k`',
                            Property( Class, "f_t,0,k" )
                          ):
                          assign
                          ( '`f__t,90,k`',
                            Property( Class, "f_t,90,k" )
                          );
                          assign
                          ( '`f__c,0,k`',
                            Property( Class, "f_c,0,k" )
                          );
                          assign
                          ( '`f__c,90,k`',
                            Property( Class, "f_c,90,k" )
                          );
                          assign
                          ( '`f__v,k`',
                            Property( Class, "f_v,k" )
                          );
                          assign
                          ( '`f__r,k`',
                            Property( Class, "f_r,k" )
                          );
                          assign
                          ( '`E__0,mean`',
                             Property( Class, "E_0,mean" )
                          );
                          assign
                          ( '`E__0,05`',
                            Property( Class, "E_0,05" )
                          );
                          assign
                          ( '`E__90,mean`',
                            Property( Class, "E_90,mean" )
                          );
                          assign
                          ( '`E__90,0,05`',
                            Property( Class, "E_90,0,05" )
                          );
                          assign
                          ( '`G__mean`',
                            Property( Class, "G_mean" )
                          );
                          assign
                          ( '`G__0,05`',
                            Property( Class, "G_0,05" )
                          );
                          assign
                          ( '`G__r,mean`',
                            Property( Class, "G_r,mean" )
                          );
                          assign
                          ( '`G__r,0,05`',
                            Property( Class, "G_r,0,05" )
                          );
                          assign
                          ( '`rho__k`',
                            Property( Class, "rho_k" )
                          );
                          assign
                          ( '`rho__mean`',
                            Property( Class, "rho_mean" )
                          );
                   end proc:
> 

 

NULL

Download oddroc.mw

 

You could save a little typing...

tomleslie 13876
March 07 2020
0 2

with the attached

> 

  restart:
  with(DETools):
> 

  M1:= Matrix( 2, 1, [ diff(y(t),t),
                       diff(x(t),t)
                     ]
             ):
  M2:= Matrix(2,2, [ [  0,    1  ],
                     [-k/m, -c/m ]
                   ]
                 ):
  M3:= Matrix( 2, 1, [ y(t),
                       x(t)
                     ]
             ):
  odesys:= M1=M2.M3 ;

(Vector(2, {(1) = diff(y(t), t), (2) = diff(x(t), t)})) = (Vector(2, {(1) = x(t), (2) = -k*y(t)/m-c*x(t)/m}))

 (1)
> 

  f:= (cVal, mVal, kVal)-> DEplot
                           ( eval
                             ( [ seq
                                 ( rhs(odesys)[j,1]=lhs(odesys)[j,1],
                                   j=1..2
                                 )
                               ],
                               [ c=cVal, m=mVal, k=kVal ]
                             ),
                             [y(t), x(t)],
                             t=-10..10,
                             x=-10..10,
                             y=-10..10
                           ):
 f(0,1,1);
 f(0.25, 1,1);
 f(1.5,1,1);
 f(2.5,1,1);

 

 

 

 
> 

 

> 

 


 

Download odeplts.mw

 

Read the help!...

tomleslie 13876
March 07 2020
1 0

and working from what you provide, I come up with the attached

> 

  restart;
  with(GroupTheory):
  ct :=  < < e | p | q | r | s | t | u | v | w | x | y | z >,
           < p | q | e | y | u | w | z | r | x | t | v | s >,
           < q | e | p | v | z | x | s | y | t | w | r | u >,
           < r | z | t | s | e | y | v | x | p | u | q | w >,
           < s | w | y | e | r | q | x | u | z | v | t | p >,
           < t | r | z | x | w | u | e | q | y | p | s | v >,
           < u | x | v | p | y | e | t | z | s | r | w | q >,
           < v | u | x | z | q | r | y | w | e | s | p | t >,
           < w | y | s | t | x | z | p | e | v | q | u | r >,
           < x | v | u | w | t | s | q | p | r | e | z | y >,
           < y | s | w | u | p | v | r | t | q | z | e | x >,
           < z | t | r | q | v | p | w | s | u | y | x | e >
         >:
> 

  G2:= Group
       ( eval
         ( ct,
           [ seq
             ( ct[1,j]=j,
               j = 1..12
             )
           ]
         ),
         labels=[ seq
                  ( ct[1,j]=j,
                    j = 1..12
                  )
                ]
      );

_m1184145664

 (1)
> 

 


 

Download genGroup.mw

Interesting...

tomleslie 13876
March 05 2020
0 0

I can't make the spacecurve and the output of the DEplot3d() command appear in the same plot either - and at the moment I'm not sure why? I have to go now, but may get back to this issue later

As a "sanity" check, I generated the ode curves a slightly different way - and nnow these can be combined with the spacecurve. See the final execution group in the attached.

Now this is really weird

Within Maple the 'spacecurve' does not appear when combined with the DEplot3d() command, but when I upload the worksheet here, the 'missing' spacecurve 'magically appears. I have no idea what is going on!!!!!

However the final execution group using odeplot()+spacecurve always seems to be OK????

``

> 

restart
> 

NULL
> 

with(plottools)
> 

with(plots)
> 

with(DEtools)
> 

interface(displayprecision = 10)
> 

PDEtools:-declare(prime = t, quiet)
> 

``
> 

NULL
> 

l := 0
> 

ross_x := diff(x(t), t) = -(3*x(t)-sqrt(3)*y(t)*l)*(1-x(t)^2)
> 

ross_y := diff(y(t), t) = (-sqrt(3)*y(t)^2*x(t)*l)*(1/2)+(3/2)*y(t)*(2-z(t))*(z(t)-y(t)^2*sqrt(1-x(t)^2))/z(t)
> 

ross_z := diff(z(t), t) = (3*(1-z(t)))*(z(t)-y(t)^2*sqrt(1-x(t)^2))
> 

NULL
> 

rossler_sys := ross_x, ross_y, ross_z
> 

Find*fixed*points
> 

sol := solve({rhs(ross_x) = 0, rhs(ross_y) = 0, rhs(ross_z) = 0}, {x(t), y(t), z(t)})
> 

fp_sol := solve({rhs(ross_x) = 0, rhs(ross_y) = 0, rhs(ross_z) = 0}, {x(t), y(t), z(t)})
> 

fp1 := fp_sol[1]
> 

fp2 := fp_sol[2]
> 

fp3 := fp_sol[3]
> 

fp4 := fp_sol[4]
> 

fp5 := fp_sol[5]
> 

fp6 := fp_sol[6]
> 

fp7 := fp_sol[7]
> 

NULL
> 

Jacobian
> 

J := frontend(Student:-VectorCalculus:-Jacobian, [map(rhs, [rossler_sys]), ([x, y, z])(t)], [{`*`, `+`, list}, {}])
> 

ev1 := Student[LinearAlgebra][Eigenvalues](eval(J, fp1))
> 

ev2 := Student[LinearAlgebra][Eigenvalues](eval(J, fp2))
> 

ev3 := Student[LinearAlgebra][Eigenvalues](eval(J, fp3))
> 

ev4 := Student[LinearAlgebra][Eigenvalues](eval(J, fp4))
> 

ev5 := Student[LinearAlgebra][Eigenvalues](eval(J, fp5))
> 

ev6 := Student[LinearAlgebra][Eigenvalues](eval(J, fp6))
> 

ev7 := Student[LinearAlgebra][Eigenvalues](eval(J, fp7))
> 

 

> 

assign(params)
> 

NULL
> 

ics1 := [z(0) = .5, x(0) = .5, y(0) = .1]; ics2 := [z(0) = .5, x(0) = -.5, y(0) = .1]; ics3 := [z(0) = .5, x(0) = .3, y(0) = .4]; ics4 := [z(0) = .5, x(0) = -.3, y(0) = .4]; ics5 := [z(0) = .35, x(0) = .1, y(0) = .5]; ics6 := [z(0) = .35, x(0) = -.1, y(0) = .5]; ics7 := [z(0) = .5, x(0) = .5, y(0) = -.1]; ics8 := [z(0) = .5, x(0) = -.5, y(0) = -.1]; ics9 := [z(0) = .5, x(0) = .3, y(0) = -.4]; ics10 := [z(0) = .5, x(0) = -.3, y(0) = -.4]; ics11 := [z(0) = .35, x(0) = .1, y(0) = -.5]; ics12 := [z(0) = .35, x(0) = -.1, y(0) = -.5]
> 

``
> 

sys := {ross_x, ross_y, ross_z}
> 

NULL
> 

a1 := DEplot3d(sys, [z(t), y(t), x(t)], t = -10 .. 10, z = 0.1e-1 .. 1, y = -1 .. 1., x = -1 .. 1, [ics1, ics2, ics3, ics4, ics5, ics6, ics7, ics8, ics9, ics10, ics11, ics12], linecolor = [aquamarine, brown, navy, orange, sienna, gray, black, magenta, blue, yellow, green, coral], stepsize = .1, thickness = 2); a2 := plots:-spacecurve([0, y, y^2], y = -1 .. 1, color = red, thickness = 3); plots; display(a1, a2)

 
> 

#
# Same thing(?) a different way
#
  ics:= [ ics1, ics2, ics3, ics4, ics5, ics6,
          ics7, ics8, ics9, ics10, ics11, ics12
        ]:
  cols:= [ aquamarine, brown, navy, orange, sienna, gray,
           black, magenta, blue, yellow, green, coral
         ]:
  display( spacecurve
           ( [0, y, y^2],
             y = -1 .. 1,
             color = red,
             thickness = 3
           ),
           seq( odeplot
                ( dsolve
                  ( [sys[], ics[j][]],
                    numeric
                  ),
                  [z(t), y(t), x(t)],
                  t = -10 .. 10,
                  color = cols[j]
                ),
                j = 1 .. 12
              )
        );

 

NULL

> 

``


 

Download plotODE.mw

 

You have been mugged...

tomleslie 13876
March 04 2020
3 2

by Maple's 2D math input style, where spaces are (sometimes) interpreted as (implied) multiplications and (sometimes) not.

It can be almost impossible to simply 'see' the difference between functioning and non-functioning 2D input - although with 1D input it is pretty obvious.

See the attached


 

> 

restart
> 

with(Physics)

[`*`, `.`, Annihilation, AntiCommutator, Antisymmetrize, Assume, Bra, Bracket, Check, Christoffel, Coefficients, Commutator, CompactDisplay, Coordinates, Creation, D_, Dagger, Decompose, Define, Dgamma, Einstein, EnergyMomentum, Expand, ExteriorDerivative, Factor, FeynmanDiagrams, FeynmanIntegral, Fundiff, Geodesics, GrassmannParity, Gtaylor, Intc, Inverse, Ket, KillingVectors, KroneckerDelta, LeviCivita, Library, LieBracket, LieDerivative, Normal, NumericalRelativity, Parameters, PerformOnAnticommutativeSystem, Projector, Psigma, Redefine, Ricci, Riemann, Setup, Simplify, SortProducts, SpaceTimeVector, StandardModel, SubstituteTensor, SubstituteTensorIndices, SumOverRepeatedIndices, Symmetrize, TensorArray, Tetrads, ThreePlusOne, ToFieldComponents, ToSuperfields, Trace, TransformCoordinates, Vectors, Weyl, `^`, dAlembertian, d_, diff, g_, gamma_]

 (1)
> 

Setup(signature = `-+++`, coordinates = (X = [t, r, theta, phi]))

(`Defined Pauli sigma matrices (Psigma): `*sigma[0]*`, `*sigma[1]*`, `)*sigma[2]*`, `*sigma[3]

  

__________________________________________________

  

`Default differentiation variables for d_, D_ and dAlembertian are:`*{X = (t, r, theta, phi)}

  

`Systems of spacetime coordinates are:`*{X = (t, r, theta, phi)}

  

_______________________________________________________

  

[coordinatesystems = {X}, signature = `- + + +`]

 (2)
> 

Setup(g_ = -dt^2+dr^2+r^2*dtheta^2+r(sin(theta))^4*dphi^2)

[metric = {(1, 1) = -1, (2, 2) = 1, (3, 3) = r^2, (4, 4) = r(sin(theta))^4}]

 (3)
> 

Setup(g_ = -dt^2+dr^2+r^2*dtheta^2+sin(theta)^2*r^2*dphi^2)

[metric = {(1, 1) = -1, (2, 2) = 1, (3, 3) = r^2, (4, 4) = sin(theta)^2*r^2}]

 (4)
> 

Setup(g_ = -dt^2+dr^2+r^2*dtheta^2+sin(theta)^2*r^2*dphi^2)

[metric = {(1, 1) = -1, (2, 2) = 1, (3, 3) = r^2, (4, 4) = sin(theta)^2*r^2}]

 (5)
> 

Setup(g_ = -exp(2*nu(r, theta))*dt^2+(exp(2*psi(r, theta)))(dphi-omega(r, theta)*dt)^2+(exp(2*mu(r, theta)))(dtheta)^2+exp(2*lambda(r, theta))*dr^2)

Error, (in Physics:-Setup) invalid subscript selector

  
> 

#
# 1-D Input version of OP's original code. This will
# still fail
#
  Setup(g_ = -exp(2*nu(r, theta))*dt^2 + exp(2*psi(r, theta))(dphi - omega(r, theta)*dt)^2 + exp(2*mu(r, theta))(dtheta)^2 + exp(2*lambda(r, theta))*dr^2)

Error, (in Physics:-Setup) invalid subscript selector

  
> 

Setup(g_ = -exp(2*nu(r, theta))*dt^2+exp(2*psi(r, theta))*(dphi-omega(r, theta)*dt)^2+exp(2*mu(r, theta))*dtheta^2+exp(2*lambda(r, theta))*dr^2)

[metric = {(1, 1) = -exp(2*nu(r, theta))+exp(2*psi(r, theta))*omega(r, theta)^2, (1, 4) = -exp(2*psi(r, theta))*omega(r, theta), (2, 2) = exp(2*lambda(r, theta)), (3, 3) = exp(2*mu(r, theta)), (4, 4) = exp(2*psi(r, theta))}]

 (6)
> 

#
# 1-D Input version of fixed code, with inserted multiplication
# signs highlighted This will succeed
#
  Setup(g_ = -exp(2*nu(r, theta))*dt^2 + exp(2*psi(r, theta))*(dphi - omega(r, theta)*dt)^2 + exp(2*mu(r, theta))*dtheta^2 + exp(2*lambda(r, theta))*dr^2)

[metric = {(1, 1) = -exp(2*nu(r, theta))+exp(2*psi(r, theta))*omega(r, theta)^2, (1, 4) = -exp(2*psi(r, theta))*omega(r, theta), (2, 2) = exp(2*lambda(r, theta)), (3, 3) = exp(2*mu(r, theta)), (4, 4) = exp(2*psi(r, theta))}]

 (7)
> 

Setup(g_ = -exp(2*nu(r, theta))*dt^2+(dphi-omega(r, theta)*dt)^2*exp(2*psi(r, theta))+exp(2*mu(r, theta))*dtheta^2+exp(2*lambda(r, theta))*dr^2)

[metric = {(1, 1) = -exp(2*nu(r, theta))+exp(2*psi(r, theta))*omega(r, theta)^2, (1, 4) = -exp(2*psi(r, theta))*omega(r, theta), (2, 2) = exp(2*lambda(r, theta)), (3, 3) = exp(2*mu(r, theta)), (4, 4) = exp(2*psi(r, theta))}]

 (8)
> 

NULL


 

Download twoD_math.mw

 

A couple of ways...

tomleslie 13876
March 02 2020
0 0

to obtain your "desired" result are in the attached

> 

restart;
simplify(f*sqrt(g/f), symbolic);
evalb( simplify(f*sqrt(g/f), symbolic) = sqrt(f)*sqrt(g) );
testeq( simplify(f*sqrt(g/f), symbolic) = sqrt(f)*sqrt(g) );

f^(1/2)*g^(1/2)

  

true

  

true

 (1)
> 

restart;
simplify(f*sqrt(g/f), assume=positive);
evalb(  simplify(f*sqrt(g/f), assume=positive)= sqrt(f)*sqrt(g) );
testeq(  simplify(f*sqrt(g/f), assume=positive)= sqrt(f)*sqrt(g) )

g^(1/2)*f^(1/2)

  

true

  

true

 (2)
> 

 

Download simpl.mw

Like Carl...

tomleslie 13876
March 01 2020
1 0

I can't find anything syntactically correct.

In fact I was able to execute the complete worksheet (see attached) once I replaced the call to the GlobalOptimization:-GlobalSolve() command with the 'equivalent' command from free DirectSearch() package. This is not a comment on the relative merits of these two optimizers - I don't have the GlobalOptimization package because I can't justify the expense, and the DirectSearch() package is free.

The relevant DirectSearch command returns a set of values for k1,..k6, which is very similar to that embedded in your original woorksheet. Your original results below

[. 219132447080011505
 , [k1 = 0.852482740113834e-3, k2 = 0.52683998680924474e-4, k3 = 0., k4 = 0.5113239298267808e-1, k5 = 0.4363021255887466e-2, k6 = 0.]
]

Results from DirectSearch()

[  0.0526758778396827,
  Vector[column](1, [0.229512260761125]),
  Vector[column](6, [0.000847342169559454, 0.0000555396381708930, 0.5549173322*10^(-10), 0.0513531153470707, 0.0630964951648806, 0.208320218411464]), 4065]

The returned values of k1..k4 are very  similar for the two optimiser. The values ''k5' and 'k6' are noticeably different, particularly the latter, but the residual error from the two optimizers is very similar.

Unfortunately, the attached will not display inline here. I'm pretty sure this is because the worksheet uses an Array of plots. I have previously noticed that'plot arrays' seem to cause problems for this site.

Note also that you will not be able to re-execute this worksheet unless you have the DirectSearch() package installed

odeOpt.mw

Maybe(?)...

tomleslie 13876
February 29 2020
1 1

the attached will help?

Otherwise you are going to  have to get a lot more specific about your problem

> 

  restart;
#
# Define the ode system
#
  sys_ode := diff(x(t), t) = a*x(t)-b*x(t)*y(t)+2*sin(t),
             diff(y(t), t) = -c*y(t)+d*x(t)*y(t);
#
# Make up two boundary conditions, since OP does not
# supply these. Two boundary/initial conditions are
# essential in order to get a numerical solution -
# ie one without arbitrary integration constants
#
  inits:=x(0)=1, y(0)=1:
#
# Generate a solution procedure with [a, b, c, d]
# identified as parameters. NB this is only possible
# when the system is an initial value problem
#
  sol:=dsolve( [sys_ode, inits],
               parameters=[a,b,c,d],
               numeric
             );

diff(x(t), t) = a*x(t)-b*x(t)*y(t)+2*sin(t), diff(y(t), t) = -c*y(t)+d*x(t)*y(t)

  

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := [a = a, b = b, c = c, d = d]; _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) = 4, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .0, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..6, {(1) = 1., (2) = 1., (3) = Float(undefined), (4) = Float(undefined), (5) = Float(undefined), (6) = Float(undefined)})), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..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..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..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..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..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] = x(t), Y[2] = y(t)]`; YP[1] := Y[3]*Y[1]-Y[4]*Y[1]*Y[2]+2*sin(X); YP[2] := Y[1]*Y[2]*Y[6]-Y[2]*Y[5]; 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] = x(t), Y[2] = y(t)]`; YP[1] := Y[3]*Y[1]-Y[4]*Y[1]*Y[2]+2*sin(X); YP[2] := Y[1]*Y[2]*Y[6]-Y[2]*Y[5]; 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..6, {(1) = 0., (2) = 1., (3) = 1., (4) = undefined, (5) = undefined, (6) = undefined}); _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) = [t, x(t), y(t)], (4) = [a = a, b = b, c = c, d = d]}); _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

 (1)
> 

#
# Set values for the parameters and plot the solution
#
  sol(parameters= [1,2,3,4]);
  plots:-odeplot( sol,
                  [[t, x(t)], [t, y(t)]],
                  t=0..1
                );

[a = 1., b = 2., c = 3., d = 4.]

  

 
> 

#
# Set different values for the parameters and plot
# the solution
#
  sol(parameters= [4, 3, 2, 1]);
  plots:-odeplot( sol,
                  [ [t, x(t)], [t, y(t)] ],
                  t=0..1
                );

[a = 4., b = 3., c = 2., d = 1.]

  

 
> 

 


 

Download ODEpar.mw

 

AFAIK...

tomleslie 13876
February 29 2020
1 0

the only reason for using the 'MTM' package is that it allows users familiar with Matlab to enter commands in Maple using  'Matlab format' and receive the same results as (s)he would get when using Matlab.

Note that there is no need to have Matlab on your system: the MTM package does not actually  'call/execute' Matlab in any way, shape or form - all it does it replicate (within Maple) what Matlab would (might??) do.

Since Matlab is not involved in the execution of any of the commands in the MTM package, this means that every command in the MTM package is executed using Maple only. So you can think of the MTM package as a utility which

  1. translates a matlab-like command to equivalent Maple syntax
  2. executes the command within Maple
  3. translates the output of the Maple command at (2) above to be the same as that which would be returned by Matlab

Since all parsing/translation/execution is actually done within Maple, my view is that the MTM package is rather pointless. Just work out the equivalent Maple command for what you want to do, and use it.

There are reasons for using Maple-Matlab links (sometimes), but the only one I might be able to justify is the ability to use Maple as a 'symbolic toolbox' from within Matlab. Saves paying Mathworks for their add-on symbolic toolbox!

When...

tomleslie 13876
February 29 2020
0 0

you issue a 'restart' command halfway through your calculation this will remove all references to the DETools() package which you previously loaded using 'with(DETools)', so the DEplot() command is now unknown to subsequent commands in your worksheet.

Two trivial solutions

  1. Remove the 'restart' command from your worksheet  (it doesn't *seem* to be doing anything useful, anyway)
  2. if you feel you have to issue a 'restart' command then reload the 'DETools' package: ie make sure that you execute with(DETools) after the 'restart'

Another way...

tomleslie 13876
February 26 2020
2 1

Having tried, it seems as if the triangle-based Koch fractal isn't actually a Lindemayer system: But I might be wrong!!

However a variant with 90 degree angles (see https://en.wikipedia.org/wiki/L-system#Example_4:_Koch_curve) is a Lindemayer system, and this is the option which Maple implements:-(

With a bit of thought one can actually use the LindeMayer system built into the Fractals:-LSystem() package ( together with a bit of rotation/translation/scaling) to produce the desired curve(s). See the attached

> 

  restart:
  with(plots):
  doKoch:= proc( ord::nonnegint )
                 uses Fractals:-LSystem, plottools, plots:
                 local state:= "F",
                       cons:= [ "F"="draw:1",
                                "+"="turn:-60",
                                "-"="turn:-240"
                              ],
                       rules:= ["F"="F+F-F+F"],
                       newstate1:= Iterate(state, rules, ord),
                       p1:= LSystemPlot
                            ( newstate1, cons ),
                       sf:= getdata(p1)[3][-1,1],
                       p2:= translate
                            ( rotate
                              ( p1, 2*Pi/3 ),
                              sf,
                              0
                            ),
                       p3:= rotate
                            ( reflect
                              ( p1,
                                [ [0,0],
                                  [sf,0]
                                ],
                                [0,0]
                              ),
                              Pi/3
                            ):
                 return display( scale( p1, 1/sf, 1/sf, [0,0] ),
                                 scale( p2, 1/sf, 1/sf, [0,0] ),
                                 scale( p3, 1/sf, 1/sf, [0,0] ),
                                 color=red
                               ):
         end proc:
> 

display( [seq( doKoch(j), j=0..4)], insequence);

 
> 

 

Download Koch.mw

How about...

tomleslie 13876
February 25 2020
1 1

the attached


 

> 

  restart:
  with(plottools):with(plots):with(VectorCalculus):
#
# Define the triangle
#
  t:=[[0,0], [1,0], [0,1]]:
#
# Plot the above triangle
#
  display(polygon(t, color=red), axes=boxed);
#
# Integrate over the above triangle
#
  int( x*y, [x,y] = Triangle(convert~(t, Vector)[])  );

 

1/24

 (1)
> 

 


 

Download tri.mw

What exactly is the problem?...

tomleslie 13876
February 24 2020
1 2

In order to obtain a numeric answer, values will have to be supplied for 'T__D' and 'T' - as in the attached.

It is also rather suspiciious that the output expression contains terms such as ln(e)5 - which really ought to evluate to 5. When you originally defined the integral, did you use ex (with a literal character 'e') rather than exp(x)??

This last point would be easy for me to check if you had uploaded a worksheet using the big green up-arrow in the MApleprimes toolbar

Anyhow the attached may help

> 

  restart:
> 

#
# Evaluate the integral under the assumption that the
# quoted upper limit (T__D/T) is greater than the the
# lower limit (1.0e-05)
#
  I1:= 0.002+(11500/T__D)*(T/T__D)^5*int( x^5/(exp(x)-1)/(1-exp(-x)), x=1.0e-05..T__D/T) assuming T__D/T > 1.0e-05;
#
# Obtaina a numerical result for specific values of
# T__D and T (which comply with the above assumption
#
  ans:= eval(I1, [T__D=1, T=0.1]);

0.2e-2-11500*(-124.4200000*exp(T__D/T)*T^5+124.4200000*T^5+(47.12388981*I)*exp(T__D/T)*T__D^4*T-32.89868134*exp(T__D/T)*T__D^3*T^2-(47.12388981*I)*T__D^4*T+120.*polylog(5, exp(T__D/T))*exp(T__D/T)*T^5+60.*polylog(3, exp(T__D/T))*exp(T__D/T)*T__D^2*T^3-120.*polylog(4, exp(T__D/T))*exp(T__D/T)*T__D*T^4+20.*exp(T__D/T)*dilog(exp(T__D/T))*T__D^3*T^2+15.*exp(T__D/T)*ln(exp(T__D/T)-1)*T__D^4*T+32.89868134*T__D^3*T^2+exp(T__D/T)*T__D^5-120.*polylog(5, exp(T__D/T))*T^5-60.*polylog(3, exp(T__D/T))*T__D^2*T^3+120.*polylog(4, exp(T__D/T))*T__D*T^4-20.*dilog(exp(T__D/T))*T__D^3*T^2-15.*ln(exp(T__D/T)-1)*T__D^4*T)/(T__D^6*(exp(T__D/T)-1.))

  

13.38439722+0.5796503975e-4*I

 (1)
> 

 


 

Download EI2.mw

Hmmm...

tomleslie 13876
February 22 2020
0 5

Why did you start a new question?

Just don't!

Use the 'reply' and 'answer' buttons in Mapleprimes to avoid duplicating your questions

The file you uploaded does not actually define any boundary conditions: check the contents of the name `BCS` in the worksheet yu have posted. A boundary condition is an equation.

If I take a guess at what you intended then I get the fina executon group in the attached whihc solves the ODE and plots solutions for all dependent functions, at least untile the value of the independent variable reaches ~0.537 at whihc the solution "explodes". Possibly a typo in the definition of the ODE system, or maybe your ODE system definition is correct and you expect(?) to generate infinite values for dependent functions?

Either way, check the final execution group in the attached

NULL

restart:

diff(diff(u(x, y), y), y) = -(diff(diff(u(x, y), x), x))+u(x, y)^2

50

diff(u(x, y), x)

u(x, y)-1

u(x, y)

diff(u(x, y), y)

1

1

10

(1/2)*(-u[m+2](y)-3*u[m](y)+4*u[m+1](y))/h

(1/2)*(u[m-2](y)+3*u[m](y)-4*u[m-1](y))/h

(u[m+1](y)-u[m-1](y))/h^2

(u[m-1](y)-2*u[m](y)+u[m+1](y))/h^2

(1/2)*(-u[2](y)-3*u[0](y)+4*u[1](y))/h

u[11](y)-1

(1/2)*(-u[2](y)-3*u[0](y)+4*u[1](y))/h

u[11](y)-1

-(1/3)*u[2](y)+(4/3)*u[1](y)

1

u[1](0), u[2](0), u[3](0), u[4](0), u[5](0), u[6](0), u[7](0), u[8](0), u[9](0), u[10](0)

0, 0, 0, 0, 0, 0, 0, 0, 0, 0

u[1](0), u[2](0), u[3](0), u[4](0), u[5](0), u[6](0), u[7](0), u[8](0), u[9](0), u[10](0), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0

1/11

Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations

> 

BCSS:=zip(`=`, [BC3], [BC4])[]:
sol:=dsolve({BCSS,seq(eq[i],i=1..N)}, numeric);
plots:-odeplot(sol, [ seq( [y, u[j](y)], j=1..10)], y=0..0.2, axes =boxed);

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) = 20, (2) = 20, (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.4171618806480616e-4, (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..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .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..20, {(1) = .1, (2) = .1, (3) = .1, (4) = .1, (5) = .1, (6) = .1, (7) = .1, (8) = .1, (9) = .1, (10) = .1, (11) = .1, (12) = .1, (13) = .1, (14) = .1, (15) = .1, (16) = .1, (17) = .1, (18) = .1, (19) = .1, (20) = .1}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, 1..20, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (1, 8) = .0, (1, 9) = .0, (1, 10) = .0, (1, 11) = .0, (1, 12) = .0, (1, 13) = .0, (1, 14) = .0, (1, 15) = .0, (1, 16) = .0, (1, 17) = .0, (1, 18) = .0, (1, 19) = .0, (1, 20) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = .0, (2, 9) = .0, (2, 10) = .0, (2, 11) = .0, (2, 12) = .0, (2, 13) = .0, (2, 14) = .0, (2, 15) = .0, (2, 16) = .0, (2, 17) = .0, (2, 18) = .0, (2, 19) = .0, (2, 20) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (3, 8) = .0, (3, 9) = .0, (3, 10) = .0, (3, 11) = .0, (3, 12) = .0, (3, 13) = .0, (3, 14) = .0, (3, 15) = .0, (3, 16) = .0, (3, 17) = .0, (3, 18) = .0, (3, 19) = .0, (3, 20) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (4, 8) = .0, (4, 9) = .0, (4, 10) = .0, (4, 11) = .0, (4, 12) = .0, (4, 13) = .0, (4, 14) = .0, (4, 15) = .0, (4, 16) = .0, (4, 17) = .0, (4, 18) = .0, (4, 19) = .0, (4, 20) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (5, 8) = .0, (5, 9) = .0, (5, 10) = .0, (5, 11) = .0, (5, 12) = .0, (5, 13) = .0, (5, 14) = .0, (5, 15) = .0, (5, 16) = .0, (5, 17) = .0, (5, 18) = .0, (5, 19) = .0, (5, 20) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (6, 8) = .0, (6, 9) = .0, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0, (6, 13) = .0, (6, 14) = .0, (6, 15) = .0, (6, 16) = .0, (6, 17) = .0, (6, 18) = .0, (6, 19) = .0, (6, 20) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (7, 8) = .0, (7, 9) = .0, (7, 10) = .0, (7, 11) = .0, (7, 12) = .0, (7, 13) = .0, (7, 14) = .0, (7, 15) = .0, (7, 16) = .0, (7, 17) = .0, (7, 18) = .0, (7, 19) = .0, (7, 20) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (8, 8) = .0, (8, 9) = .0, (8, 10) = .0, (8, 11) = .0, (8, 12) = .0, (8, 13) = .0, (8, 14) = .0, (8, 15) = .0, (8, 16) = .0, (8, 17) = .0, (8, 18) = .0, (8, 19) = .0, (8, 20) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (9, 8) = .0, (9, 9) = .0, (9, 10) = .0, (9, 11) = .0, (9, 12) = .0, (9, 13) = .0, (9, 14) = .0, (9, 15) = .0, (9, 16) = .0, (9, 17) = .0, (9, 18) = .0, (9, 19) = .0, (9, 20) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (10, 8) = .0, (10, 9) = .0, (10, 10) = .0, (10, 11) = .0, (10, 12) = .0, (10, 13) = .0, (10, 14) = .0, (10, 15) = .0, (10, 16) = .0, (10, 17) = .0, (10, 18) = .0, (10, 19) = .0, (10, 20) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (11, 8) = .0, (11, 9) = .0, (11, 10) = .0, (11, 11) = .0, (11, 12) = .0, (11, 13) = .0, (11, 14) = .0, (11, 15) = .0, (11, 16) = .0, (11, 17) = .0, (11, 18) = .0, (11, 19) = .0, (11, 20) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (12, 8) = .0, (12, 9) = .0, (12, 10) = .0, (12, 11) = .0, (12, 12) = .0, (12, 13) = .0, (12, 14) = .0, (12, 15) = .0, (12, 16) = .0, (12, 17) = .0, (12, 18) = .0, (12, 19) = .0, (12, 20) = .0, (13, 1) = .0, (13, 2) = .0, (13, 3) = .0, (13, 4) = .0, (13, 5) = .0, (13, 6) = .0, (13, 7) = .0, (13, 8) = .0, (13, 9) = .0, (13, 10) = .0, (13, 11) = .0, (13, 12) = .0, (13, 13) = .0, (13, 14) = .0, (13, 15) = .0, (13, 16) = .0, (13, 17) = .0, (13, 18) = .0, (13, 19) = .0, (13, 20) = .0, (14, 1) = .0, (14, 2) = .0, (14, 3) = .0, (14, 4) = .0, (14, 5) = .0, (14, 6) = .0, (14, 7) = .0, (14, 8) = .0, (14, 9) = .0, (14, 10) = .0, (14, 11) = .0, (14, 12) = .0, (14, 13) = .0, (14, 14) = .0, (14, 15) = .0, (14, 16) = .0, (14, 17) = .0, (14, 18) = .0, (14, 19) = .0, (14, 20) = .0, (15, 1) = .0, (15, 2) = .0, (15, 3) = .0, (15, 4) = .0, (15, 5) = .0, (15, 6) = .0, (15, 7) = .0, (15, 8) = .0, (15, 9) = .0, (15, 10) = .0, (15, 11) = .0, (15, 12) = .0, (15, 13) = .0, (15, 14) = .0, (15, 15) = .0, (15, 16) = .0, (15, 17) = .0, (15, 18) = .0, (15, 19) = .0, (15, 20) = .0, (16, 1) = .0, (16, 2) = .0, (16, 3) = .0, (16, 4) = .0, (16, 5) = .0, (16, 6) = .0, (16, 7) = .0, (16, 8) = .0, (16, 9) = .0, (16, 10) = .0, (16, 11) = .0, (16, 12) = .0, (16, 13) = .0, (16, 14) = .0, (16, 15) = .0, (16, 16) = .0, (16, 17) = .0, (16, 18) = .0, (16, 19) = .0, (16, 20) = .0, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = .0, (17, 5) = .0, (17, 6) = .0, (17, 7) = .0, (17, 8) = .0, (17, 9) = .0, (17, 10) = .0, (17, 11) = .0, (17, 12) = .0, (17, 13) = .0, (17, 14) = .0, (17, 15) = .0, (17, 16) = .0, (17, 17) = .0, (17, 18) = .0, (17, 19) = .0, (17, 20) = .0, (18, 1) = .0, (18, 2) = .0, (18, 3) = .0, (18, 4) = .0, (18, 5) = .0, (18, 6) = .0, (18, 7) = .0, (18, 8) = .0, (18, 9) = .0, (18, 10) = .0, (18, 11) = .0, (18, 12) = .0, (18, 13) = .0, (18, 14) = .0, (18, 15) = .0, (18, 16) = .0, (18, 17) = .0, (18, 18) = .0, (18, 19) = .0, (18, 20) = .0, (19, 1) = .0, (19, 2) = .0, (19, 3) = .0, (19, 4) = .0, (19, 5) = .0, (19, 6) = .0, (19, 7) = .0, (19, 8) = .0, (19, 9) = .0, (19, 10) = .0, (19, 11) = .0, (19, 12) = .0, (19, 13) = .0, (19, 14) = .0, (19, 15) = .0, (19, 16) = .0, (19, 17) = .0, (19, 18) = .0, (19, 19) = .0, (19, 20) = .0, (20, 1) = .0, (20, 2) = .0, (20, 3) = .0, (20, 4) = .0, (20, 5) = .0, (20, 6) = .0, (20, 7) = .0, (20, 8) = .0, (20, 9) = .0, (20, 10) = .0, (20, 11) = .0, (20, 12) = .0, (20, 13) = .0, (20, 14) = .0, (20, 15) = .0, (20, 16) = .0, (20, 17) = .0, (20, 18) = .0, (20, 19) = .0, (20, 20) = .0}, datatype = float[8], order = C_order), Array(1..20, 1..20, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (1, 8) = .0, (1, 9) = .0, (1, 10) = .0, (1, 11) = .0, (1, 12) = .0, (1, 13) = .0, (1, 14) = .0, (1, 15) = .0, (1, 16) = .0, (1, 17) = .0, (1, 18) = .0, (1, 19) = .0, (1, 20) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = .0, (2, 9) = .0, (2, 10) = .0, (2, 11) = .0, (2, 12) = .0, (2, 13) = .0, (2, 14) = .0, (2, 15) = .0, (2, 16) = .0, (2, 17) = .0, (2, 18) = .0, (2, 19) = .0, (2, 20) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (3, 8) = .0, (3, 9) = .0, (3, 10) = .0, (3, 11) = .0, (3, 12) = .0, (3, 13) = .0, (3, 14) = .0, (3, 15) = .0, (3, 16) = .0, (3, 17) = .0, (3, 18) = .0, (3, 19) = .0, (3, 20) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (4, 8) = .0, (4, 9) = .0, (4, 10) = .0, (4, 11) = .0, (4, 12) = .0, (4, 13) = .0, (4, 14) = .0, (4, 15) = .0, (4, 16) = .0, (4, 17) = .0, (4, 18) = .0, (4, 19) = .0, (4, 20) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (5, 8) = .0, (5, 9) = .0, (5, 10) = .0, (5, 11) = .0, (5, 12) = .0, (5, 13) = .0, (5, 14) = .0, (5, 15) = .0, (5, 16) = .0, (5, 17) = .0, (5, 18) = .0, (5, 19) = .0, (5, 20) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (6, 8) = .0, (6, 9) = .0, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0, (6, 13) = .0, (6, 14) = .0, (6, 15) = .0, (6, 16) = .0, (6, 17) = .0, (6, 18) = .0, (6, 19) = .0, (6, 20) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (7, 8) = .0, (7, 9) = .0, (7, 10) = .0, (7, 11) = .0, (7, 12) = .0, (7, 13) = .0, (7, 14) = .0, (7, 15) = .0, (7, 16) = .0, (7, 17) = .0, (7, 18) = .0, (7, 19) = .0, (7, 20) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (8, 8) = .0, (8, 9) = .0, (8, 10) = .0, (8, 11) = .0, (8, 12) = .0, (8, 13) = .0, (8, 14) = .0, (8, 15) = .0, (8, 16) = .0, (8, 17) = .0, (8, 18) = .0, (8, 19) = .0, (8, 20) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (9, 8) = .0, (9, 9) = .0, (9, 10) = .0, (9, 11) = .0, (9, 12) = .0, (9, 13) = .0, (9, 14) = .0, (9, 15) = .0, (9, 16) = .0, (9, 17) = .0, (9, 18) = .0, (9, 19) = .0, (9, 20) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (10, 8) = .0, (10, 9) = .0, (10, 10) = .0, (10, 11) = .0, (10, 12) = .0, (10, 13) = .0, (10, 14) = .0, (10, 15) = .0, (10, 16) = .0, (10, 17) = .0, (10, 18) = .0, (10, 19) = .0, (10, 20) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (11, 8) = .0, (11, 9) = .0, (11, 10) = .0, (11, 11) = .0, (11, 12) = .0, (11, 13) = .0, (11, 14) = .0, (11, 15) = .0, (11, 16) = .0, (11, 17) = .0, (11, 18) = .0, (11, 19) = .0, (11, 20) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (12, 8) = .0, (12, 9) = .0, (12, 10) = .0, (12, 11) = .0, (12, 12) = .0, (12, 13) = .0, (12, 14) = .0, (12, 15) = .0, (12, 16) = .0, (12, 17) = .0, (12, 18) = .0, (12, 19) = .0, (12, 20) = .0, (13, 1) = .0, (13, 2) = .0, (13, 3) = .0, (13, 4) = .0, (13, 5) = .0, (13, 6) = .0, (13, 7) = .0, (13, 8) = .0, (13, 9) = .0, (13, 10) = .0, (13, 11) = .0, (13, 12) = .0, (13, 13) = .0, (13, 14) = .0, (13, 15) = .0, (13, 16) = .0, (13, 17) = .0, (13, 18) = .0, (13, 19) = .0, (13, 20) = .0, (14, 1) = .0, (14, 2) = .0, (14, 3) = .0, (14, 4) = .0, (14, 5) = .0, (14, 6) = .0, (14, 7) = .0, (14, 8) = .0, (14, 9) = .0, (14, 10) = .0, (14, 11) = .0, (14, 12) = .0, (14, 13) = .0, (14, 14) = .0, (14, 15) = .0, (14, 16) = .0, (14, 17) = .0, (14, 18) = .0, (14, 19) = .0, (14, 20) = .0, (15, 1) = .0, (15, 2) = .0, (15, 3) = .0, (15, 4) = .0, (15, 5) = .0, (15, 6) = .0, (15, 7) = .0, (15, 8) = .0, (15, 9) = .0, (15, 10) = .0, (15, 11) = .0, (15, 12) = .0, (15, 13) = .0, (15, 14) = .0, (15, 15) = .0, (15, 16) = .0, (15, 17) = .0, (15, 18) = .0, (15, 19) = .0, (15, 20) = .0, (16, 1) = .0, (16, 2) = .0, (16, 3) = .0, (16, 4) = .0, (16, 5) = .0, (16, 6) = .0, (16, 7) = .0, (16, 8) = .0, (16, 9) = .0, (16, 10) = .0, (16, 11) = .0, (16, 12) = .0, (16, 13) = .0, (16, 14) = .0, (16, 15) = .0, (16, 16) = .0, (16, 17) = .0, (16, 18) = .0, (16, 19) = .0, (16, 20) = .0, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = .0, (17, 5) = .0, (17, 6) = .0, (17, 7) = .0, (17, 8) = .0, (17, 9) = .0, (17, 10) = .0, (17, 11) = .0, (17, 12) = .0, (17, 13) = .0, (17, 14) = .0, (17, 15) = .0, (17, 16) = .0, (17, 17) = .0, (17, 18) = .0, (17, 19) = .0, (17, 20) = .0, (18, 1) = .0, (18, 2) = .0, (18, 3) = .0, (18, 4) = .0, (18, 5) = .0, (18, 6) = .0, (18, 7) = .0, (18, 8) = .0, (18, 9) = .0, (18, 10) = .0, (18, 11) = .0, (18, 12) = .0, (18, 13) = .0, (18, 14) = .0, (18, 15) = .0, (18, 16) = .0, (18, 17) = .0, (18, 18) = .0, (18, 19) = .0, (18, 20) = .0, (19, 1) = .0, (19, 2) = .0, (19, 3) = .0, (19, 4) = .0, (19, 5) = .0, (19, 6) = .0, (19, 7) = .0, (19, 8) = .0, (19, 9) = .0, (19, 10) = .0, (19, 11) = .0, (19, 12) = .0, (19, 13) = .0, (19, 14) = .0, (19, 15) = .0, (19, 16) = .0, (19, 17) = .0, (19, 18) = .0, (19, 19) = .0, (19, 20) = .0, (20, 1) = .0, (20, 2) = .0, (20, 3) = .0, (20, 4) = .0, (20, 5) = .0, (20, 6) = .0, (20, 7) = .0, (20, 8) = .0, (20, 9) = .0, (20, 10) = .0, (20, 11) = .0, (20, 12) = .0, (20, 13) = .0, (20, 14) = .0, (20, 15) = .0, (20, 16) = .0, (20, 17) = .0, (20, 18) = .0, (20, 19) = .0, (20, 20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, 1..20, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (1, 8) = .0, (1, 9) = .0, (1, 10) = .0, (1, 11) = .0, (1, 12) = .0, (1, 13) = .0, (1, 14) = .0, (1, 15) = .0, (1, 16) = .0, (1, 17) = .0, (1, 18) = .0, (1, 19) = .0, (1, 20) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = .0, (2, 9) = .0, (2, 10) = .0, (2, 11) = .0, (2, 12) = .0, (2, 13) = .0, (2, 14) = .0, (2, 15) = .0, (2, 16) = .0, (2, 17) = .0, (2, 18) = .0, (2, 19) = .0, (2, 20) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (3, 8) = .0, (3, 9) = .0, (3, 10) = .0, (3, 11) = .0, (3, 12) = .0, (3, 13) = .0, (3, 14) = .0, (3, 15) = .0, (3, 16) = .0, (3, 17) = .0, (3, 18) = .0, (3, 19) = .0, (3, 20) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (4, 8) = .0, (4, 9) = .0, (4, 10) = .0, (4, 11) = .0, (4, 12) = .0, (4, 13) = .0, (4, 14) = .0, (4, 15) = .0, (4, 16) = .0, (4, 17) = .0, (4, 18) = .0, (4, 19) = .0, (4, 20) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (5, 8) = .0, (5, 9) = .0, (5, 10) = .0, (5, 11) = .0, (5, 12) = .0, (5, 13) = .0, (5, 14) = .0, (5, 15) = .0, (5, 16) = .0, (5, 17) = .0, (5, 18) = .0, (5, 19) = .0, (5, 20) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (6, 8) = .0, (6, 9) = .0, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0, (6, 13) = .0, (6, 14) = .0, (6, 15) = .0, (6, 16) = .0, (6, 17) = .0, (6, 18) = .0, (6, 19) = .0, (6, 20) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (7, 8) = .0, (7, 9) = .0, (7, 10) = .0, (7, 11) = .0, (7, 12) = .0, (7, 13) = .0, (7, 14) = .0, (7, 15) = .0, (7, 16) = .0, (7, 17) = .0, (7, 18) = .0, (7, 19) = .0, (7, 20) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (8, 8) = .0, (8, 9) = .0, (8, 10) = .0, (8, 11) = .0, (8, 12) = .0, (8, 13) = .0, (8, 14) = .0, (8, 15) = .0, (8, 16) = .0, (8, 17) = .0, (8, 18) = .0, (8, 19) = .0, (8, 20) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (9, 8) = .0, (9, 9) = .0, (9, 10) = .0, (9, 11) = .0, (9, 12) = .0, (9, 13) = .0, (9, 14) = .0, (9, 15) = .0, (9, 16) = .0, (9, 17) = .0, (9, 18) = .0, (9, 19) = .0, (9, 20) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (10, 8) = .0, (10, 9) = .0, (10, 10) = .0, (10, 11) = .0, (10, 12) = .0, (10, 13) = .0, (10, 14) = .0, (10, 15) = .0, (10, 16) = .0, (10, 17) = .0, (10, 18) = .0, (10, 19) = .0, (10, 20) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (11, 8) = .0, (11, 9) = .0, (11, 10) = .0, (11, 11) = .0, (11, 12) = .0, (11, 13) = .0, (11, 14) = .0, (11, 15) = .0, (11, 16) = .0, (11, 17) = .0, (11, 18) = .0, (11, 19) = .0, (11, 20) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (12, 8) = .0, (12, 9) = .0, (12, 10) = .0, (12, 11) = .0, (12, 12) = .0, (12, 13) = .0, (12, 14) = .0, (12, 15) = .0, (12, 16) = .0, (12, 17) = .0, (12, 18) = .0, (12, 19) = .0, (12, 20) = .0, (13, 1) = .0, (13, 2) = .0, (13, 3) = .0, (13, 4) = .0, (13, 5) = .0, (13, 6) = .0, (13, 7) = .0, (13, 8) = .0, (13, 9) = .0, (13, 10) = .0, (13, 11) = .0, (13, 12) = .0, (13, 13) = .0, (13, 14) = .0, (13, 15) = .0, (13, 16) = .0, (13, 17) = .0, (13, 18) = .0, (13, 19) = .0, (13, 20) = .0, (14, 1) = .0, (14, 2) = .0, (14, 3) = .0, (14, 4) = .0, (14, 5) = .0, (14, 6) = .0, (14, 7) = .0, (14, 8) = .0, (14, 9) = .0, (14, 10) = .0, (14, 11) = .0, (14, 12) = .0, (14, 13) = .0, (14, 14) = .0, (14, 15) = .0, (14, 16) = .0, (14, 17) = .0, (14, 18) = .0, (14, 19) = .0, (14, 20) = .0, (15, 1) = .0, (15, 2) = .0, (15, 3) = .0, (15, 4) = .0, (15, 5) = .0, (15, 6) = .0, (15, 7) = .0, (15, 8) = .0, (15, 9) = .0, (15, 10) = .0, (15, 11) = .0, (15, 12) = .0, (15, 13) = .0, (15, 14) = .0, (15, 15) = .0, (15, 16) = .0, (15, 17) = .0, (15, 18) = .0, (15, 19) = .0, (15, 20) = .0, (16, 1) = .0, (16, 2) = .0, (16, 3) = .0, (16, 4) = .0, (16, 5) = .0, (16, 6) = .0, (16, 7) = .0, (16, 8) = .0, (16, 9) = .0, (16, 10) = .0, (16, 11) = .0, (16, 12) = .0, (16, 13) = .0, (16, 14) = .0, (16, 15) = .0, (16, 16) = .0, (16, 17) = .0, (16, 18) = .0, (16, 19) = .0, (16, 20) = .0, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = .0, (17, 5) = .0, (17, 6) = .0, (17, 7) = .0, (17, 8) = .0, (17, 9) = .0, (17, 10) = .0, (17, 11) = .0, (17, 12) = .0, (17, 13) = .0, (17, 14) = .0, (17, 15) = .0, (17, 16) = .0, (17, 17) = .0, (17, 18) = .0, (17, 19) = .0, (17, 20) = .0, (18, 1) = .0, (18, 2) = .0, (18, 3) = .0, (18, 4) = .0, (18, 5) = .0, (18, 6) = .0, (18, 7) = .0, (18, 8) = .0, (18, 9) = .0, (18, 10) = .0, (18, 11) = .0, (18, 12) = .0, (18, 13) = .0, (18, 14) = .0, (18, 15) = .0, (18, 16) = .0, (18, 17) = .0, (18, 18) = .0, (18, 19) = .0, (18, 20) = .0, (19, 1) = .0, (19, 2) = .0, (19, 3) = .0, (19, 4) = .0, (19, 5) = .0, (19, 6) = .0, (19, 7) = .0, (19, 8) = .0, (19, 9) = .0, (19, 10) = .0, (19, 11) = .0, (19, 12) = .0, (19, 13) = .0, (19, 14) = .0, (19, 15) = .0, (19, 16) = .0, (19, 17) = .0, (19, 18) = .0, (19, 19) = .0, (19, 20) = .0, (20, 1) = .0, (20, 2) = .0, (20, 3) = .0, (20, 4) = .0, (20, 5) = .0, (20, 6) = .0, (20, 7) = .0, (20, 8) = .0, (20, 9) = .0, (20, 10) = .0, (20, 11) = .0, (20, 12) = .0, (20, 13) = .0, (20, 14) = .0, (20, 15) = .0, (20, 16) = .0, (20, 17) = .0, (20, 18) = .0, (20, 19) = .0, (20, 20) = .0}, datatype = float[8], order = C_order), Array(1..20, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (13, 1) = .0, (13, 2) = .0, (13, 3) = .0, (13, 4) = .0, (13, 5) = .0, (13, 6) = .0, (14, 1) = .0, (14, 2) = .0, (14, 3) = .0, (14, 4) = .0, (14, 5) = .0, (14, 6) = .0, (15, 1) = .0, (15, 2) = .0, (15, 3) = .0, (15, 4) = .0, (15, 5) = .0, (15, 6) = .0, (16, 1) = .0, (16, 2) = .0, (16, 3) = .0, (16, 4) = .0, (16, 5) = .0, (16, 6) = .0, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = .0, (17, 5) = .0, (17, 6) = .0, (18, 1) = .0, (18, 2) = .0, (18, 3) = .0, (18, 4) = .0, (18, 5) = .0, (18, 6) = .0, (19, 1) = .0, (19, 2) = .0, (19, 3) = .0, (19, 4) = .0, (19, 5) = .0, (19, 6) = .0, (20, 1) = .0, (20, 2) = .0, (20, 3) = .0, (20, 4) = .0, (20, 5) = .0, (20, 6) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 0, (20) = 0}, datatype = integer[8]), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..40, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0, (21) = .0, (22) = .0, (23) = .0, (24) = .0, (25) = .0, (26) = .0, (27) = .0, (28) = .0, (29) = .0, (30) = .0, (31) = .0, (32) = .0, (33) = .0, (34) = .0, (35) = .0, (36) = .0, (37) = .0, (38) = .0, (39) = .0, (40) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 0, (20) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order), Array(1..20, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = -121.0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..20, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (1, 8) = .0, (1, 9) = .0, (1, 10) = .0, (1, 11) = .0, (1, 12) = .0, (1, 13) = .0, (1, 14) = .0, (1, 15) = .0, (1, 16) = .0, (1, 17) = .0, (1, 18) = .0, (1, 19) = .0, (1, 20) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = .0, (2, 9) = .0, (2, 10) = .0, (2, 11) = .0, (2, 12) = .0, (2, 13) = .0, (2, 14) = .0, (2, 15) = .0, (2, 16) = .0, (2, 17) = .0, (2, 18) = .0, (2, 19) = .0, (2, 20) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (3, 8) = .0, (3, 9) = .0, (3, 10) = .0, (3, 11) = .0, (3, 12) = .0, (3, 13) = .0, (3, 14) = .0, (3, 15) = .0, (3, 16) = .0, (3, 17) = .0, (3, 18) = .0, (3, 19) = .0, (3, 20) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (4, 8) = .0, (4, 9) = .0, (4, 10) = .0, (4, 11) = .0, (4, 12) = .0, (4, 13) = .0, (4, 14) = .0, (4, 15) = .0, (4, 16) = .0, (4, 17) = .0, (4, 18) = .0, (4, 19) = .0, (4, 20) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (5, 8) = .0, (5, 9) = .0, (5, 10) = .0, (5, 11) = .0, (5, 12) = .0, (5, 13) = .0, (5, 14) = .0, (5, 15) = .0, (5, 16) = .0, (5, 17) = .0, (5, 18) = .0, (5, 19) = .0, (5, 20) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (6, 8) = .0, (6, 9) = .0, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0, (6, 13) = .0, (6, 14) = .0, (6, 15) = .0, (6, 16) = .0, (6, 17) = .0, (6, 18) = .0, (6, 19) = .0, (6, 20) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = u[1](y), Y[2] = diff(u[1](y),y), Y[3] = u[2](y), Y[4] = diff(u[2](y),y), Y[5] = u[3](y), Y[6] = diff(u[3](y),y), Y[7] = u[4](y), Y[8] = diff(u[4](y),y), Y[9] = u[5](y), Y[10] = diff(u[5](y),y), Y[11] = u[6](y), Y[12] = diff(u[6](y),y), Y[13] = u[7](y), Y[14] = diff(u[7](y),y), Y[15] = u[8](y), Y[16] = diff(u[8](y),y), Y[17] = u[9](y), Y[18] = diff(u[9](y),y), Y[19] = u[10](y), Y[20] = diff(u[10](y),y)]`; YP[2] := -(242/3)*Y[3]+(242/3)*Y[1]+Y[1]^2; YP[4] := Y[3]^2-121*Y[1]+242*Y[3]-121*Y[5]; YP[6] := Y[5]^2-121*Y[3]+242*Y[5]-121*Y[7]; YP[8] := Y[7]^2-121*Y[5]+242*Y[7]-121*Y[9]; YP[10] := Y[9]^2-121*Y[7]+242*Y[9]-121*Y[11]; YP[12] := Y[11]^2-121*Y[9]+242*Y[11]-121*Y[13]; YP[14] := Y[13]^2-121*Y[11]+242*Y[13]-121*Y[15]; YP[16] := Y[15]^2-121*Y[13]+242*Y[15]-121*Y[17]; YP[18] := Y[17]^2-121*Y[15]+242*Y[17]-121*Y[19]; YP[20] := Y[19]^2-121*Y[17]+242*Y[19]-121; YP[1] := Y[2]; YP[3] := Y[4]; YP[5] := Y[6]; YP[7] := Y[8]; YP[9] := Y[10]; YP[11] := Y[12]; YP[13] := Y[14]; YP[15] := Y[16]; YP[17] := Y[18]; YP[19] := Y[20]; 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] = u[1](y), Y[2] = diff(u[1](y),y), Y[3] = u[2](y), Y[4] = diff(u[2](y),y), Y[5] = u[3](y), Y[6] = diff(u[3](y),y), Y[7] = u[4](y), Y[8] = diff(u[4](y),y), Y[9] = u[5](y), Y[10] = diff(u[5](y),y), Y[11] = u[6](y), Y[12] = diff(u[6](y),y), Y[13] = u[7](y), Y[14] = diff(u[7](y),y), Y[15] = u[8](y), Y[16] = diff(u[8](y),y), Y[17] = u[9](y), Y[18] = diff(u[9](y),y), Y[19] = u[10](y), Y[20] = diff(u[10](y),y)]`; YP[2] := -(242/3)*Y[3]+(242/3)*Y[1]+Y[1]^2; YP[4] := Y[3]^2-121*Y[1]+242*Y[3]-121*Y[5]; YP[6] := Y[5]^2-121*Y[3]+242*Y[5]-121*Y[7]; YP[8] := Y[7]^2-121*Y[5]+242*Y[7]-121*Y[9]; YP[10] := Y[9]^2-121*Y[7]+242*Y[9]-121*Y[11]; YP[12] := Y[11]^2-121*Y[9]+242*Y[11]-121*Y[13]; YP[14] := Y[13]^2-121*Y[11]+242*Y[13]-121*Y[15]; YP[16] := Y[15]^2-121*Y[13]+242*Y[15]-121*Y[17]; YP[18] := Y[17]^2-121*Y[15]+242*Y[17]-121*Y[19]; YP[20] := Y[19]^2-121*Y[17]+242*Y[19]-121; YP[1] := Y[2]; YP[3] := Y[4]; YP[5] := Y[6]; YP[7] := Y[8]; YP[9] := Y[10]; YP[11] := Y[12]; YP[13] := Y[14]; YP[15] := Y[16]; YP[17] := Y[18]; YP[19] := Y[20]; 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..20, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0., (7) = 0., (8) = 0., (9) = 0., (10) = 0., (11) = 0., (12) = 0., (13) = 0., (14) = 0., (15) = 0., (16) = 0., (17) = 0., (18) = 0., (19) = 0., (20) = 0.}); _vmap := array( 1 .. 20, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4), ( 5 ) = (5), ( 6 ) = (6), ( 7 ) = (7), ( 9 ) = (9), ( 8 ) = (8), ( 11 ) = (11), ( 10 ) = (10), ( 13 ) = (13), ( 12 ) = (12), ( 15 ) = (15), ( 14 ) = (14), ( 18 ) = (18), ( 19 ) = (19), ( 16 ) = (16), ( 17 ) = (17), ( 20 ) = (20) ] ); _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) = [y, u[1](y), diff(u[1](y), y), u[2](y), diff(u[2](y), y), u[3](y), diff(u[3](y), y), u[4](y), diff(u[4](y), y), u[5](y), diff(u[5](y), y), u[6](y), diff(u[6](y), y), u[7](y), diff(u[7](y), y), u[8](y), diff(u[8](y), y), u[9](y), diff(u[9](y), y), u[10](y), diff(u[10](y), y)], (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

  

 

NULL


 

Download odeProb.mw

Maybe the attached...

tomleslie 13876
February 21 2020
0 1

will meet your requirements - but you should bear in mind

  1. I fixed several syntax-type errors
  2. I reorganised the calculation, so that answers could be obtaained with generating intermediate expressions whihc were too long to even be displayed
  3. I deleted everything which I'm pretty sure you don't need, and which was just producing 'clutter' and adding to execution time

In any of the above steps, it is possible (likely?)  that I have misinterpreted your intention, so you should check the attached very carefully

> 

restart
> 

"# #` Parameters whihc don't seem to change` # A:=0.2500000001: Pr:=6.2: B:=0: rho[s]:=5200: K[nf]:=0.6842: sin(omega):=-(1)/(sqrt(2)): rho[fl]:=997.1: rhobeta[fl]:=20939.1: lambda[T]:=0: phi:=0.05: "
> 

"F(0):=0: F(1):=0: F(2):=A: delta(0):=0: delta(1):=1:Theta(0):=0: Theta(1):=B: for k from 1 to 20 do delta(k):=0: od: "
> 

d := [[M = 0, alpha = 0], [M = 1, alpha = (1/4)*Pi], [M = 2, alpha = 70*Pi*(1/180)], [M = 4, alpha = (1/2)*Pi]]; for zz to 4 do for k from 0 to 20 do F(k+3) := eval(-((1/2)*(1-phi)^2.5*(1-phi+phi*rho[s]/rho[fl])*add((m+1)*(m+2)*F(m+2)*(cos(omega)*delta(k-m+1)+sin(omega)*F(k-m)), m = 0 .. k)+(1-phi)^2.5*M*sin(alpha)^2*(1-(k+1)*F(k+1))+(1-phi)^2.5*(1-phi+phi*`&rho;&beta;`[s]/`&rho;&beta;`[fl])*lambda[T]*Theta(k))/((k+1)*(k+2)*(k+3)), d[zz]); Theta(k+2) := eval(-(1/2)*Pr*add((m+1)*Theta(m+1)*(cos(omega)*delta(k-m+1)+sin(omega)*F(k-m)), m = 0 .. k)/((k+1)*(k+2)*K[nf]), d[zz]) end do; f := unapply(evalf(add(F(i)*eta^i, i = 0 .. 18)), eta); p[zz] := diff(f(eta), `$`(eta, 1)) end do

proc (eta) options operator, arrow; .2500000001*eta^2+0.7844756746e-3*eta^5+0.4835302095e-5*eta^8+0.3134856919e-7*eta^11+0.2010398018e-9*eta^14+0.1269383510e-11*eta^17 end proc

  

.5000000002*eta+0.3922378373e-2*eta^4+0.3868241676e-4*eta^7+0.3448342611e-6*eta^10+0.2814557225e-8*eta^13+0.2157951967e-10*eta^16

  

proc (eta) options operator, arrow; .2500000001*eta^2-0.7330401580e-1*eta^3-0.9163001970e-2*eta^4-0.8157969525e-2*eta^5-0.4259579497e-2*eta^6-0.2179525247e-2*eta^7-0.1379189228e-2*eta^8-0.8990709844e-3*eta^9-0.6230154706e-3*eta^10-0.4514213406e-3*eta^11-0.3374293302e-3*eta^12-0.2590468799e-3*eta^13-0.2032414607e-3*eta^14-0.1624093735e-3*eta^15-0.1318406735e-3*eta^16-0.1084964562e-3*eta^17-0.9035990296e-4*eta^18 end proc

  

.5000000002*eta-.2199120474*eta^2-0.3665200788e-1*eta^3-0.4078984762e-1*eta^4-0.2555747698e-1*eta^5-0.1525667673e-1*eta^6-0.1103351382e-1*eta^7-0.8091638860e-2*eta^8-0.6230154706e-2*eta^9-0.4965634747e-2*eta^10-0.4049151962e-2*eta^11-0.3367609439e-2*eta^12-0.2845380450e-2*eta^13-0.2436140602e-2*eta^14-0.2109450776e-2*eta^15-0.1844439755e-2*eta^16-0.1626478253e-2*eta^17

  

proc (eta) options operator, arrow; .2500000001*eta^2-.2589162997*eta^3-0.3236453744e-1*eta^4-0.4521844937e-1*eta^5-0.1624666455e-1*eta^6-0.8552037902e-2*eta^7-0.5183881866e-2*eta^8-0.3127148195e-2*eta^9-0.2192623648e-2*eta^10-0.1577563635e-2*eta^11-0.1180719689e-2*eta^12-0.9087173806e-3*eta^13-0.7133440500e-3*eta^14-0.5704741414e-3*eta^15-0.4632747183e-3*eta^16-0.3813496501e-3*eta^17-0.3176761165e-3*eta^18 end proc

  

.5000000002*eta-.7767488991*eta^2-.1294581498*eta^3-.2260922468*eta^4-0.9747998730e-1*eta^5-0.5986426531e-1*eta^6-0.4147105493e-1*eta^7-0.2814433376e-1*eta^8-0.2192623648e-1*eta^9-0.1735319998e-1*eta^10-0.1416863627e-1*eta^11-0.1181332595e-1*eta^12-0.9986816700e-2*eta^13-0.8557112121e-2*eta^14-0.7412395493e-2*eta^15-0.6482944052e-2*eta^16-0.5718170097e-2*eta^17

  

proc (eta) options operator, arrow; .2500000001*eta^2-.5864321264*eta^3-0.7330401576e-1*eta^4-.1610295286*eta^5-0.4159950620e-1*eta^6-0.2700577342e-1*eta^7-0.1321116627e-1*eta^8-0.6667182731e-2*eta^9-0.4901942912e-2*eta^10-0.3305156270e-2*eta^11-0.2568387529e-2*eta^12-0.1997823260e-2*eta^13-0.1580207564e-2*eta^14-0.1273768591e-2*eta^15-0.1035785272e-2*eta^16-0.8541271071e-3*eta^17-0.7119515940e-3*eta^18 end proc

  

.5000000002*eta-1.759296379*eta^2-.2932160630*eta^3-.8051476430*eta^4-.2495970372*eta^5-.1890404139*eta^6-.1056893302*eta^7-0.6000464458e-1*eta^8-0.4901942912e-1*eta^9-0.3635671897e-1*eta^10-0.3082065035e-1*eta^11-0.2597170238e-1*eta^12-0.2212290590e-1*eta^13-0.1910652886e-1*eta^14-0.1657256435e-1*eta^15-0.1452016082e-1*eta^16-0.1281512869e-1*eta^17

 (1)
> 

plot([p[1], p[2], p[3], p[4]], eta = 0 .. 2, 0 .. 1.5)

 
> 

plot([p[1], p[2], p[3], p[4]], eta = 0 .. 2, -4 .. 1.5)

 
> 

 

Download DTM3.mw

First 82 83 84 85 86 87 88 Last Page 84 of 207