Advanced Search

mmcdara

6650 Reputation

18 Badges

8 years, 133 days
Contact mmcdara

MaplePrimes Activity


These are answers submitted by mmcdara

@hajnayeb Hi again, This reply...

mmcdara 6650
August 05 2020
1 7

@hajnayeb 

Hi again, 

This reply doesn't account for your last one (I had prepared mine and I was about to send it when I read yours).
So do not take it as a final answer to your problem but as a quick survey of what can be done.

PS: I did not use a white noise but a slightly colored one : the third argument in procedure KG (choosen to 0.1) is the correlation length of the random process [should be 0 for a true white noise but this choice will not do correct results with the attached file]).
 

> 

restart:
> 

with(LinearAlgebra):
with(plots):
with(Statistics):
> 

# Step 1: compute the "impulse response function" (IRF) defined as the solution of
# the ode with a Dirac source term at time t.
# To simplify I assume here the ic for the original problem are x(0)=D(x)(0)=0.
#
# It can be shown this IRF is the solution of an homogeneous ode complemented with
# non homogeneous ic:

homogeneous_edo    := m*diff(x(t),t$2)+c*diff(x(t),t)+k*x(t)=0;
non_homogeneous_ic := x(0)=0, D(x)(0)=1/m;
IRF                := unapply( rhs(dsolve([homogeneous_edo, non_homogeneous_ic],x(t))), t);

m*(diff(diff(x(t), t), t))+c*(diff(x(t), t))+k*x(t) = 0

  

x(0) = 0, (D(x))(0) = 1/m

  

proc (t) options operator, arrow; exp((1/2)*(-c+(c^2-4*k*m)^(1/2))*t/m)/(c^2-4*k*m)^(1/2)-exp(-(1/2)*(c+(c^2-4*k*m)^(1/2))*t/m)/(c^2-4*k*m)^(1/2) end proc

 (1)
> 

# Step 2: Let F(t) the random excitation and let mu__F(t) its mean and R__F(t__1, t__2) it's
# autocorrelation function.
# Here are some results

# Response X(t) of the system under the random loading F(t)

X(t) = Int(F(tau)*'IRF'(t-tau), tau=-infinity..+infinity);

X(t) = Int(F(tau)*IRF(t-tau), tau = -infinity .. infinity)

 (2)
> 

# Mean value of the response

mu__X(t) = Int(mu__F(tau)*'IRF'(t-tau), tau=-infinity..+infinity);
mu__X(t) = Int(mu__F(t-tau)*'IRF'(tau), tau=-infinity..+infinity);  #equivalently

# if F(t) is stationnary so is X(t)

mu__X = mu__F*Int('IRF'(tau), tau=-infinity..+infinity);

mu__X(t) = Int(mu__F(tau)*IRF(t-tau), tau = -infinity .. infinity)

  

mu__X(t) = Int(mu__F(t-tau)*IRF(tau), tau = -infinity .. infinity)

  

mu__X = mu__F*(Int(IRF(tau), tau = -infinity .. infinity))

 (3)
> 

# Autocorrelation function of X(t)

R__X(t__1, t__2) = Int(Int(R__X(tau__1, tau__2)*'IRF'(t__1-tau__1)*'IRF'(t__2-tau__2), tau__1=-infinity..+infinity), tau__1=-infinity..+infinity);

R__X(t__1, t__2) = Int(Int(R__X(tau__1, tau__2)*IRF(t__1-tau__1)*IRF(t__2-tau__2), tau__1 = -infinity .. infinity), tau__1 = -infinity .. infinity)

 (4)
> 

# Spectral density S__X(xi)
# Let H(xi) the Fourier transform of IRF(tau) and S__F(xi) the spectral density of F(t) then

S__X(xi) = conjugate(H(xi))*H(xi)*S__F(xi);

S__X(xi) = conjugate(H(xi))*H(xi)*S__F(xi)

 (5)
> 

# How to draw trajectories: a notional example
# F(t) is a stationnary gaussian random process (look how its point realizations F are constructed),
# with gaussian correlation function (look the expression of the variable K in procedure KG)
# (be free to ask me all precisions you would need).

# Procedure KG generates a random realization of F(t).
# For each such realization there exist a trajectory X(t), solution of the ode with random
# excitation KG(...).
# Nb_traj=10 such trajectories are drawn


KG := proc(X, Y, psi, sigma)
  local NX, DX, K, mu, k, y:
  NX := numelems(X);
  DX := < seq(Vector[row](NX, X), k=1..NX) >^+:
  K  := (sigma, psi) -> evalf( sigma^2 *~ exp~(-((DX - DX^+) /~ psi)^~2) ):
  mu := add(Y) / NX;
  k  := (t, sigma, psi) -> evalf( convert(sigma^2 *~ exp~(-((t -~ X ) /~ psi)^~2), Vector[row]) ):
  y  := mu + k(t, sigma, psi) . (K(sigma, psi))^(-1) . convert((Y -~ mu), Vector[column]):
  return y
end proc:
> 

NX := 10:
T  := [$1..NX]:

Nb_traj := 10:

for traj from 1 to Nb_traj do
  F  := convert( Statistics:-Sample(Normal(0, 1), NX), list):

  X := dsolve(
         {
           m*diff(x(t),t$2)+c*diff(x(t),t)+k*x(t)=KG(T, F, 0.1, 1), x(0)=1, D(x)(0)=0
         },
         numeric, parameters=[m, c, k]
       ):

  X(parameters=[1, 1, 1]):
  pl||traj := odeplot(X, [t, x(t)], t=0..10, color=ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]))
end do:

display(seq(pl||traj, traj=1..Nb_traj))
 

 
> 

 


 

Download RandomSourceTerm.mw

@Rouben Rostamian  @EarlA slight si...

mmcdara 6650
August 04 2020
2 2

@Rouben Rostamian  @Earl

A slight simplification which avoids using the parallel axis theorem
Ic := int(int(((x-x__c)^2+y^2), y=-f(x)..f(x)), x=0..1);

And a funny (maybe?) different point of view 


 

> 

restart;
> 

with(Statistics):
> 

f := sqrt(x^(2/3) - x^2)

(x^(2/3)-x^2)^(1/2)

 (1)
> 

# Normalization constant

C := int(x, x=0..1);

1/2

 (2)
> 

# identify f/C to the pdf of some random variable X

pdf := unapply(piecewise(x<0, 0, x < 1, f/C, 0), x);
X   := Distribution(PDF = pdf);

proc (x) options operator, arrow; piecewise(x < 0, 0, x < 1, 2*(x^(2/3)-x^2)^(1/2), 0) end proc

  

_m4548517376

 (3)
> 

# Note that the position of the center of mass doesn't change if you consider only
# the part of the lamina defined by y >=0.
# So, in Rouben's notations

x__c := Mean(X);
evalf(x__c);

(2/5)*GAMMA(3/4)^2*2^(1/2)/Pi^(1/2)

  

