Maple 2015 Questions and Posts

These are Posts and Questions associated with the product, Maple 2015

Hi, 

I often need to draw several curves on the same figure and the question of which colors to choose to distinguish them is crucial. 
In particular when these curves depend on a parameter I would like to color them by using a "regular" change of color.
The trick I use is this one
 

# here I want to plot N=10 different curves

N := 10:

# step 1: generate a list of colours using colorscheme
# (I was hoping for 10 levels of colors but I got 19)

plot(x, x=0..1, colorscheme=["Gold", "Blue"], numpoints=N):
MyColors := op([1,2, 2], %);
M := numelems(MyColors[..,1]);

# step 2: do the plots of interest, for instance

plots:-display(
  seq(
    plot(m*x^2, x=0..1, color=ColorTools:-Color([entries(MyColors[round(m/N*M)], nolist)])), 
    m=1..N
  )
);

My question is: How can we make that MyColors be a matrix whose the number of rows equal N (this to avoid the round(m/N*M) operation)?
(I will accept any other strategy)

Thanks in advance

A fascinating race is presently running (even if the latest results seem  to have put an end to it).
I'm talking of course about the US presidential elections.

My purpose is not to do politics but to discuss of a point of detail that really left me puzzled: the possibility of an electoral college tie.
I guess that this possibility seems as an aberration for a lot of people living in democratic countries. Just because almost everywhere at World electoral colleges contain an odd number of members to avoid such a situation!