.4792560936

 (4)
> 

# The variance of X is equal to the inertia of the lamina around the vertical axis
# defined by x = x__c

I__x__c := Variance(X);
evalf(%)

(1/800)*(-256*GAMMA(3/4)^4+75*Pi^2)/Pi

  

0.6483790821e-1

 (5)
> 

# To obtain the inertia around the center of mass of the lamina, whose position is (x__c, 0)
# you need to add a term I__y__c wich corresponds to Rouben's int(int(y^2, y=-f(x)..f(x)), x=0..1).
# In the statistical framework int(y^2, y=-f(x)..f(x)) is the variance V__y(x) of a random variable
# uniformly distributed with support [-f(x), f(x)].
# This variance can be interpreted as a function of the random variable whose distribution is X.
# Takink its mean just gives the value of I__y__c


g       := unapply(f, x);
V__y    := eval(Variance(Uniform(-g(x), g(x))), x=RandomVariable(X)):
I__y__c := Mean(V__y);

proc (x) options operator, arrow; (x^(2/3)-x^2)^(1/2) end proc

  

(1/32)*Pi

 (6)
> 

I__c := I__x__c + I__y__c;
evalf(%);

(1/800)*(-256*GAMMA(3/4)^4+75*Pi^2)/Pi+(1/32)*Pi

  

.1630126786

 (7)
> 

 


 

Download A_statistical_point_of_view.mw

Hi, Customize the code below as you...

mmcdara 6650
July 31 2020
1 1

Hi, 

Customize the code below as you want

> 

restart:
> 

interface(rtablesize=10):
> 

MyErrorPlot := proc(x, y, Ey, symb, symbsize, symbcol, errwidth, errcol)
  # symb     : symbol corresponding tocouples (x[i], y[i])
  # symbsize : its size
  # symbcol  : its color
  # errwidth : width of the horizontal bars
  # errcol   : color of the error bars
  plots:-display(
    Statistics:-ScatterPlot(x, y, symbol=symb, symbolsize=symbsize, color=symbcol),
    seq(plot([[x[n], y[n]-Ey[n]], [x[n], y[n]+Ey[n]]], color=errcol), n=1..numelems(x)),
    seq(plot([[x[n]-errwidth, y[n]-Ey[n]], [x[n]+errwidth, y[n]-Ey[n]]], color=errcol), n=1..numelems(x)),
    seq(plot([[x[n]-errwidth, y[n]+Ey[n]], [x[n]+errwidth, y[n]+Ey[n]]], color=errcol), n=1..numelems(x))
  )
end proc:
> 

A := [$1..10]:
B := [seq(sqrt(i), i = 1..10)]:
C := [seq(i/10, i = 1..10)]:

MyErrorPlot(A, B, C, circle, 15, red, 0.2, green)

 
> 

 


 

Download MyErrorPlot.mw

 

 

 

Hi, Here is another porposition.In ...

mmcdara 6650
July 29 2020
0 0

Hi, 

Here is another proposition.
In this one the circle is drawn by drawing two disjoint arcs: onre is outside the "hill", and thus visible, the other is "inside" the hill and invisible while the "hill" is opaque.
These two arc are drawn with separate options (either spacecurve with different values of  linestyle or either tubeplot of different colors).
Of course it's less immediate that all that had been presented previously, but this can show you how you can customize your own plots.

To make the graphics to appear set to true the variable EnableGraphics at the beginning of the worksheet.

Split_Circle.mw



PS : I used the "hill" function defined by Kitonum

Hu,  Probably lengthy, crude and n...

mmcdara 6650
July 19 2020
1 0

Hu, 

Probably lengthy, crude and not top programming... but it seems to work


 

> 

restart:
> 

rewrite:=proc(expr)
  local terms, REW, rewritten, eq, eqloc, eqs, deq, der, expo, der0, ind, D_deg, deg:

  if type(expr, `+`) then
    terms := [op(expr)]:
  else
    terms := expr:
  end if:


  REW := 0:
  for eq in terms do
    if type(eq, `*`) then
      eqs := op(eq)
    else
      eqs := copy(eq)
    end if:

    rewritten := 1:
    for eqloc in [eqs] do  
      deq       := convert(eqloc, 'D'):
      if has(deq, 'D') then
        if type(deq, `^`) then
          der, expo := op(deq);
        else
          der, expo := copy(deq), 1
        end if;
        der0      := op(0, der);
        ind       := op(op(der0));
        D_deg     := op(op(0, der0));
        if numelems([D_deg]) = 1 then deg:=1 else deg := D_deg[2] end if;

        rewritten := rewritten * (f[ind-deg](x))^expo:
      else
        rewritten := rewritten * deq:
      end if:
    end do;

    REW := REW + rewritten
  end do:
  return REW:
end proc:

  
> 

eq := f[4](x)*diff(f[3](x),x$2)^5 + diff(f[-2](x), x)*diff(f[2](x), x$3);

rewrite(eq);

f[4](x)*(diff(diff(f[3](x), x), x))^5+(diff(f[-2](x), x))*(diff(diff(diff(f[2](x), x), x), x))

  

f[4](x)*f[1](x)^5+f[-3](x)*f[-1](x)

 (1)
> 

 


 

Download proposal.mw

@Violet Just a little error in your...

mmcdara 6650
July 19 2020
1 2

@Violet 

Just a little error in your R code:

library("deSolve", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")

spirs <- function(t, x, parms) {
  with(as.list(c(parms, x)), {
    dS=L-(1-u)*b*S*P+th*R-nu*S
    dI=(1-u)*b*S*P-(nu+sg+v+m+a-ph)*I
    dP= m*I-g1*P
    dR=(v+a)*I-(nu+th)*R
# the order of the unknonws to return must be the same that the one in x
    list( c(dS, dP, dI, dR) )  
  })
}

x       <- c(S=1000,P=100,I=150,R=0)
parms <- c(L=0.05,b=0.002,ph=0.01,a=0.01,m=0.01,nu=0.001,th=0.005,u=0.1, v=0.1,sg=0.01, g1=0.03)
times <- 1:60
run_d <- ode(x, time, spirs, parms, method = rkMethod("rk45f"))

 

 

PS : I've done years ago a lot of comparisons between  R deSolve and Maple dsolve[numeric] (18 and 2015). and the results, as well as computational times were very close.
nevertheless, if you want to use the lsode-lsode suite of methods, I advice you to use R instead of Maple.

Here is what you want. (PS: the help pa...

mmcdara 6650
June 30 2020
1 1

Here is what you want.

(PS: the help pages about events are not the clearer pages in Maple, in my own opinion)

> 

restart
> 

ode2 := diff(varphi(t), t, t)+omega^2*sin(varphi(t));

p0 := evalf((10*(1/180))*Pi);
te := 6;
event1 := [[diff(varphi(t), t), halt]];

ld := dsolve([eval(ode2, omega = 2*Pi), varphi(0) = p0, (D(varphi))(0) = 0], numeric, range = 0 .. te, events = event1);

diff(diff(varphi(t), t), t)+omega^2*sin(varphi(t))

  

.1745329252

  

6

  

[[diff(varphi(t), t), halt]]

  

Warning, cannot evaluate the solution further right of .50095360, event #1 triggered a halt

  

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, ( "right" ) = 6., ( "left" ) = 0. ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([Array(1..2, 1..21, {(1, 1) = 1.0, (1, 2) = .0, (1, 3) = 1.0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = 1.0, (1, 8) = undefined, (1, 9) = .0, (1, 10) = 1.0, (1, 11) = undefined, (1, 12) = undefined, (1, 13) = undefined, (1, 14) = undefined, (1, 15) = undefined, (1, 16) = undefined, (1, 17) = undefined, (1, 18) = undefined, (1, 19) = undefined, (1, 20) = undefined, (1, 21) = undefined, (2, 1) = 1.0, (2, 2) = .0, (2, 3) = 100.0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = undefined, (2, 9) = undefined, (2, 10) = 0.10e-6, (2, 11) = undefined, (2, 12) = .0, (2, 13) = undefined, (2, 14) = .0, (2, 15) = .0, (2, 16) = undefined, (2, 17) = undefined, (2, 18) = undefined, (2, 19) = undefined, (2, 20) = undefined, (2, 21) = undefined}, datatype = float[8], order = C_order), proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[2]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, Array(1..1, 1..2, {(1, 1) = undefined, (1, 2) = undefined}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 1, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 1, (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}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .5009536044156753, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.7363088496439619e-3, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..2, {(1) = .1745329252, (2) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174532275382763, (2) = 0.3164380549290857e-2}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0.20166160408230382e-16, (2) = 6.855358404355688}, 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..6, {(1, 1) = 0.20166160408230382e-16, (1, 2) = -0.3645741252411073e-2, (1, 3) = -0.12134604818129788e-2, (1, 4) = 0.9449858716279857e-2, (1, 5) = 0.10945612556221545e-1, (1, 6) = 0.12190325843795934e-2, (2, 1) = 6.855358404355688, (2, 2) = 6.855387810726147, (2, 3) = 6.855354230575938, (2, 4) = 6.855105093491461, (2, 5) = 6.8550184393577585, (2, 6) = 6.8553542038169235}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0.3112864904113937e-16, (2) = 6.855358404355692}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = .1745329252, (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) = -6.855355274192771}, 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] = varphi(t), Y[2] = diff(varphi(t),t)]`; YP[2] := -39.4784176043574*sin(Y[1]); YP[1] := Y[2]; 0 end proc, -1, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[2]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = varphi(t), Y[2] = diff(varphi(t),t)]`; YP[2] := -39.4784176043574*sin(Y[1]); YP[1] := Y[2]; 0 end proc, -1, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[2]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )), ( 3 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([Array(1..2, 1..21, {(1, 1) = 1.0, (1, 2) = .0, (1, 3) = 1.0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = 1.0, (1, 8) = .5009536044156753, (1, 9) = 0.20166160408230382e-16, (1, 10) = 1.0, (1, 11) = -0.36456738311832493e-2, (1, 12) = 1.0, (1, 13) = 0.12186516886744188e-2, (1, 14) = 1.0, (1, 15) = 0.6082953353156279e-2, (1, 16) = 1.0, (1, 17) = 0.10947135945331345e-1, (1, 18) = 1.0, (1, 19) = 0.20166160408230382e-16, (1, 20) = 1.0, (1, 21) = .5009536044156753, (2, 1) = 1.0, (2, 2) = .0, (2, 3) = 100.0, (2, 4) = 1.0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = .5009536044156753, (2, 9) = undefined, (2, 10) = 0.10e-6, (2, 11) = undefined, (2, 12) = 1.0, (2, 13) = undefined, (2, 14) = .0, (2, 15) = .0, (2, 16) = undefined, (2, 17) = undefined, (2, 18) = undefined, (2, 19) = undefined, (2, 20) = undefined, (2, 21) = .5009536044156753}, datatype = float[8], order = C_order), proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[2]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, Array(1..1, 1..2, {(1, 1) = undefined, (1, 2) = undefined}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 101, (10) = 1, (11) = 93, (12) = 93, (13) = 0, (14) = 0, (15) = 0, (16) = 1, (17) = 1, (18) = 178, (19) = 30000, (20) = 1, (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) = 1, (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}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .5009536044156753, (2) = 0.10e-5, (3) = 0.22706126977408814e-1, (4) = 0.500001e-14, (5) = .0, (6) = 0.7363088496439619e-3, (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) = .5025505036867043, (20) = 0.22706126977408814e-1, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..2, {(1) = .1745329252, (2) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174532275382763, (2) = 0.3164380549290857e-2}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0.20166160408230382e-16, (2) = 6.855358404355688}, 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..6, {(1, 1) = 0.20166160408230382e-16, (1, 2) = -0.3645741252411073e-2, (1, 3) = -0.12134604818129788e-2, (1, 4) = 0.9449858716279857e-2, (1, 5) = 0.10945612556221545e-1, (1, 6) = 0.12190325843795934e-2, (2, 1) = 6.855358404355688, (2, 2) = 6.855387810726147, (2, 3) = 6.855354230575938, (2, 4) = 6.855105093491461, (2, 5) = 6.8550184393577585, (2, 6) = 6.8553542038169235}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0.3112864904113937e-16, (2) = 6.855358404355692}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = -.174533005711098, (2) = 0.20166160408230382e-16}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0.20166160408230382e-16, (2) = 6.855358404355688}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .5009536044156753, (1, 2) = -.174533005711098, (2, 0) = -.174533005711098, (2, 1) = 0.20166160408230382e-16, (2, 2) = .500022579464815, (3, 0) = .500022579464815, (3, 1) = -.17453003458153782, (3, 2) = -0.6382473873356864e-2, (4, 0) = -0.6382473873356864e-2, (4, 1) = .5003329211151017, (4, 2) = -.17453168520700998, (5, 0) = -.17453168520700998, (5, 1) = -0.4254995859433421e-2, (5, 2) = .5006432627653885, (6, 0) = .5006432627653885, (6, 1) = -.17453267558476632, (6, 2) = -0.2127501912837157e-2}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = varphi(t), Y[2] = diff(varphi(t),t)]`; YP[2] := -39.4784176043574*sin(Y[1]); YP[1] := Y[2]; 0 end proc, -1, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[2]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, 0, 0]), ( 13 ) = (), ( 12 ) = (Array(1..93, 0..2, {(1, 1) = .0, (1, 2) = .1745329252, (2, 0) = .1745329252, (2, 1) = .0, (2, 2) = 0.18407721241099048e-3, (3, 0) = 0.18407721241099048e-3, (3, 1) = .17453280905514362, (3, 2) = -0.12619144118906497e-2, (4, 0) = -0.12619144118906497e-2, (4, 1) = 0.36815442482198096e-3, (4, 2) = .17453246062072755, (5, 0) = .17453246062072755, (5, 1) = -0.2523827161359515e-2, (5, 2) = 0.5522316372329715e-3, (6, 0) = 0.5522316372329715e-3, (6, 1) = .17453187989721072, (6, 2) = -0.37857365859869344e-2, (7, 0) = -0.37857365859869344e-2, (7, 1) = 0.7363088496439619e-3, (7, 2) = .1745310668853582, (8, 0) = .1745310668853582, (8, 1) = -0.5047641023358316e-2, (8, 2) = 0.5316744108339218e-2, (9, 0) = 0.5316744108339218e-2, (9, 1) = .17443604137810215, (9, 2) = -0.36441499288392076e-1, (10, 0) = -0.36441499288392076e-1, (10, 1) = 0.9897179367034475e-2, (10, 2) = .17419727674301455, (11, 0) = .17419727674301455, (11, 1) = -0.6780563023559488e-1, (11, 2) = 0.1447761462572973e-1, (12, 0) = 0.1447761462572973e-1, (12, 1) = .17381496756579243, (12, 2) = -0.9911444690216951e-1, (13, 0) = -0.9911444690216951e-1, (13, 1) = 0.19058049884424985e-1, (13, 2) = .17328942567803401, (14, 0) = .17328942567803401, (14, 1) = -.13034241375896907, (14, 2) = 0.24110111233386848e-1, (15, 0) = 0.24110111233386848e-1, (15, 1) = .17254417523942078, (15, 2) = -.16466145604758078, (16, 0) = -.16466145604758078, (16, 1) = 0.2916217258234871e-1, (16, 2) = .17162594173155682, (17, 0) = .17162594173155682, (17, 1) = -.19881705207234882, (17, 2) = 0.3421423393131057e-1, (18, 0) = 0.3421423393131057e-1, (18, 1) = .17053563643790248, (18, 2) = -.23277526397700765, (19, 0) = -.23277526397700765, (19, 1) = 0.3926629528027244e-1, (19, 2) = .1692743415533087, (20, 0) = .1692743415533087, (20, 1) = -.2665023497538466, (20, 2) = 0.44940093082728484e-1, (21, 0) = 0.44940093082728484e-1, (21, 1) = .16765553624781762, (21, 2) = -.3040629406279465, (22, 0) = -.3040629406279465, (22, 1) = 0.5061389088518453e-1, (22, 2) = .16582467913627932, (23, 0) = .16582467913627932, (23, 1) = -.3412425456526048, (23, 2) = 0.5628768868764058e-1, (24, 0) = 0.5628768868764058e-1, (24, 1) = .16378406427193246, (24, 2) = -.3779944413829568, (25, 0) = -.3779944413829568, (25, 1) = 0.61961486490096626e-1, (25, 2) = .1615362486213279, (26, 0) = .1615362486213279, (26, 1) = -.41427242828499644, (26, 2) = 0.6811600578387549e-1, (27, 0) = 0.6811600578387549e-1, (27, 1) = .1588669885255837, (27, 2) = -.453035224354377, (28, 0) = -.453035224354377, (28, 1) = 0.7427052507765435e-1, (28, 2) = .15596119261200517, (29, 0) = .15596119261200517, (29, 1) = -.4911291584416978, (29, 2) = 0.804250443714332e-1, (30, 0) = 0.804250443714332e-1, (30, 1) = .15282315170324726, (30, 2) = -.5284976437938735, (31, 0) = -.5284976437938735, (31, 1) = 0.8657956366521206e-1, (31, 2) = .1494574995825724, (32, 0) = .1494574995825724, (32, 1) = -.5650851392367926, (32, 2) = 0.9315390814609589e-1, (33, 0) = 0.9315390814609589e-1, (33, 1) = .14561645758183764, (33, 2) = -.6032441918562621, (34, 0) = -.6032441918562621, (34, 1) = 0.9972825262697972e-1, (34, 2) = .14152785980868893, (35, 0) = .14152785980868893, (35, 1) = -.6403849163333629, (35, 2) = .10630259710786355, (36, 0) = .10630259710786355, (36, 1) = .13719861013811915, (36, 2) = -.6764439617474145, (37, 0) = -.6764439617474145, (37, 1) = .11287694158874738, (37, 2) = .13263601861915097, (38, 0) = .13263601861915097, (38, 1) = -.711359775906317, (38, 2) = .119845172899741, (39, 0) = .119845172899741, (39, 1) = .12755393105479498, (39, 2) = -.747053013272364, (40, 0) = -.747053013272364, (40, 1) = .12681340421073464, (40, 2) = .12222803812035075, (41, 0) = .12222803812035075, (41, 1) = -.7813260342615114, (41, 2) = .13378163552172825, (42, 0) = .13378163552172825, (42, 1) = .11666846858581881, (42, 2) = -.814112668627815, (43, 0) = -.814112668627815, (43, 1) = .1407498668327219, (43, 2) = .1108857945860264, (44, 0) = .1108857945860264, (44, 1) = -.8453495705545896, (44, 2) = .14811880982683875, (45, 0) = .14811880982683875, (45, 1) = .10454007728404852, (45, 2) = -.8766297942144323, (46, 0) = -.8766297942144323, (46, 1) = .15548775282095562, (46, 2) = 0.979707090288822e-1, (47, 0) = 0.979707090288822e-1, (47, 1) = -.906041321689006, (47, 2) = .16285669581507248, (48, 0) = .16285669581507248, (48, 1) = 0.9119169943208455e-1, (48, 2) = -.9335201399465292, (49, 0) = -.9335201399465292, (49, 1) = .17022563880918934, (49, 2) = 0.8421750146782593e-1, (50, 0) = 0.8421750146782593e-1, (50, 1) = -.9590064212629931, (50, 2) = .17804355821708154, (51, 0) = .17804355821708154, (51, 1) = 0.7662159618078511e-1, (51, 2) = -.983805265666594, (52, 0) = -.983805265666594, (52, 1) = .18586147762497374, (52, 2) = 0.6884103757579298e-1, (53, 0) = 0.6884103757579298e-1, (53, 1) = -1.0062376907723325, (53, 2) = .1936793970328659, (54, 0) = .1936793970328659, (54, 1) = 0.60894548536701315e-1, (54, 2) = -1.0262483290074216, (55, 0) = -1.0262483290074216, (55, 1) = .2014973164407581, (55, 2) = 0.5280124093257794e-1, (56, 0) = 0.5280124093257794e-1, (56, 1) = -1.0437878150898912, (56, 2) = .20894558539088942, (57, 0) = .20894558539088942, (57, 1) = 0.4497186952165007e-1, (57, 2) = -1.058159599240839, (58, 0) = -1.058159599240839, (58, 1) = .2163938543410207, (58, 2) = 0.3704406311146387e-1, (59, 0) = 0.3704406311146387e-1, (59, 1) = -1.0702166404374696, (59, 2) = .22384212329115202, (60, 0) = .22384212329115202, (60, 1) = 0.2903516958079011e-1, (60, 2) = -1.0799317760685898, (61, 0) = -1.0799317760685898, (61, 1) = .23129039224128334, (61, 2) = 0.20962701162870183e-1, (62, 0) = 0.20962701162870183e-1, (62, 1) = -1.0872831547929174, (62, 2) = .23830523221322775, (63, 0) = .23830523221322775, (63, 1) = 0.13317688199163954e-1, (63, 2) = -1.092030411266159, (64, 0) = -1.092030411266159, (64, 1) = .24532007218517216, (64, 2) = 0.5646814430533518e-2, (65, 0) = 0.5646814430533518e-2, (65, 1) = -1.0946567753827166, (65, 2) = .25233491215711656, (66, 0) = .25233491215711656, (66, 1) = -0.2035018872931373e-2, (66, 2) = -1.0951569754203698, (67, 0) = -1.0951569754203698, (67, 1) = .259349752129061, (67, 2) = -0.9712902985745047e-2, (68, 0) = -0.9712902985745047e-2, (68, 1) = -1.0935300176706553, (68, 2) = .2660678480939827, (69, 0) = .2660678480939827, (69, 1) = -0.1704851596954728e-1, (69, 2) = -1.0899807796199925, (70, 0) = -1.0899807796199925, (70, 1) = .27278594405890444, (70, 2) = -0.24353753975034524e-1, (71, 0) = -0.24353753975034524e-1, (71, 1) = -1.084490031224452, (71, 2) = .2795040400238261, (72, 0) = .2795040400238261, (72, 1) = -0.31615603084683634e-1, (72, 2) = -1.077067834464908, (73, 0) = -1.077067834464908, (73, 1) = .28622213598874785, (73, 2) = -0.3882114098584355e-1, (74, 0) = -0.3882114098584355e-1, (74, 1) = -1.0677277904389917, (74, 2) = .2930313501725285, (75, 0) = .2930313501725285, (75, 1) = -0.4605380363613945e-1, (75, 2) = -1.0563215827662154, (76, 0) = -1.0563215827662154, (76, 1) = .2998405643563091, (76, 2) = -0.53202206356303613e-1, (77, 0) = -0.53202206356303613e-1, (77, 1) = -1.042984215240012, (77, 2) = .3066497785400898, (78, 0) = .3066497785400898, (78, 1) = -0.6025327970271975e-1, (78, 2) = -1.0277407394066185, (79, 0) = -1.0277407394066185, (79, 1) = .31345899272387046, (79, 2) = -0.6719414990810116e-1, (80, 0) = -0.6719414990810116e-1, (80, 1) = -1.0106197444303908, (80, 2) = .32059167071093014, (81, 0) = .32059167071093014, (81, 1) = -0.743327598585143e-1, (81, 2) = -.9907070155189648, (82, 0) = -.9907070155189648, (82, 1) = .32772434869798983, (82, 2) = -0.8132223263390685e-1, (83, 0) = -0.8132223263390685e-1, (83, 1) = -.9688103368834072, (83, 2) = .3348570266850496, (84, 0) = .3348570266850496, (84, 1) = -0.8814857087685264e-1, (84, 2) = -.944974607823269, (85, 0) = -.944974607823269, (85, 1) = .34198970467210926, (85, 2) = -0.9479812742083839e-1, (86, 0) = -0.9479812742083839e-1, (86, 1) = -.9192485905889495, (86, 2) = .3494494187932333, (87, 0) = .3494494187932333, (87, 1) = -.10154901554387676, (87, 2) = -.8903778843165101, (88, 0) = -.8903778843165101, (88, 1) = .35690913291435744, (88, 2) = -.10807723297578052, (89, 0) = -.10807723297578052, (89, 1) = -.8595616002422856, (89, 2) = .3643688470354815, (90, 0) = .3643688470354815, (90, 1) = -.11436851385107127, (90, 2) = -.8268683615385499, (91, 0) = -.8268683615385499, (91, 1) = .37182856115660556, (91, 2) = -.12040913825656249, (92, 0) = -.12040913825656249, (92, 1) = -.7923707609023437, (92, 2) = .3796454147407636, (93, 0) = .3796454147407636, (93, 1) = -.12645567515768255, (93, 2) = -.7543689666976098}, datatype = float[8], order = C_order)), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = varphi(t), Y[2] = diff(varphi(t),t)]`; YP[2] := -39.4784176043574*sin(Y[1]); YP[1] := Y[2]; 0 end proc, -1, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[2]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )), ( 4 ) = (3) ] ); _y0 := Array(0..2, {(1) = 0., (2) = .1745329252}); _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); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _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]; _dat[4][26] := _EnvDSNumericSaveDigits; _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, varphi(t), diff(varphi(t), t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _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)
> 