So strange a situation that I did a few things to pass the time (of course with the earphones on the head so I don't miss a thing).
This is done with Maple 2015 and I believe that the amazing Iterator package (that I can't use thanks to the teleworking :-( ) could be used to do much more interesting things.

 

restart:

with(Statistics):

ElectoralCollege := Matrix(51, 2, [

Alabama,        9,        Kentucky,        8,        North_Dakota,        3,

Alaska,        3,        Louisiana,        8,        Ohio,        18,

Arizona,        11,        Maine,        4,        Oklahoma,        7,

Arkansas,        6,        Maryland,        10,        Oregon,        7,

California,        55,        Massachusetts,        11,        Pennsylvania,        20,

Colorado,        9,        Michigan,        16,        Rhode_Island,        4,

Connecticut,        7,        Minnesota,        10,        South_Carolina,        9,

Delaware,        3,        Mississippi,        6,        South_Dakota,        3,

District_of_Columbia,        3,        Missouri,        10,        Tennessee,        11,

Florida,        29,        Montana,        3,        Texas,        38,

Georgia,        16,        Nebraska,        5,        Utah,        6,

Hawaii,        4,        Nevada,        6,        Vermont,        3,

Idaho,        4,        New_Hampshire,        4,        Virginia,        13,

Illinois,        20,        New_Jersey,        14,        Washington,        12,

Indiana,        11,        New_Mexico,        5,        West_Virginia,        5,

Iowa,        6,        New_York,        29,        Wisconsin,        10,

Kansas,        6,        North_Carolina,        15,        Wyoming,        3
]):
 

ElectoralCollege := Vector(4, {(1) = ` 51 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(1)

add(ElectoralCollege[..,2]):
tie := %/2;

269

(2)

ec := convert(ElectoralCollege, listlist):

# Sets of states that form an electoral college tie

R      := 10^5:
nbties := 0:
states := NULL:
for r from 1 to R do
  poll  := combinat:-randperm(ec):
  cpoll := CumulativeSum(op~(2, poll)):
  if tie in cpoll then
    nbties := nbties+1;
    place  := ListTools:-Search(tie, cpoll);
    states := states, op~(1, poll)[1..place]:   # see below
  end if:
end do:

# electoral college tie is not so rare an event
# (prob of occurrence about 9.4 %).
#
# Why the hell the US constitution did not decide to have an odd
# number or electors to avoid ths kind of situation instead of
# introducing a complex mechanism when tie appears????

nbties;
evalf(nbties/R);

states := [states]:

9397

 

0.9397000000e-1

(3)

# What states participate to the tie?

names := sort(ElectoralCollege[..,1]):

all_states_in_ties := [op(op~(states))]:

howoften := Vector(
                    51,
                    i -> ListTools:-Occurrences(names[i], all_states_in_ties)
            ):

ScatterPlot(Vector(51, i->i), howoften);

 

# All the states seem to appear equally likely in an electoral college tie.
# Why? Does someone have a guess?
#
# The reason is obvious, as each state must appear in the basket of a candidate,
# then in case of a tie each state is either in op~(1, poll)[1..place] (candidate 1)
# or either in op~(1, poll)[place+1..51] (candidate 2);
# So, as we obtained 9397 ties, each states appears exactly 9397 times (with
# different occurences in the baskets of candidate 1 and 2).

 

# Lengths of the configurations that lead to a tie.
#
# Pleas refer to the answer above to understand why Histogram(lengths) should be
# symmetric.
lengths := map(i -> numelems(states[i]), [$1..nbties]):
sort(Tally(lengths))

[14 = 1, 15 = 2, 16 = 7, 17 = 36, 18 = 78, 19 = 179, 20 = 341, 21 = 507, 22 = 652, 23 = 849, 24 = 1015, 25 = 1041, 26 = 1056, 27 = 997, 28 = 862, 29 = 657, 30 = 515, 31 = 300, 32 = 158, 33 = 95, 34 = 41, 35 = 6, 36 = 2]

(4)

Histogram(lengths, range=min(lengths)..max(lengths), discrete=true)

 

ShortestConfigurations := map(i -> if lengths[i]=min(lengths) then states[i] end if, [$1..nbties]):
print~(ShortestConfigurations):

[New_York, Wisconsin, Illinois, Kentucky, Florida, New_Jersey, Mississippi, Indiana, Virginia, Maryland, California, Massachusetts, North_Carolina, Texas]

(5)

LargestConfigurations := map(i -> if lengths[i]=max(lengths) then states[i] end if, [$1..nbties]):
print~(LargestConfigurations):

[Alaska, Tennessee, North_Carolina, South_Carolina, District_of_Columbia, Colorado, Minnesota, Georgia, South_Dakota, New_Hampshire, Wyoming, Ohio, Rhode_Island, Arizona, Delaware, Montana, West_Virginia, Vermont, Michigan, Kentucky, Louisiana, Arkansas, Maine, Missouri, New_Mexico, Virginia, Maryland, Oregon, Wisconsin, Iowa, Kansas, Connecticut, North_Dakota, Nevada, Hawaii, Oklahoma]

 

[West_Virginia, Maryland, Massachusetts, Colorado, South_Dakota, Kentucky, Kansas, Wyoming, North_Dakota, Indiana, Michigan, Utah, Louisiana, Ohio, Alabama, Nebraska, Connecticut, Illinois, Oklahoma, Alaska, New_Jersey, District_of_Columbia, Oregon, Nevada, Missouri, Delaware, Washington, New_Hampshire, Arizona, Maine, South_Carolina, Hawaii, Vermont, Montana, Rhode_Island, Idaho]

(6)

# What could be the largest composition of a basket in case of a tie?
# (shortest composition is the complementary of the largest one)

ecs   := sort(ec, key=(x-> x[2]));
csecs := CumulativeSum(op~(2, ecs)):

# Where would the break locate?

tieloc := ListTools:-BinaryPlace(csecs, tie);

csecs[tieloc..tieloc+1]

[[North_Dakota, 3], [Alaska, 3], [Delaware, 3], [South_Dakota, 3], [District_of_Columbia, 3], [Montana, 3], [Vermont, 3], [Wyoming, 3], [Maine, 4], [Rhode_Island, 4], [Hawaii, 4], [Idaho, 4], [New_Hampshire, 4], [Nebraska, 5], [New_Mexico, 5], [West_Virginia, 5], [Arkansas, 6], [Mississippi, 6], [Utah, 6], [Nevada, 6], [Iowa, 6], [Kansas, 6], [Oklahoma, 7], [Oregon, 7], [Connecticut, 7], [Kentucky, 8], [Louisiana, 8], [Alabama, 9], [Colorado, 9], [South_Carolina, 9], [Maryland, 10], [Minnesota, 10], [Missouri, 10], [Wisconsin, 10], [Arizona, 11], [Massachusetts, 11], [Tennessee, 11], [Indiana, 11], [Washington, 12], [Virginia, 13], [New_Jersey, 14], [North_Carolina, 15], [Michigan, 16], [Georgia, 16], [Ohio, 18], [Pennsylvania, 20], [Illinois, 20], [Florida, 29], [New_York, 29], [Texas, 38], [California, 55]]

 

40

 

Array(%id = 18446744078888202358)

(7)

# This 40  states coniguration is not a tie.
#
# But list all the states in basket of candidate 1 and look to the 41th state (which is
# in the basket of candidate 2)

ecs[1..tieloc];
print():
ecs[tieloc+1]

[[North_Dakota, 3], [Alaska, 3], [Delaware, 3], [South_Dakota, 3], [District_of_Columbia, 3], [Montana, 3], [Vermont, 3], [Wyoming, 3], [Maine, 4], [Rhode_Island, 4], [Hawaii, 4], [Idaho, 4], [New_Hampshire, 4], [Nebraska, 5], [New_Mexico, 5], [West_Virginia, 5], [Arkansas, 6], [Mississippi, 6], [Utah, 6], [Nevada, 6], [Iowa, 6], [Kansas, 6], [Oklahoma, 7], [Oregon, 7], [Connecticut, 7], [Kentucky, 8], [Louisiana, 8], [Alabama, 9], [Colorado, 9], [South_Carolina, 9], [Maryland, 10], [Minnesota, 10], [Missouri, 10], [Wisconsin, 10], [Arizona, 11], [Massachusetts, 11], [Tennessee, 11], [Indiana, 11], [Washington, 12], [Virginia, 13]]

 

 

[New_Jersey, 14]

(8)

# It appears that exchanging Virginia and New_Jersey increases by 1 unit the college of candidate 1
# and produces a tie.

LargestBasketEver := [ ecs[1..tieloc-1][], ecs[tieloc+1] ];

add(op~(2, LargestBasketEver))

[[North_Dakota, 3], [Alaska, 3], [Delaware, 3], [South_Dakota, 3], [District_of_Columbia, 3], [Montana, 3], [Vermont, 3], [Wyoming, 3], [Maine, 4], [Rhode_Island, 4], [Hawaii, 4], [Idaho, 4], [New_Hampshire, 4], [Nebraska, 5], [New_Mexico, 5], [West_Virginia, 5], [Arkansas, 6], [Mississippi, 6], [Utah, 6], [Nevada, 6], [Iowa, 6], [Kansas, 6], [Oklahoma, 7], [Oregon, 7], [Connecticut, 7], [Kentucky, 8], [Louisiana, 8], [Alabama, 9], [Colorado, 9], [South_Carolina, 9], [Maryland, 10], [Minnesota, 10], [Missouri, 10], [Wisconsin, 10], [Arizona, 11], [Massachusetts, 11], [Tennessee, 11], [Indiana, 11], [Washington, 12], [New_Jersey, 14]]

 

269

(9)

# The largest electoral college tie contains 40 states (the shortest 11)


 

Download ElectoralCollegeTie.mw

Hi, 

I have an ODE system with a free parameter P, that I can only solve numerically.
To simplify let's suppose this parameter has only 2 values P1 and P2.
Let SOL1 the solution of this system when P=P1 and SOL2 its solution for P=P2:

Again for a sake of simplification let's suppose the independent variable is t and the dependent one is x(t).

How can I plot the difference x1(t)-x2(t) where x1(t) is the solution corresponding to SOL1 and x2(t) the solution corresponding to SOL2
(sorry I screwed something up and changed the config of my keyboard making the interrogation point unaccessible)

Here is a detailed example

Thanks in advance
 

restart:

with(plots):

# This is a 1D problem.
# A body with unit mass is submitted to an acceleration and to a solid-solid friction resistance
# The friction model is a simple Coulomb's

sys := {
         diff(x(t), t) = v(t),
         diff(v(t), t) = t - piecewise(v(t) > 0, f, v(t) < 0, -f, 0),
         x(0) = 0,
         v(0) = 0
       };

sol__0 := dsolve(sys, numeric, parameters=[f], maxfun=10^6)

sys := {diff(v(t), t) = t-piecewise(0 < v(t), f, v(t) < 0, -f, 0), diff(x(t), t) = v(t), v(0) = 0, x(0) = 0}

 

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := [f = f]; _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 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 1, (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) = 1000000, (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) = .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..3, {(1) = 0., (2) = 0., (3) = 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..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..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .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] = v(t), Y[2] = x(t)]`; YP[1] := X-piecewise(0 < Y[1], Y[3], Y[1] < 0, -Y[3], 0); YP[2] := Y[1]; 0 end proc, -1, 0, 0, 0, 0, 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] = v(t), Y[2] = x(t)]`; YP[1] := X-piecewise(0 < Y[1], Y[3], Y[1] < 0, -Y[3], 0); YP[2] := Y[1]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..3, {(1) = 0., (2) = 0., (3) = 0.}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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, v(t), x(t)], (4) = [f = f]}); _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 system doesn't seem to be solved formally

sol__exact := dsolve(sys, {v(t), x(t)})

# Attempt of solving it numerically.
#
# Let's observe that the "v" solution must be an increasing function of t
# so the friction model could be replaced by piecewise(v(t) > 0, f, 0) for v(t) >=0.
# But this is a particular case (think of a sine acceleration for instance) and
# I keep using the initial Coulomb's model.
#
# The main problem, as it appears below is the discontinuity of the friction effect at v(t)=0
# (thus at the initial time).

sol := dsolve(sys, numeric, parameters=[f], maxfun=10^6):
sol(parameters=[1]):
sol(0.1);

Error, (in sol) cannot evaluate the solution further right of 0.15837836e-1, maxfun limit exceeded (see ?dsolve,maxfun for details)

 

# Clearly increasing maxfum will help to simulate on a longer time but it's not the
# most efficient way to handle this problem.
#
# The classical way is to use a regularized Coulomb's model ; for instance

sys := {
             diff(x(t), t) = v(t),
             diff(v(t), t) = t - f *tanh(v(t)/c),
             x(0) = 0,
             v(0) = 0
           };

sol__1 := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8)

{diff(v(t), t) = t-f*tanh(v(t)/c), diff(x(t), t) = v(t), v(0) = 0, x(0) = 0}

 

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := [f = f, c = c]; _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 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 2, (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) = 100000000, (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) = .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..4, {(1) = 0., (2) = 0., (3) = Float(undefined), (4) = 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..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..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..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] = v(t), Y[2] = x(t)]`; YP[1] := X-Y[3]*tanh(Y[1]/Y[4]); YP[2] := Y[1]; 0 end proc, -1, 0, 0, 0, 0, 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] = v(t), Y[2] = x(t)]`; YP[1] := X-Y[3]*tanh(Y[1]/Y[4]); YP[2] := Y[1]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..4, {(1) = 0., (2) = 0., (3) = 0., (4) = 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); 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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, v(t), x(t)], (4) = [f = f, c = c]}); _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

(2)

# The smaller c the closer the solution to the one invoking the Coulomb's model.
#
# The question is: how can we measure as close to sol__exact we are?
# The idea is to choose a "small" value of c, lets say c__0, and to declare that the solution
# obtained with c__0 is "the exact solution".
#
# Studying the convergence of the solution as c__0 becomes smaller and smaller  will help
# choosing c__0.
# Of course the computational time is an increasing function of 1/c__0 : the smaller c__0
# the higher this time.

c__0 := 1e-8:  # arbitrary choice

sol__ref := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8):
sol__ref(parameters=[1, c__0]):

tmax := 10:

CodeTools:-Usage(sol__ref(tmax)):
numfun__0 := sol__ref(numfun);

odeplot(sol__ref, [t, v(t)], t=0..tmax)

memory used=33.20KiB, alloc change=0 bytes, cpu time=34.76s, real time=36.19s, gc time=0ns

 

50831166

 

 

# Assuming c__0 = 1e-8 is a value small enough to consider the solution plotted above
# is "the solution for the true Coulomb's model".
#
# Here are a few solutions for larger values of c:

cs   := 10^~[$-6..-1]:
tmax := 10:

printf("c = 1e-08  number of calls = %d\n", numfun__0):

for i from 1 to nops(cs) do
  sol__1(parameters=[1,cs[i]]):
  sol__1(tmax):
  couleur := ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]):
  Xplot__||i := odeplot(
                         sol__1, [t, x(t)], t=0..tmax,
                         legend=cs[i],
                         color=couleur
                       ):
  printf("c = %1.0e  number of calls = %d\n", cs[i], sol__1(numfun))
end do:


LY  := plottools:-transform((x, y) -> [x, log(y)]):
LXY := plottools:-transform((x, y) -> [log(x), log(y)]):

display(
         seq(Xplot__||i , i=1..nops(cs)),
         view=[0..3, 0..2], labels=[t, x], title="x(t) (lin-lin)"
       );

c = 1e-08  number of calls = 50831166

c = 1e-06  number of calls = 1308504

c = 1e-05  number of calls = 126848
c = 1e-04  number of calls = 12825
c = 1e-03  number of calls = 1513
c = 1e-02  number of calls = 456
c = 1e-01  number of calls = 246

 

 

# From this plot it seems that values of c less than 1e-4 give roughly the same curves.
#
# Remembering that the cimputational time is an increasing function of 1/c, we face the
# practical dilema of getting a solution "very close" to the ideal one (Coulomb's model)
# in a reasonable amount of time.
# Thus we need a small value of c nut not too small.
#
# The natural idea is to plot the difference between x(t) obtained with some value of c
# and x(t) obtained with our "reference value" c = c__0.
#
# How can we do this efficiently?
# (unfortunately piecewise and array/Array outputs are not allowed for parameterized solutions)
#
# For instance this is very slow


gap := proc(c, s1, s2, n)
  local s, h, k:
  sol__1(parameters=[1, c]):
  sol__1(s2):
  h := Matrix(n, 2):
  k := 1:
  for s from s1 to s2 by (s2-s1)/(n-1) do
    h[k, 1] := s:
    h[k, 2] := eval(x(t), sol__1(s)) - eval(x(t), sol__ref(s)):
    k       := k+1:
  end do:
  return h
end proc:


GAP := CodeTools:-Usage( gap(1e-5, 0, 1, 11) ):

Statistics:-ScatterPlot(GAP[..,1], GAP[..,2], symbol=solidcircle, color=blue)

memory used=399.35KiB, alloc change=0 bytes, cpu time=34.34s, real time=34.81s, gc time=0ns

 

 

 


 

Download odeplot_differences.mw


 

 

Hi everyone! Currently I'm studying magnetism, and I was thinking that maybe seeing represented the helix movement of an atom with a vector v parallel to the Magnetic Field B subject to the Force of Lorentz and the Hysteresis cycle created in the magnetazation of a ferromagnetic material could help me understand it more. I tried to create the plots in Maple 2015 but I couldn't.. anyone can help me by creating those two plots?

Dear Users! 

I want to evaluate following integral (f[i]'s) for different values of rho as shown. But find some problem and can't find f[3], f[4], f[5] and f[6]. If the integral from direct way is hard to calculate please guide me how I can evaluate it using numerical integration?

restart; HAM := [1, 2, 3, 4, 5, 6]; U := y^6*sin(x); alpha := 1/3;

for i while i <= nops(HAM) do

rho := op(i, HAM);

f[i] := int(int(U/((x^rho-p^rho)^alpha*(y^rho-q^rho)^alpha), q = 0 .. y), p = 0 .. x)

end do

Thanks in advance!

Special request to @acer @Carl Love @Kitonum @Preben Alsholm

Hi, 

I want to simulate a decision rule (let's say a production control strategy to fix the ideas) based on the outcomes of a binomial random variable.

One stage of this simulation involves operations which correspond to the NaiveProc described below. This procedure is a straightforward and rather naive translation in Maple of these operations.
It is given here to help you to understand what is to be done.

This coding being extremely slow, I implemented the required operations differently.
Given the fact that I'm only concerned by binomial samplings with outcomes 0 and 1, I basically replaced these samplings by three operations  U:=... , Got01:=... , Got0:= ... .
Each new implementation is about 200 times faster than the naive one.

My question is: Can we do still better?
For information the same simulation realized with R runs in less than 15s when REP = 10^8 (more than 100 times faster than what I'm capable to do with Maple)

Thanks for your help.


 

restart:

with(Statistics):

# This is the naive coding of the problem
#
# It has the advantage to explain clearly (at least I hope so) what is to be done.

NaiveProc := proc(REP, a, b, N)
  local roll, rep, p, K, N0, N1:
  roll := rand(a .. b);
  N0   := 0:
  N1   := 0:
  for rep from 1 to REP do
    p := roll();
    K := Sample(Binomial(N, p), 1)[1]:
    if K=0 then
      N0 := N0+1:
    elif K=1 then
      N1 := N1+1:
    end if:
  end do:
  N0, N1:
end proc:
    
CodeTools:-Usage(NaiveProc(10^3, 0, 5e-4, 25000));

memory used=78.92MiB, alloc change=68.00MiB, cpu time=1.21s, real time=6.30s, gc time=24.80ms

 

93, 80

(1)

# As the previous code is very slow I implemented several variants based on the fact
# that I care only on values of K equal to 0 and 1 and that it is then useless to sample
# a Binomial random variable to do this.
#
# Each of them is much faster than the naive coding... but can we even faster?
# For comparison, the same simulation realized in R takes less than 15s to run with REP=10^8

pi := (n, p, k) -> binomial(n, k) * (1-p)^(n-k) * p^k;


f1 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := Select(t -> is(t <= 0), U-P01):
  Got0  := Select(t -> is(t <= 0), U-P0 );
  N01   := numelems(Got01):
  N0    := numelems(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f1(10^5, 0, 5e-4, 25000));

print():

f2 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := select[flatten](`<=`, U-P01, 0):
  Got0  := select[flatten](`<=`, U-P0 , 0):
  N01   := numelems(Got01):
  N0    := numelems(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f2(10^5, 0, 5e-4, 25000));

print():

f3 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := Vector(REP, i -> `if`(U[i] <= P01[i], 1, 0)):
  Got0  := Vector(REP, i -> `if`(U[i] <= P0 [i], 1, 0)):
  N01   := add(Got01):
  N0    := add(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f3(10^5, 0, 5e-4, 25000));

proc (n, p, k) options operator, arrow; binomial(n, k)*(1-p)^(n-k)*p^k end proc

 

memory used=376.83MiB, alloc change=296.77MiB, cpu time=4.12s, real time=4.03s, gc time=430.38ms

 

8118, 7979

 

 

memory used=173.26MiB, alloc change=0 bytes, cpu time=2.30s, real time=1.98s, gc time=461.22ms

 

7993, 7793

 

 

memory used=154.19MiB, alloc change=0 bytes, cpu time=2.28s, real time=2.01s, gc time=432.40ms

 

7989, 7869

(2)



 

 


 

Download MakeItFaster.mw

Dear Users!

Hope you would be fine. I have some problem in execution the last loops (highlighted as red) where sumation is present. When NN>3 it takes alot of time more than 12 hours. Is there any alternative command to reduce the query. I am waiting for your response. Thanks in advance. 

restart; with(LinearAlgebra); Digits := 30; NN := 2; nu := 1; M1 := NN; M2 := NN; M3 := NN;

for k1 from 0 while k1 <= M1-1 do for k2 from 0 while k2 <= M2-1 do for k3 from 0 while k3 <= M3-1 do

SGP[M3*(M2*k1+k2)+k3+1] := simplify(sum((-1)^(k1-i1)*GAMMA(k1+i1+2*nu)*x^i1*(sum((-1)^(k2-i2)*GAMMA(k2+i2+2*nu)*y^i2*(sum((-1)^(k3-i3)*GAMMA(k3+i3+2*nu)*z^i3/(GAMMA(i3+nu+1/2)*factorial(k3-i3)*factorial(i3)), i3 = 0 .. k3))/(GAMMA(i2+nu+1/2)*factorial(k2-i2)*factorial(i2)), i2 = 0 .. k2))/(GAMMA(i1+nu+1/2)*factorial(k1-i1)*factorial(i1)), i1 = 0 .. k1)) end do end do end do;

SGPxyz := `<,>`(seq(seq(seq(SGP[M3*(M2*(i-1)+j-1)+k], k = 1 .. M3), j = 1 .. M2), i = 1 .. M1));

Lambda := `<,>`(seq(seq(seq(chi[M3*(M2*(i-1)+j-1)+k], k = 1 .. M3), j = 1 .. M2), i = 1 .. M1));

for i while i <= NN^3 do for j while j <= NN^3 do for k while k <= NN^3 do

q[i, j, k] := int(int(int(SGP[i]*SGP[j]*SGP[k]*(-x^2+x)^(nu-1/2)*(-y^2+y)^(nu-1/2)*(-z^2+z)^(nu-1/2), z = 0 .. 1), y = 0 .. 1), x = 0 .. 1) end do end do end do;

U := Matrix(NN^3, NN^3, 0);

for j while j <= NN^3 do for k while k <= NN^3 do U[j, k] := simplify(sum(chi[i1]*q[i1, j, k], i1 = 1 .. NN^3)) end do end do;

F := simplify(evalm(U));

Special request to @acer @Carl Love @Kitonum @Preben Alsholm

Hi,

I do not understand why this simple procedure evaluates so differently depending on the type of its second parameter?

A typo somewhere or a bug?


 

restart:

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

KL := (a, b) -> (1/4)*(2*ln(a+b)*a^2+4*ln(a+b)*b*a+2*ln(a+b)*b^2-2*ln(b)*b^2-a^2-2*a*b)/a

proc (a, b) options operator, arrow; (1/4)*(2*ln(a+b)*a^2+4*ln(a+b)*b*a+2*ln(a+b)*b^2-2*ln(b)*b^2-a^2-2*a*b)/a end proc

(2)

evalf(KL(1e-10, 1/2))

-.2

(3)

evalf(KL(1e-10, 0.5))

-0.2500000000e-10

(4)

evalf(KL(1e-10, convert(0.5, rational)))

-.2

(5)

limit(KL(a, 1/2), a=0, left);
limit(KL(a, 1/2), a=0, right)

-(1/2)*ln(2)

 

-(1/2)*ln(2)

(6)

limit(KL(a, 0.5), a=0, left);
limit(KL(a, 0.5), a=0, right)

Float(-infinity)

 

Float(infinity)

(7)

 


 

Download I_am_lost.mw

I have an integral of a sum

Int(Sum(-(Nb*t-n*t__rev)*exp((1/2)*(L+sqrt(-4*C*L*R^2+L^2))*t/(L*R*C)-(1/2)*(Nb*t-n*t__rev)^2/(Nb^2*sigma__b^2)), n = 0 .. Nb-1), t);

that I want to convert into a sum of integrals since Maple will not integrate the sum but it will integrate the summands. I tried IntegrationTools, but none of these tools seem to do it. It is a part of a larger expression so doing this by hand is not a real option.

Mac Dude.

Hi!

The algorithm in this PDF is written in Fortran, but unfortunately I do not known this programming languaje (actually, I am not an "expert" in progamation).

Algoritmo_Fortram.pdf  [removed by moderator. © Institute of Mathematics AS CR, 1980]

 

https://dml.cz/handle/10338.dmlcz/103855

Could someone please write this algorithm in Maple? Or, at least, indicate me how to do.

Many thanks in advance for your comments.

Hello Everyone
My maple is evaluating the following provided in figure. I have  attached MAPLE file. Kindly help me to solve this. Thanks in advance.

question.mw

Hi, 

While "numpoints=..." is said to be an admissible option of SurfaceOfRevolution (or VolumeOfRevolution, default value set to 50), the graphical results doesn't account for the value of numpoints.

Is this a known issue?

TIA

Surface0fRevolution.mw



 

Hi, 

I would like to share this work I've done. 
No big math here, just a demonstrator of Maple's capabilities in 3D visualization.

All the plots in the file have been discarded to reduce the size of this post. Here is a screen capture to give you an idea of what is inside the file.

Download 3D_Visualization.mw

Hi, 

I ran this command
c := NeutralSpread(Color("RGB", "Yellow"), 10);
and I wonder how I could get the  NearestNameColor of all elements in list  c ?


Thanks in advance

Hello

I wanted to create a group of interactive components that, when given an eqaution in a math container of an ellepses (in the form of ((x/A)^2)+((y/B)^2)=1) would pick A and B and give back two formulas of the type

x(t)=A*cos(t)

y(t)=b*sin(t)

So, is there a way to extract A and B from the initial formula?

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