# This first step is necessary.
# Use any value of T larger that the supposed time that fired the event

T :=10:
ld(T):

# now catch the time when the event fired

FiredTime := op( ld( eventfired=[1] ) );

Warning, cannot evaluate the solution further right of .50095360, event #1 triggered a halt

  

HFloat(0.5009536044156753)

 (2)
> 

plots:-odeplot(ld, [t, D(varphi)(t)], t=0..FiredTime)

 
> 

 


 

Download eventfired.mw

Is it this that you want? (the example ...

mmcdara 6650
June 28 2020
0 1

Is it this that you want?

(the example comes from the help pages)
 

> 

restart:
> 

with(CurveFitting):
> 

S := Spline([[0,0],[1,1],[2,4],[3,3]], v);

S := piecewise(v < 1, (4/5)*v^3+(1/5)*v, v < 2, -2*v^3+(42/5)*v^2-(41/5)*v+14/5, (6/5)*v^3-(54/5)*v^2+(151/5)*v-114/5)

 (1)
> 

dS := NULL:
for n from 1 to nops(S)-1 do
  if n::odd then
    dS := dS, op(n, S)
  else
    dS := dS, diff(op(n, S), v)
  end if:
end do:

dS := piecewise(dS, diff(op(-1, S), v));

dS := piecewise(v < 1, (12/5)*v^2+1/5, v < 2, -6*v^2+(84/5)*v-41/5, (18/5)*v^2-(108/5)*v+151/5)

 (2)
> 

 


 

Download Spline_derivative.mw

Probably closer to what is donne in Matl...

mmcdara 6650
June 14 2020
1 5

Probably closer to what is donne in Matlab

 

> 

restart:
> 

test:= proc(x):
> 

             return x+1, x^2:
> 

        end proc:
> 

 

> 

A, B := test(y):
A, B

y+1, y^2

 (1)
> 

# unassign A and B

A:='A': B:='B': A, B;
``, B := test(y):
A, B

A, B

  

A, y^2

 (2)
> 

# unassign A and B

A:='A': B:='B': A, B;
A, `` := test(y):
A, B

A, B

  

y+1, B

 (3)
> 

 


 

Download CloserToMatlab.mw

 

Is this what you want?  ...

mmcdara 6650
June 11 2020
0 2

Is this what you want?
(assumming derivatives with respect to x and y do commute)
 

> 

restart:
> 

alias(V=V(x, y)):
alias(alpha=alpha(x,y)):
alias(beta=beta(x,y)):
> 

cond_1 := diff(V, x)=V*alpha:
cond_2 := diff(V, y)=V*beta:
> 

eval(diff(V, x$2), cond_1);
eval(diff(V, x$2), cond_2);

(diff(V, x))*alpha+V*(diff(alpha, x))

  

diff(diff(V, x), x)

 (1)
> 

U := copy(V):
for k from 1 to 2 do
  U := eval(diff(U, y), cond_2);
end do;

V*beta

  

V*beta^2+V*(diff(beta, y))

 (2)
> 

der := proc(r, s)
  local fx, fy:

  fx := proc(W, r)
    local U, k:
    U := copy(W);
    for k from 1 to r do
      U := eval(diff(U, y), cond_2);
    end do;
    return U:
  end proc:

  fy := proc(W, s)
    local U, k:
    U := copy(W);
    for k from 1 to s do
      U := eval(diff(U, y), cond_2);
    end do;
    return U:
  end proc:

  if r=0 then
    if s=0 then
      V
    else
      Diff(V, y$s) = fy(V, s)
    end if:
  else
    if s=0 then
      Diff(V, y$r) = fx(V, r)
    else
      Diff(Diff(V, x$r), y$s) = fx(fy(V, s), r); # or  fy(fx(V, r), s) assuming derivatives commute
    end if:
  end if:
end proc

proc (r, s) local fx, fy; fx := proc (W, r) local U, k; U := copy(W); for k to r do U := eval(diff(U, y), cond_2) end do; return U end proc; fy := proc (W, s) local U, k; U := copy(W); for k to s do U := eval(diff(U, y), cond_2) end do; return U end proc; if r = 0 then if s = 0 then V else Diff(V, `$`(y, s)) = fy(V, s) end if else if s = 0 then Diff(V, `$`(y, r)) = fx(V, r) else Diff(Diff(V, `$`(x, r)), `$`(y, s)) = fx(fy(V, s), r) end if end if end proc

 (3)
> 

printlevel := 2:
for r from 0 to 2 do
  for s from 0 to 2 do
    der(r, s)
  end do:
end do;

V

  

Diff(V, y) = V*beta

  

Diff(V, y, y) = V*beta^2+V*(diff(beta, y))

  

Diff(V, y) = V*beta

  

Diff(Diff(V, x), y) = V*beta^2+V*(diff(beta, y))

  

Diff(Diff(V, x), y, y) = V*beta^3+3*V*beta*(diff(beta, y))+V*(diff(diff(beta, y), y))

  

Diff(V, y, y) = V*beta^2+V*(diff(beta, y))

  

Diff(Diff(V, x, x), y) = V*beta^3+3*V*beta*(diff(beta, y))+V*(diff(diff(beta, y), y))

  

Diff(Diff(V, x, x), y, y) = V*beta^4+6*V*beta^2*(diff(beta, y))+3*V*(diff(beta, y))^2+4*V*beta*(diff(diff(beta, y), y))+V*(diff(diff(diff(beta, y), y), y))

 (4)
> 

 


 

Download derivatives.mw

 

Your function is extremely simple for it...

mmcdara 6650
June 10 2020
1 5

Your function is extremely simple for it can be rewritten this way

ff = C + f (xx[1]) + ... f (xx[N]) 

where C is a constant and f(x) = x^2 - 10 * cos(2*Pi*x)  (A and B are two constants)
Note that the hessien of ff iis thus a diagonal matrix

Each component of the gradient of ff is defined by G := u -> 2*u + 20*Pi * sin(2*Pi*u) .
Thus each component of the gradient of ff  is null for exactly the same values of xx[n].
In the range xx[n]=-5..5 (for each n) the function u -> 2*u + 20*Pi * sin(2*Pi*u)  has 21 zeroes.
To find their values use this pice of code
Z := NULL:
z := -5:
k := 1;
while k <= 21 do
  z := RootFinding:-NextZero(u->2*u+20*Pi*sin(2*Pi*u), z);
  Z := Z, z:
  k := k+1:
end do:
Z;

The 2nd derivative of u -> 2*u - 10 * cos(2*Pi*u) is  the function H := u -> 2 + 20*Pi * sin(2*Pi*u)
Note that the hessien of ff is a diagonal matrix whose diagonal elements are just H(xx[n])
The values of H(u) for each valus of Z are calculated by
Hval := map(u -> 2+20*Pi*sin(2*Pi*u), [Z])

Among these 21 values, 13 of them are poitive :
numelems(select[flatten](`>`, Hval, 0))

Then ff has a minimum if all the components (1000) of its hessian are negative at the points where the gradient of ff is null, and has a maximum if all the components of its hessian are positve at the points where the gradient of ff is null,
As 13 values of Hval are positive, one might expect 13^1000 minima (and 8^1000 maxima) for ff.
The number of saddle points is 21^1000 - 13^1000 - 8^1000

The largest value of ff at a point where its gradient is 0, is C + 1000 * V where V is the highest value of 2*u - 10 * cos(2*Pi*u) (u in [Z]).
This value V is about 18.94... and is obtained for z = 4.52299.. (twentieth zero).
It appears that H(4.52299..) is negative.
Thus ff reaches a maximum if xx[1] = ... = xx[1000] = 4.52299..  and then takes the value 
C + 18.94.. * 1000 (which is the highest value of ff within [-5, 5]^1000

If I'm not mistaken

Here is an example (data XY come from th...

mmcdara 6650
June 03 2020
0 0

Here is an example (data XY come from the CurveFitting:-Spline help page)
The trap would be to search the extrema of each pice of the global function without accounting for the interval this piece is defined on. Another trap would be to search for extrema outsise the interval defined by the data.

 

> 

restart:
> 

with(CurveFitting):
> 

XY := [[0,0],[1,1],[2,4],[3,3]];
S := Spline([[0,0],[1,1],[2,4],[3,3]], v);

XY := [[0, 0], [1, 1], [2, 4], [3, 3]]

  

piecewise(v < 1, (4/5)*v^3+(1/5)*v, v < 2, -2*v^3+(42/5)*v^2-(41/5)*v+14/5, (6/5)*v^3-(54/5)*v^2+(151/5)*v-114/5)

 (1)
> 

# ranges

bounds := min(op~(1, XY)), seq(rhs(op(k, S)), k in [seq(1..nops(S)-1, 2)]), max(op~(1, XY));
ranges := seq(bounds[k]..bounds[k+1], k=1..numelems([bounds])-1)

0, 1, 2, 3

  

0 .. 1, 1 .. 2, 2 .. 3

 (2)
> 

# 1st derivative of the S

dS := seq(diff(op(k, S), v), k in [seq(2..nops(S), 2)]), diff(op(-1, S), v);

(12/5)*v^2+1/5, -6*v^2+(84/5)*v-41/5, (18/5)*v^2-(108/5)*v+151/5

 (3)
> 

# extrema

sol := seq([fsolve(dS[k], v=ranges[k])], k=1..numelems([dS]))

[], [], [2.218264040]

 (4)
> 

# reconstruction as a piecewise solution

piecewise(seq(op([v >= op(1, ranges[k]) and v < op(2, ranges[k]), sol[k]]), k=1..numelems([sol])))

piecewise(0 <= v and v < 1, [], 1 <= v and v < 2, [], 2 <= v and v < 3, [2.218264040])

 (5)
> 

 


 

Download extrema.mw

@JAMET  @KitonumHi JAMET, I do...

mmcdara 6650
June 02 2020
0 0

@JAMET  @Kitonum

Hi JAMET, 

I don't know if the solution Kitonum has provided suits you (I understood you were looking for something more formal?).
If I'm mistaken forget what comes below, otherwise this worksheet might interest you.
From a line L defined by its cartesian equation it contains the step by step procedure to:

  1. construct the equation of L in the complex plane
  2. construct the equation of the line orthogonal to L and which passes through (0, 0)
  3. construct the expression of the intersection of these two lines

All of this is done formally, two examples follow with drawings

 

> 

restart:
> 

alias(conj = conjugate):
> 

assumptions := (a::real, b::real, c::real, x::real, y::real)

a::real, b::real, c::real, x::real, y::real

 (1)
> 

# Cartesian equation of line L

L := a*x + b*y - c

a*x+b*y-c

 (2)
> 

# let Z a complex number

rel_1 := Z = x + I*y

Z = x+I*y

 (3)
> 

# let Z bar its conjugate

rel_2 := `#mrow(mover(mo(Z),mo("&#x305;")))` = eval(conj(Z), rel_1) assuming assumptions

`#mrow(mover(mo(Z),mo("&#x305;")))` = x-I*y

 (4)
> 

# explain x and y in terms of Z and Z bar

interface(warnlevel=0):
rel_3 := op( solve([rel_1, rel_2], [x, y]) )

[x = (1/2)*Z+(1/2)*`#mrow(mover(mo(Z),mo("&#x305;")))`, y = -((1/2)*I)*(Z-`#mrow(mover(mo(Z),mo("&#x305;")))`)]

 (5)
> 

# rewrite L by using the above equalities

eval(L, rel_3);

a*((1/2)*Z+(1/2)*`#mrow(mover(mo(Z),mo("&#x305;")))`)-((1/2)*I)*b*(Z-`#mrow(mover(mo(Z),mo("&#x305;")))`)-c

 (6)
> 

# take the numerator of this expression

numer(%);

-I*b*Z+I*b*`#mrow(mover(mo(Z),mo("&#x305;")))`+a*Z+a*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c

 (7)
> 

# collect the result according to Z and Z bar

`#mrow(mo("&#8466;"))`  := collect(%, [Z, `#mrow(mover(mo(Z),mo("&#x305;")))`])

(-I*b+a)*Z+(I*b+a)*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c

 (8)
> 

# define the complex number v from a and b

rel_4 := v = a+I*b

v = I*b+a

 (9)
> 

# let v bar its conjugate

rel_5 := `#mrow(mover(mo(v),mo("&#x305;")))` = eval(conj(v), rel_4)  assuming assumptions;

`#mrow(mover(mo(v),mo("&#x305;")))` = -I*b+a

 (10)
> 

# reverse rel_4 and rel_5

rel_6 := (rhs=~lhs)~([rel_4, rel_5])

[I*b+a = v, -I*b+a = `#mrow(mover(mo(v),mo("&#x305;")))`]

 (11)
> 

# here is the equation of L in the complex plane

`#mrow(mo("&#8466;"))` := eval(`#mrow(mo("&#8466;"))`, rel_6)

Z*`#mrow(mover(mo(v),mo("&#x305;")))`+v*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c

 (12)
> 

# find the complex equation of the line orthogonal to L and passing through the origin [0, 0]
#
# firstly define w as a rotation by Pi/2 of the vector whose affix is [a, b] in the complex plane

rel_7 := w = v*I

w = I*v

 (13)
> 

# let w bar its conjugate

rel_8 := `#mrow(mover(mo(w),mo("&#x305;")))` = eval(conj(w), rel_7);
rel_8 := eval(rel_8, conj(v) = `#mrow(mover(mo(v),mo("&#x305;")))`);

`#mrow(mover(mo(w),mo("&#x305;")))` = -I*conj(v)

  

`#mrow(mover(mo(w),mo("&#x305;")))` = -I*`#mrow(mover(mo(v),mo("&#x305;")))`

 (14)
> 

# here is the equation of the line orthogonal to L, and which passes through the origin

`#mrow(mo("&#8459;"))` := eval(
                                `#mrow(mo("&#8466;"))`,
                                [
                                  v=w,
                                  `#mrow(mover(mo(v),mo("&#x305;")))` = `#mrow(mover(mo(w),mo("&#x305;")))`,
                                  c=0
                                ]
                              )

Z*`#mrow(mover(mo(w),mo("&#x305;")))`+w*`#mrow(mover(mo(Z),mo("&#x305;")))`

 (15)
> 

# where do lines L and H intersect ?
#
# let's form this system in Z and Z bar (brute force)

sys := [ `#mrow(mo("&#8466;"))` = 0, eval(`#mrow(mo("&#8459;"))`, [rel_7, rel_8]) = 0 ]:

print~(sys):

Z*`#mrow(mover(mo(v),mo("&#x305;")))`+v*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c = 0

  

-I*Z*`#mrow(mover(mo(v),mo("&#x305;")))`+I*v*`#mrow(mover(mo(Z),mo("&#x305;")))` = 0

 (16)
> 

# where do lines L and H intersect ?

sol := op( solve( sys, [ Z, `#mrow(mover(mo(Z),mo("&#x305;")))` ] ) )

[Z = c/`#mrow(mover(mo(v),mo("&#x305;")))`, `#mrow(mover(mo(Z),mo("&#x305;")))` = c/v]

 (17)
> 

# use the expression of v to explicit the solution (obviously Z and Z bar are conjugate)

sol := eval(sol, [rel_4, rel_5]);

[Z = c/(-I*b+a), `#mrow(mover(mo(Z),mo("&#x305;")))` = c/(I*b+a)]

 (18)
> 

# transform the solution into a prettier form

sol := lhs~(sol) =~ numer~(rhs~(sol)) *~ conj~(denom~(rhs~(sol)))
                    /~
                    expand~(denom~(rhs~(sol)) *~ conj~(denom~(rhs~(sol))) )  assuming assumptions

[Z = -c*(-I*b-a)/(a^2+b^2), `#mrow(mover(mo(Z),mo("&#x305;")))` = c*(-I*b+a)/(a^2+b^2)]

 (19)
> 

# this complex number -represents the affix of the intersection point of the two lines

`#mrow(mo("&#8484;"))` := simplify(eval(Z, sol));  # asuming a^2+b^2 <> 0

c*(I*b+a)/(a^2+b^2)

 (20)
> 

# example

data := [a=3, b=-2, c=1, x=X, y=Y]:

with(plots):

f := eval(a*X + b*Y - c, data);
`#mrow(msub(mo(f),mo("&#8869;")))` := simplify(
                                        eval(
                                          eval(
                                            eval(`#mrow(mo("&#8459;"))`, [rel_1, rel_2, rel_7, rel_8]),
                                            [rel_4, rel_5]
                                          ),
                                          data
                                        )
                                      );
__Z := eval(`#mrow(mo("&#8484;"))`, data);

display(
  implicitplot(f, X=-1..1, Y=-1..1, color=blue),
  implicitplot(`#mrow(msub(mo(f),mo("&#8869;")))`, X=-1..1, Y=-1..1, color=red),
  pointplot([[Re(__Z), Im(__Z)]], symbol=circle, symbolsize=20),
  gridlines=true,
  scaling=constrained
);

3*X-2*Y-1

  

6*Y+4*X

  

3/13-(2/13)*I

  

 
> 

# just another example

data := [a=0, b=-3, c=-2, x=X, y=Y]:

with(plots):

f := eval(a*X + b*Y - c, data);
`#mrow(msub(mo(f),mo("&#8869;")))` := simplify(
                                        eval(
                                          eval(
                                            eval(`#mrow(mo("&#8459;"))`, [rel_1, rel_2, rel_7, rel_8]),
                                            [rel_4, rel_5]
                                          ),
                                          data
                                        )
                                      );
__Z := eval(`#mrow(mo("&#8484;"))`, data);

display(
  implicitplot(f, X=-1..1, Y=-1..1, color=blue),
  implicitplot(`#mrow(msub(mo(f),mo("&#8869;")))`, X=-1..1, Y=-1..1, color=red),
  pointplot([[Re(__Z), Im(__Z)]], symbol=circle, symbolsize=20),
  gridlines=true,
  scaling=constrained
);

-3*Y+2

  

6*X

  

(2/3)*I

  

 
> 

 


 

Download formal_solution_in_complex_plane.mw

 

Hi,  Have you try to use the geomet...

mmcdara 6650
June 01 2020
0 3

Hi, 
Have you try to use the geometry package?
Here are some of its capabilities that could interest you.


 

> 

restart:
> 

with(geometry):
with(RealDomain):
> 

local D, O:
> 

_EnvHorizontalName := x: _EnvVerticalName := y:
point(O, [0, 0]);

O

 (1)
> 

# D is assumed neither horizontal nor vertical

line(D, a*x+b*y+c) assuming a^2 + b^2 > 0, a<>0, b<>0:
PerpendicularLine(Perp, O, D):
projection(H, O, D):
coordinates(H);

[-a*c/(a^2+b^2), -b*c/(a^2+b^2)]

 (2)
> 

# D is assumed vertical

line(D, b*y+c) assuming b<>0, b^2 > 0:
PerpendicularLine(Perp, O, D):
projection(H, O, D):
coordinates(H);

[0, -c/b]

 (3)
> 

# D is assumed horizontal

line(D, a*x+c) assuming a<>0, a^2 > 0:
PerpendicularLine(Perp, O, D):
projection(H, O, D):
coordinates(H);

[-c/a, 0]

 (4)
> 

 


 

Download Pied_de_nez.mw

 

There are probably a lot of ways to that...

mmcdara 6650
June 01 2020
0 1

There are probably a lot of ways to do that.
Here is one of them
(the set of points I construct are in some sense the "shortest" for adding them any other point from the initial list obviously creates a new set of points whose sum id larger than the target sum)

 

> 

restart:
> 

randomize():
> 

M    := 10.:
roll := rand(0. .. M):
N    := 10:
L    := [seq(roll(), k=1..N)]

[3.829591035, 3.239717830, 8.629534219, 7.550798661, 7.992442854, 7.267252468, 2.715869968, 3.288246716, 3.345034572, 4.004192509]

 (1)
> 

ROLL := rand(0..N*M):
S    := ROLL();  # target sum

8.585001451

 (2)
> 

# First way

c := Statistics:-CumulativeSum(L):
if c[-1] < S then
  printf("The sum of all the numbers (%a) is less than rge desired sum (%a)\n", c[-1], S)
else
  select[flatten](`<`, c, S);
  selected_numbers := L[1..numelems(%)+1];

  # verify
  is(add(selected_numbers) > S)
end if;

Array([3.82959103500000, 7.06930886500000])

  

[3.829591035, 3.239717830, 8.629534219]

  

true

 (3)
> 

# Repeat this many (R) times.

R := 10:
for r from 1 to R do
  rL := combinat:-randperm(L):
  c := Statistics:-CumulativeSum(rL):
  if c[-1] < S then
    printf("The sum of all the numbers (%a) is less than rge desired sum (%a)\n", c[-1], S)
  else
    select[flatten](`<`, c, S);
    selected_numbers := rL[1..numelems(%)+1];
    printf("%a\n", selected_numbers);
  end if;
end do:

[2.715869968, 8.629534219]
[4.004192509, 7.267252468]
[8.629534219]
[3.829591035, 3.288246716, 8.629534219]
[2.715869968, 3.829591035, 7.267252468]
[8.629534219]
[7.550798661, 7.992442854]
[3.829591035, 7.267252468]
[7.267252468, 7.992442854]
[8.629534219]

  
> 

 

> 

 

> 

 

> 

 

> 

 

> 

# select the smallest number in wich is larger than 3.14

L[ListTools:-BinaryPlace(L, 3.14)+1];
min(select[flatten](`>`, L, 3.14))

4

  

4

 (4)
> 

# select the smallest number in wich is larger or equal than 4

min(select[flatten](`>`, L, 4));
min(select[flatten](`>=`, L, 4));

5

  

4

 (5)
> 

# all the points in L that are larger or equal than 4


LL := select[flatten](`>=`, L, 4)
 

[4, 5, 6, 7, 8, 9, 10]

 (6)
> 

# randomly select some of them, let's say 3 of them

combinat:-randperm(LL)[1..3]

[5, 4, 10]

 (7)
> 

 

> 

 

> 

 


 

Download search.mw

Without any constraints on the selcted numbers, all the solutions are equivalent.
But maybe you are looking for the selected numbers whose sum is the closest to the target sum?
Or the minimum set of numbers whose sum is larger than this target sum?

 

First 43 44 45 46 47 48 49 Last Page 45 of 57