acer

32343 Reputation

29 Badges

19 years, 327 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

It seems that the solution obtained at lower values of n can provide a useful initial point for a higher value of n.

Below it's done in two ways.

The first uses this bootstrapping approach for values of n=7,13,17,18,19,20.

The second does it in a loop, using all n from 2 to 20. (slower overall, but perhaps more robust!?)

restart

interface(rtablesize = 10)

kernelopts(version)

`Maple 2019.0, X86 64 LINUX, Mar 9 2019, Build ID 1384062`

r[0] := 3.12*10^(-5)

r[1] := 1.00*10^(-4)

r[2] := 1.25*10^(-4)

r[3] := 1.87*10^(-4)

r[4] := 2.50*10^(-4)

e1 := t[0] = (-p[0, fail]^7+1)/((1-p[0, fail])*((-p[0, fail]^7+1)/(1-p[0, fail])+(17/2+(17/2)*p[0, fail]+(33/2)*p[0, fail]^2+(33/2)*p[0, fail]^3+(65/2)*p[0, fail]^4+(65/2)*p[0, fail]^5+(65/2)*p[0, fail]^6)/p[0, idle]+(1-r[0])*(p[0, fail]^6+(-p[0, fail]^6+1)*p[0, succ]/(1-p[0, fail]))/r[0])); e2 := t[1] = (-p[1, fail]^7+1)/((1-p[1, fail])*((-p[1, fail]^7+1)/(1-p[1, fail])+(9/2+(9/2)*p[1, fail]+(17/2)*p[1, fail]^2+(17/2)*p[1, fail]^3+(33/2)*p[1, fail]^4+(33/2)*p[1, fail]^5+(33/2)*p[1, fail]^6)/p[1, idle]+(1-r[1])*(p[1, fail]^6+(-p[1, fail]^6+1)*p[1, succ]/(1-p[1, fail]))/r[1])); e3 := t[2] = (-p[2, fail]^7+1)/((1-p[2, fail])*((-p[2, fail]^7+1)/(1-p[2, fail])+(5/2+(5/2)*p[2, fail]+(9/2)*p[2, fail]^2+(9/2)*p[2, fail]^3+(17/2)*p[2, fail]^4+(17/2)*p[2, fail]^5+(17/2)*p[2, fail]^6)/p[2, idle]+(1-r[2])*(p[2, fail]^6+(-p[2, fail]^6+1)*p[2, succ]/(1-p[2, fail]))/r[2])); e4 := t[3] = (-p[3, fail]^7+1)/((1-p[3, fail])*((-p[3, fail]^7+1)/(1-p[3, fail])+(3/2+(3/2)*p[3, fail]+(5/2)*p[3, fail]^2+(5/2)*p[3, fail]^3+(9/2)*p[3, fail]^4+(9/2)*p[3, fail]^5+(9/2)*p[3, fail]^6)/p[3, idle]+(1-r[3])*(p[3, fail]^6+(-p[3, fail]^6+1)*p[3, succ]/(1-p[3, fail]))/r[3])); e5 := t[4] = (-p[4, fail]^7+1)/((1-p[4, fail])*((-p[4, fail]^7+1)/(1-p[4, fail])+(1+p[4, fail]+(3/2)*p[4, fail]^2+(3/2)*p[4, fail]^3+(5/2)*p[4, fail]^4+(5/2)*p[4, fail]^5+(5/2)*p[4, fail]^6)/p[4, idle]+(1-r[4])*(p[4, fail]^6+(-p[4, fail]^6+1)*p[4, succ]/(1-p[4, fail]))/r[4])); e6 := p[0, fail] = n*t[0]*(1-(1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n)+0.1590000000e-1*n*t[0]*(1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e7 := p[1, fail] = n*t[1]*(1-(1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n)+0.1590000000e-1*n*t[1]*(1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e8 := p[2, fail] = n*t[2]*(1-(1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n*(1-t[3])^n*(1-t[4])^n)+0.1590000000e-1*n*t[2]*(1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e9 := p[3, fail] = (1/2)*n*t[3]*(1-(1-t[3])^(n-1)*(1-t[4])^n)+(1/2)*n*t[3]*(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n)+(0.1590000000e-1*((1/2)*n*t[3]*(1-t[3])^(n-1)*(1-t[4])^n+(1/2)*n*t[3]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e10 := p[4, fail] = (1/2)*n*t[4]*(1-(1-t[4])^(n-1)*(1-t[3])^n)+(1/2)*n*t[4]*(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n)+(0.1590000000e-1*((1/2)*n*t[4]*(1-t[4])^(n-1)*(1-t[3])^n+(1/2)*n*t[4]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e11 := p[0, succ] = .9841000000*n*t[0]*(1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e12 := p[1, succ] = .9841000000*n*t[1]*(1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e13 := p[2, succ] = .9841000000*n*t[2]*(1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n*(1-t[3])^n*(1-t[4])^n/(1-(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e14 := p[3, succ] = (.9841000000*((1/2)*n*t[3]*(1-t[3])^(n-1)*(1-t[4])^n+(1/2)*n*t[3]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e15 := p[4, succ] = (.9841000000*((1/2)*n*t[4]*(1-t[4])^(n-1)*(1-t[3])^n+(1/2)*n*t[4]*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n))/(1-(1/2)*(1-t[3])^n*(1-t[4])^n-(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^n*(1-t[4])^n); e16 := p[0, idle] = (1-t[0])^(n-1)*(1-t[1])^n*(1-t[2])^n; e17 := p[1, idle] = (1-t[1])^(n-1)*(1-t[0])^n*(1-t[2])^n; e18 := p[2, idle] = (1-t[2])^(n-1)*(1-t[1])^n*(1-t[0])^n; e19 := p[3, idle] = (1/2)*(1-t[3])^(n-1)*(1-t[4])^n+(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[3])^(n-1)*(1-t[4])^n; e20 := p[4, idle] = (1/2)*(1-t[4])^(n-1)*(1-t[3])^n+(1/2)*(1-t[0])^n*(1-t[1])^n*(1-t[2])^n*(1-t[4])^(n-1)*(1-t[3])^n

sys := [seq(eval(cat(e, i)), i = 1 .. 20)]

conds:=map(`=`,indets(sys,And(name,Non(constant))) minus {n},0..1);

{p[0, fail] = 0 .. 1, p[0, idle] = 0 .. 1, p[0, succ] = 0 .. 1, p[1, fail] = 0 .. 1, p[1, idle] = 0 .. 1, p[1, succ] = 0 .. 1, p[2, fail] = 0 .. 1, p[2, idle] = 0 .. 1, p[2, succ] = 0 .. 1, p[3, fail] = 0 .. 1, p[3, idle] = 0 .. 1, p[3, succ] = 0 .. 1, p[4, fail] = 0 .. 1, p[4, idle] = 0 .. 1, p[4, succ] = 0 .. 1, t[0] = 0 .. 1, t[1] = 0 .. 1, t[2] = 0 .. 1, t[3] = 0 .. 1, t[4] = 0 .. 1}

sol[7] := fsolve(eval(sys, n = 7), conds);

{p[0, fail] = 0.1725732060e-2, p[0, idle] = .9903415166, p[0, succ] = .1043549027, p[1, fail] = 0.3088989313e-2, p[1, idle] = .9905749067, p[1, succ] = .1868442149, p[2, fail] = 0.3454963251e-2, p[2, idle] = .9906375843, p[2, succ] = .2089969336, p[3, fail] = 0.4860178835e-2, p[3, idle] = .9862434767, p[3, succ] = .2970648958, p[4, fail] = 0.5619842175e-2, p[4, idle] = .9863404432, p[4, succ] = .3435277781, t[0] = 0.2981348398e-3, t[1] = 0.5336754095e-3, t[2] = 0.5969115356e-3, t[3] = 0.6286120788e-3, t[4] = 0.7268596334e-3}

sol[13] := fsolve(eval(sys, n = 13), sol[7], conds);

{p[0, fail] = 0.1811745261e-2, p[0, idle] = .9817474126, p[0, succ] = .1035504405, p[1, fail] = 0.3242216921e-2, p[1, idle] = .9819805650, p[1, succ] = .1854040568, p[2, fail] = 0.3626151289e-2, p[2, idle] = .9820431831, p[2, succ] = .2073876078, p[3, fail] = 0.4984792186e-2, p[3, idle] = .9739267569, p[3, succ] = .2952561294, p[4, fail] = 0.5763517604e-2, p[4, idle] = .9740230984, p[4, succ] = .3414368147, t[0] = 0.3004375377e-3, t[1] = 0.5377969671e-3, t[2] = 0.6015257528e-3, t[3] = 0.6324491804e-3, t[4] = 0.7312975182e-3}

sol[17] := fsolve(eval(sys, n = 17), sol[13], conds);

{p[0, fail] = 0.1902958391e-2, p[0, idle] = .9759851872, p[0, succ] = .1030103109, p[1, fail] = 0.3405000403e-2, p[1, idle] = .9762181798, p[1, succ] = .1844371111, p[2, fail] = 0.3808095785e-2, p[2, idle] = .9762807580, p[2, succ] = .2063070876, p[3, fail] = 0.5120884499e-2, p[3, idle] = .9657202245, p[3, succ] = .2940419271, p[4, fail] = 0.5920604827e-2, p[4, idle] = .9658161483, p[4, succ] = .3400331814, t[0] = 0.3020036099e-3, t[1] = 0.5406000512e-3, t[2] = 0.6046639560e-3, t[3] = 0.6350512932e-3, t[4] = 0.7343070663e-3}

sol[18] := fsolve(eval(sys, n = 18), sol[17], conds);

{p[0, fail] = 0.1930046058e-2, p[0, idle] = .9745404430, p[0, succ] = .1028747923, p[1, fail] = 0.3453365646e-2, p[1, idle] = .9747733954, p[1, succ] = .1841945044, p[2, fail] = 0.3862160540e-2, p[2, idle] = .9748359636, p[2, succ] = .2060359860, p[3, fail] = 0.5161607349e-2, p[3, idle] = .9636691261, p[3, succ] = .2937373132, p[4, fail] = 0.5967623643e-2, p[4, idle] = .9637649452, p[4, succ] = .3396810444, t[0] = 0.3023990953e-3, t[1] = 0.5413079232e-3, t[2] = 0.6054564644e-3, t[3] = 0.6357074542e-3, t[4] = 0.7350659702e-3}

sol[19] := fsolve(eval(sys, n = 19), sol[18], conds);

{p[0, fail] = 0.1958858820e-2, p[0, idle] = .9730940045, p[0, succ] = .1027390773, p[1, fail] = 0.3504818935e-2, p[1, idle] = .9733269167, p[1, succ] = .1839515461, p[2, fail] = 0.3919679479e-2, p[2, idle] = .9733894748, p[2, succ] = .2057644920, p[3, fail] = 0.5205026377e-2, p[3, idle] = .9616182321, p[3, succ] = .2934322698, p[4, fail] = 0.6017760081e-2, p[4, idle] = .9617139466, p[4, succ] = .3393284112, t[0] = 0.3027961900e-3, t[1] = 0.5420186757e-3, t[2] = 0.6062522009e-3, t[3] = 0.6363658969e-3, t[4] = 0.7358275140e-3}

sol[20] := fsolve(eval(sys, n = 20), sol[19], conds);

{p[0, fail] = 0.1989402444e-2, p[0, idle] = .9716458621, p[0, succ] = .1026031649, p[1, fail] = 0.3559370593e-2, p[1, idle] = .9718787340, p[1, succ] = .1837082346, p[2, fail] = 0.3980664162e-2, p[2, idle] = .9719412822, p[2, succ] = .2054926035, p[3, fail] = 0.5251149844e-2, p[3, idle] = .9595675377, p[3, succ] = .2931267948, p[4, fail] = 0.6071023693e-2, p[4, idle] = .9596631475, p[4, succ] = .3389752792, t[0] = 0.3031949056e-3, t[1] = 0.5427323292e-3, t[2] = 0.6070511885e-3, t[3] = 0.6370266357e-3, t[4] = 0.7365917143e-3}

eval( eval(map(lhs-rhs,sys), n=20), sol[20] );

[0., 0., 0., -0.2e-12, 0.3e-12, -0.35e-10, -0.65e-10, -0.75e-10, -0.129e-9, -0.148e-9, -0.29e-8, -0.49e-8, -0.57e-8, -0.88e-8, -0.103e-7, -0.6e-9, -0.5e-9, -0.4e-9, -0.11e-8, -0.12e-8]

#
# Automating it, but incrementing n by only 1 each iteration.
#
S[2] := fsolve(eval(sys, n = 2), conds):
for i from 3 to 20 do
  st := time();
  S[i] := fsolve(eval(sys, n = i), S[i-1], conds);
  elapsed := time()-st;
  if type(eval(S[i],1), 'specfunc(fsolve)') then
    print("failed at i=%1",i);
  else
    print(sprintf("succeeded at i=%a in %a seconds",i,elapsed));
  end if;
end do:

"succeeded at i=3 in .72e-1 seconds"

"succeeded at i=4 in .81e-1 seconds"

"succeeded at i=5 in .100 seconds"

"succeeded at i=6 in .141 seconds"

"succeeded at i=7 in .234 seconds"

"succeeded at i=8 in 1.525 seconds"

"succeeded at i=9 in .537 seconds"

"succeeded at i=10 in .839 seconds"

"succeeded at i=11 in 1.271 seconds"

"succeeded at i=12 in 1.867 seconds"

"succeeded at i=13 in 2.681 seconds"

"succeeded at i=14 in 3.874 seconds"

"succeeded at i=15 in 5.707 seconds"

"succeeded at i=16 in 9.394 seconds"

"succeeded at i=17 in 9.906 seconds"

"succeeded at i=18 in 17.662 seconds"

"succeeded at i=19 in 23.167 seconds"

"succeeded at i=20 in 36.280 seconds"

S[20];

{p[0, fail] = 0.1989402444e-2, p[0, idle] = .9716458621, p[0, succ] = .1026031649, p[1, fail] = 0.3559370593e-2, p[1, idle] = .9718787340, p[1, succ] = .1837082346, p[2, fail] = 0.3980664162e-2, p[2, idle] = .9719412822, p[2, succ] = .2054926035, p[3, fail] = 0.5251149844e-2, p[3, idle] = .9595675377, p[3, succ] = .2931267948, p[4, fail] = 0.6071023693e-2, p[4, idle] = .9596631475, p[4, succ] = .3389752792, t[0] = 0.3031949056e-3, t[1] = 0.5427323292e-3, t[2] = 0.6070511885e-3, t[3] = 0.6370266357e-3, t[4] = 0.7365917143e-3}

eval( eval(map(lhs-rhs,sys), n=20), S[20] );

[0., 0., 0., -0.2e-12, 0.3e-12, -0.35e-10, -0.65e-10, -0.75e-10, -0.129e-9, -0.148e-9, -0.29e-8, -0.49e-8, -0.57e-8, -0.88e-8, -0.103e-7, -0.6e-9, -0.5e-9, -0.4e-9, -0.11e-8, -0.12e-8]

 

 

Download solveSYS_ac.mw

If we plot the values for each variable, in the solution as n changes, we might observe that all the curves appear linear except for p[j,fail] , j=1..4, which are monotomic and could be interpolated. I would guess that reasonably close initial points could be generated for the higher n, based on a smaller number of the solutions for lower n (which compute more quickly).
solveSYS_ac2.mw
I would not be surprised if symbolic examination of the equations could lead to a non-numeric or hybrid symbolic-numeric approach where the solutions (dependent on n) were constructed very quickly. I don't have time for that.

restart;

Digits := 12;

12

ee := (-69/50)*sin((32/25)*x) = (4/5)*exp((-23/50)*x):
 

rts := Student:-Calculus1:-Roots(ee, x=0..10, numeric);

[2.59251746013, 4.86028632939, 7.37831403654, 9.81251449999]

evalf[10](rts[1]);

2.592517460

# another way
RootFinding:-NextZero(unapply((lhs-rhs)(ee),x), 0.0);
evalf[10](%);

2.59251746013

2.592517460

Download nz.mw

It turns out that increasing Digits up front may not be necessary here. But I include it anyway, since it may give you an idea in the future.

The help page for parse describes that the offset option can be supplied alongside the statement option, and that the lastread option allows for a running update to a specified name to hold the position least read.

Doing that in a loop allows for multiple statements to be parsed from a single string. Or you can construct a procedure to execute such a loop.

That is more robust than simply splitting the string at colons and/or semicolons, since those could appear in locations where the parsing should not be split.

An extra choice is to try and discern (and return) information about what (previously unassigned) names might have been newly assigned by this process. I realize that this may be of less relevance to your purpose, since your script that constructs the strings likely already knows to what names it will assign.

restart;

#
# This technique works even if the statement separators
# include a mix of colon and semicolon.
#
# It also works if colon or semicolon appears in places other
# than as statement separators.
#
# Since this use of `anames(user)` only detect assignments to
# names which were not yet assigned you may choose to have
# the procedure return NULL and get rid of the anames stuff.
#

parseall:=proc(s::string)
  local res, n, A, temp;
  A := {anames(':-user')};
  n:=0; res := [];
  while n<length(s) do
    temp := parse(s,':-lastread'='n',':-offset'=n,':-statement');
    res := [res[],[{anames(':-user')} minus A,eval(temp,1)]];
    A := {anames(':-user')};
  end do;
  res;
end proc:

 

data1 := "x:='x';y:='y';sol:=dsolve(diff(y(x),x)=1,y(x));":

out := parseall(data1);

[[{}, x], [{}, y], [{sol}, y(x) = x+_C1]]

x,y,sol;

x, y, y(x) = x+_C1

# What names were just newly assigned?
map[2](op,1,eval(out,1));

[{}, {}, {sol}]

data2 := "x:='x':f:=proc(a::integer) a^2; end proc;y:='y':`g;h`;sol:=dsolve(diff(y(x),x)=1,y(x));":

out := parseall(data2);

[[{}, x], [{f}, proc (a::integer) a^2 end proc], [{}, y], [{}, `g;h`], [{}, y(x) = x+_C1]]

x,y,sol;

x, y, y(x) = x+_C1

# What names were just newly assigned?
# Notice that `sol` was already assigned, so didn't show up.
map[2](op,1,eval(out,1));

[{}, {f}, {}, {}, {}]

sol := 'sol': f := 'f':
out := parseall(data2);

[[{}, x], [{f}, proc (a::integer) a^2 end proc], [{}, y], [{}, `g;h`], [{sol}, y(x) = x+_C1]]

# What names were just newly assigned?
map[2](op,1,eval(out,1));

[{}, {f}, {}, {}, {sol}]

# This first example succeeds, by splitting the string.
# There is no indication of what was assigned.
#
eval~(parse~(StringTools:-Split(data1, ";")));

[x, y, y(x) = x+_C1]

# This example fails, by splitting the string. It does not
# correctly interpret colons and semicolons used within
# procedures as statement separators. It also does not
# correctly interpret colons and semicolons appearing other
# than as statement separators.
#
eval~(parse~(StringTools:-Split(data2, ";")));

Warning, extra characters at end of parsed string

Error, incorrect syntax in parse: reserved word `end` unexpected (near 4th character of parsed string)

map(lprint, StringTools:-Split(data2, ";")):

"x:='x':f:=proc(a::integer) a^2"
" end proc"
"y:='y':`g"
"h`"
"sol:=dsolve(diff(y(x),x)=1,y(x))"
""

Download parsestatements.mw

By "newly assigned" I mean unassigned before the call. I don't mean re-assigned.

The error message describes the situation: the Int call is not supported by evalhf.

Did you just want to make it faster? Note that the quadrature scheme can utilize evalhf for the evaluations of the integrand.

restart; f1 := proc (XX::Array) local GG, i; GG := 0; for i to 500 do GG := GG+Int(z^2/(1000*XX[i]-z)^2, z = 0 .. 100) end do; return GG end proc; A1 := Array(1 .. 500, [seq(-i^2+100-1, i = 1 .. 500)]); CodeTools:-Usage(evalf(f1(A1)))

memory used=34.63MiB, alloc change=6.01MiB, cpu time=180.00ms, real time=188.00ms, gc time=12.68ms

.2915102677

restart; f1 := proc (XX::Array) local i; add(evalf(Int(z^2/(1000*XX[i]-z)^2, z = 0 .. 100, method = _d01ajc)), i = 1 .. 500) end proc; A1 := Array(1 .. 500, [seq(-i^2+100-1, i = 1 .. 500)]); CodeTools:-Usage(f1(A1))

memory used=13.34MiB, alloc change=2.00MiB, cpu time=73.00ms, real time=81.00ms, gc time=0ns

.2915102681

infolevel[`evalf/int`] := 2; forget(evalf); evalf(Int(z^2/(1000*A1[2]-z)^2, z = 0 .. 100, method = _d01ajc)); infolevel[`evalf/int`] := 0

2

evalf/int/control: integrating on 0 .. 100 the integrand

z^2/(95000-z)^2

Control: Entering NAGInt
trying d01ajc (nag_1d_quad_gen)
d01ajc: trying evalhf callbacks
d01ajc: result=.369928326560882854e-4
d01ajc: abserr=.205351472804473049e-18; num_subint=1; fun_count=21
result=.369928326560882854e-4

0.3699283266e-4

0

 

Download evalhf_int_ac.mw

Or do you have other computations to go inside procedure f1, which you also want done under evalhf?

The subs command does syntactic substituion without evaluation (unless the command name is indexed by the key eval).

The eval command does a more mathematical substitution, in the sense that it has additional knowledge of the mathematics germane to some substitutions.

Often the two produce the same result. Sometimes one specifically wants to avoid evaluation prior to being able to do subsequent manipulations of the result.

Sometimes subs is a bit faster. Relying on that only makes sense up to the situations where one doesn't need the mathematical knowledge.

subs(x=1.2, sin(x));
                         sin(1.2)

subs[eval](x=1.2, sin(x));
                       0.9320390860

eval(sin(x), x=1.2);
                       0.9320390860

There are also differences in how substitutions are done. The subs command allows you to do them together or in sequence. Eg,

subs([a=b,b=c], [a,b]);
                          [b, c]

subs(a=b,b=c, [a,b]);
                          [c, c]

subs(b=c,a=b, [a,b]);
                          [b, c]

eval([a,b], [a=b,b=c]);
                          [b, c]

I'm not sure which aspects are most relevant to your query. There are other subtleties or examples which might interest you.

I ran the code below on an Intel i7-7700 3.60GHz which has 4 physical cores with hyperthreading (8 threads or logical cores). Since Maple doesn't recognize this as having "new" improved hyperthreading it sets kernelopts(numcpus) to 4 by default. So in some of the examples below I set kernelopts(numcpus) to 7 and kernelopts(gcmaxthreads) to 7 or 8 (which seems optimal on this host for some examples).

I also ran it on an Intel i5-7400 3.00GHz which has 4 physical cores without hyperthreading (4 threads or logical cores). The relative timing behavior amongst the example was generally the same, except timings were not as good. I left kernelopts(numcpus) and kernelopts(gcmaxthreads) alone, at their default value of 4, since raising those resulted in worse performance.

I did all that using 64bit Maple 2019.0 on Linux.

I observed that almost half of the total computation time is spent in memory management (garbage collection). I found that by reducing the setting of kernelopts(gcthreadmemorysize) from its default I could improve the total timing. This benefit seemed to happen across other setting variants.

I made sure that the machines were not otherwise significantly loaded, since the timing measurements are in real (wall clock) time. I made sure that for size nelems=10^7 the machines had at least 2.5GB of free RAM and so would not swap out the computation.

One interesting result is that simple use of Threads:-Seq attained almost as good performance as the more complicated Threads:-Task set up, under the tweaked kernelopts settings.

For the Threads:-Task set up I tried to ensure that few temporary Arrays were created (ie. without concatenation). The exception to that was that I did not find improvement for replacing seq calls with inplace map calls or straight Array construction with initializer.

restart;
'numcpus'=kernelopts(numcpus), 'gcmaxthreads'=kernelopts(gcmaxthreads),
'gcthreadmemorysize'=kernelopts(gcthreadmemorysize);

nelems := 10000000:
n := 374894756873546859847556:

str,sgctr := time[real](), kernelopts(gctotaltime[real]):
A := Array(1 .. 4, 1 .. nelems):
A[1,1..nelems] := Array([seq( i^10*n, i=1..nelems)]):
A[2,1..nelems] := Array([seq( length(A[1,i])-1, i=1..nelems)]):
A[3,1..nelems] := Array([seq(iquo(A[1,i], 10^(A[2,i]-2)),i=1..nelems)]):
A[4,1..nelems] := Array([seq(irem(A[1,i],1000),i=1..nelems)]):
RT,GCRT := time[real]()-str, kernelopts(gctotaltime[real])-sgctr:
print(sprintf("total real time: %a secs   gc real time: %a secs",
              evalf[4](RT), evalf[4](GCRT)));

numcpus = 4, gcmaxthreads = numcpus, gcthreadmemorysize = 67108864

"total real time: 53.82 secs   gc real time: 29.47 secs"

restart;
kernelopts(numcpus=7): kernelopts(gcmaxthreads=7): kernelopts(gcthreadmemorysize=2^17):
'numcpus'=kernelopts(numcpus), 'gcmaxthreads'=kernelopts(gcmaxthreads),
'gcthreadmemorysize'=kernelopts(gcthreadmemorysize);

nelems := 10000000:
n := 374894756873546859847556:

str,sgctr := time[real](), kernelopts(gctotaltime[real]):
A := Array(1 .. 4, 1 .. nelems):
A[1,1..nelems] := Array([seq(i^10*n, i=1..nelems)]):
A[2,1..nelems] := Array([seq(length(A[1,i])-1, i=1..nelems)]):
A[3,1..nelems] := Array([seq(iquo(A[1,i], 10^(A[2,i]-2)),i=1..nelems)]):
A[4,1..nelems] := Array([seq(irem(A[1,i],1000),i=1..nelems)]):
RT,GCRT := time[real]()-str, kernelopts(gctotaltime[real])-sgctr:
print(sprintf("total real time: %a secs   gc real time: %a secs",
              evalf[4](RT), evalf[4](GCRT)));

numcpus = 7, gcmaxthreads = 7, gcthreadmemorysize = 131072

"total real time: 44.62 secs   gc real time: 20.62 secs"

restart;
kernelopts(numcpus=7): kernelopts(gcmaxthreads=7): kernelopts(gcthreadmemorysize=2^17):
'numcpus'=kernelopts(numcpus), 'gcmaxthreads'=kernelopts(gcmaxthreads),
'gcthreadmemorysize'=kernelopts(gcthreadmemorysize);

nelems := 10000000:
n := 374894756873546859847556:

str,sgctr := time[real](), kernelopts(gctotaltime[real]):
A := Array(1 .. 4, 1 .. nelems):
A[1,1..nelems] := Array([Threads:-Seq(i^10*n, i=1..nelems)]):
A[2,1..nelems] := Array([Threads:-Seq(ilog10(A[1,i]), i=1..nelems)]):
A[3,1..nelems] := Array([Threads:-Seq(iquo(A[1,i], 10^(A[2,i]-2)),i=1..nelems)]):
A[4,1..nelems] := Array([Threads:-Seq(irem(A[1,i],1000),i=1..nelems)]):
RT,GCRT := time[real]()-str, kernelopts(gctotaltime[real])-sgctr:
print(sprintf("total real time: %a secs   gc real time: %a secs",
              evalf[4](RT), evalf[4](GCRT)));

numcpus = 7, gcmaxthreads = 7, gcthreadmemorysize = 131072

"total real time: 18.00 secs   gc real time: 9.000 secs"

restart;
kernelopts(numcpus=7): kernelopts(gcmaxthreads=7): kernelopts(gcthreadmemorysize=2^17):
'numcpus'=kernelopts(numcpus), 'gcmaxthreads'=kernelopts(gcmaxthreads),
'gcthreadmemorysize'=kernelopts(gcthreadmemorysize);

subtask := proc(lo,hi,AA,nn) local i;
  AA[1,lo..hi] := Array([seq(i^10*n, i=lo..hi)]);
  AA[2,lo..hi] := Array([seq(length(AA[1,i])-1, i=lo..hi)]);
  AA[3,lo..hi] := Array([seq(iquo(AA[1,i],10^(AA[2,i]-2)), i=lo..hi)]);
  AA[4,lo..hi] := Array([seq(irem(AA[1,i],1000), i=lo..hi)]);
  return NULL;
end proc:                  
fulltask := proc(num::posint,n::posint,N::posint) local A,i;
  A := Array(1..4,1..num,order=C_order);
  Threads:-Task:-Start(null, Task=[subtask, 1, ceil(num/N), A, n],
                       seq(Task=[subtask, ceil((i-1)*num/N)+1, ceil(i*num/N), A, n],
                           i=2..N));
  return A;
end proc:

nelems := 10000000:
n := 374894756873546859847556:
N := kernelopts(numcpus):
A,RT,GCRT := CodeTools:-Usage( Threads:-Task:-Start( fulltask, nelems, n, N),
                               output=[output, realtime, gcrealtime], quiet ):
print(sprintf("total real time: %a secs   gc real time: %a secs",
              evalf[4](RT), evalf[4](GCRT)));

numcpus = 7, gcmaxthreads = 7, gcthreadmemorysize = 131072

"total real time: 17.54 secs   gc real time: 9.306 secs"

 

Download threadthingm.mw

On these Linux machines the best performing variant above was about 3.1 times as fast as the first (unoptimzed) variant. The timings do change slightly upon re-execution, however. It was also faster than other versions posted so far for this Question. It would be interesting to see results for a 4-core (modern) hyperthreaded host on 64bit MS-Windows.

Using Maple versions 17.02 (2013) through to 2019.0,

solve([x^2+y^2+z^2-4*x+6*y-2*z-11 = 0, 2*x+2*y-z = -18],
      [x, y, z], parametric, real);

           [[x = -4/3, y = -19/3, z = 8/3]]

 

  • Widget driven interface (displayed inline in worksheet in Maple GUI or Maple Cloud or Maple Player, created via source or mouse layout):
       Embedded Components, a newer alternative to Maplets
    See Help topics:
       - EmbeddedComponents
       - examples,ProgrammaticContentGeneration
       - DocumentTools,Layout
  • "Clickable Math" access to commands applied to 2D output or input (right-click popup menu in Maple 2017 and earlier, and side-panel menu in Maple 2018 and later):
       customized Context-Menu entries
    See Help topics:
       - ContextMenu
       - examples,ContextMenu
    See more ideas:
       - package-specific context-menu additions
       - context-menu additions for "phasor" package (see ModuleLoad)
       - context driven help

Here are a few more workarounds.

restart;
solve({x^2 = 2, x^3 = sqrt(2)^3}, [x,W]):
subs((W=W)=NULL, %);

                         [[     (1/2)]]
                         [[x = 2     ]]

Of course, these next two are restricted to a particular kind of problem or domain.

restart;
SolveTools:-PolynomialSystem({x^2 = 2, x^3 = sqrt(2)^3});

                          /     (1/2)\ 
                         { x = 2      }
                          \          / 

restart;
solve({x^2 = 2, x^3 = sqrt(2)^3}, [x], real);

                         [[     (1/2)]]
                         [[x = 2     ]]

When solve behaves weirdly I usually look at eliminate as well.

Call the randomize command first with the same seed value.

restart:

kernelopts(version);

`Maple 2019.0, X86 64 LINUX, Mar 9 2019, Build ID 1384062`

(1)

with(DETools):

a1:=m^2*((1-N)/(2-N))*(r/2)*Dp:
a2:=((1-N)*m^2/(2-N))*((G[r]*(Nb-Nt)/64)*(r^5/6-h^2*r/2)-(B[r]/4)*(r^3/4-h^2*r/2)*(Nt/Nb)):
a3:=a1-a2:
a4:=r^2*a3:
DE1:=r^2*diff(v(r),r,r)+r*diff(v(r),r)-(m^2*r^2+1)*v(r)=a4:

b1:=dsolve(DE1,v(r));

v(r) = BesselI(1, m*r)*_C2+BesselK(1, m*r)*_C1-(1/128)*(-1+N)*((G[r]*(-(1/3)*r^4+h^2)*Nb^2+(-(-(1/3)*r^4+h^2)*Nt*G[r]+64*Dp)*Nb-16*B[r]*(h^2-(1/2)*r^2)*Nt)*m^4+(-8*Nb^2*r^2*G[r]+8*Nb*Nt*r^2*G[r]+64*Nt*B[r])*m^2-64*Nb*G[r]*(Nb-Nt))*r/(Nb*m^4*(-2+N))

(2)

sort(N-1,order=plex(N)):
sort(N-2,order=plex(N)):

 

F:=proc(u) local hu,ho;
     if u::`*` then
       hu,ho:=selectremove(t->depends(t,[r,m]) or t::rational,u);
       sort(expand(`*`(op(hu))),order=plex(r))*`*`(op(ho));
     else u; end if;
end proc:

 

collect(rhs(b1), [BesselI, BesselK, G[r], B[r]], u->F(simplify(u,size)));

BesselI(1, m*r)*_C2+BesselK(1, m*r)*_C1+((1/384)*r^5+(1/16)*r^3/m^2-(1/128)*h^2*r+(1/2)*r/m^4)*(N-1)*(Nb-Nt)*G[r]/(N-2)+(-(1/16)*r^3+(1/8)*h^2*r-(1/2)*r/m^2)*(N-1)*Nt*B[r]/(Nb*(N-2))-(N-1)*Dp*r/(2*N-4)

(3)

 

Download collect_example.mw


[edit] That F procedure is only there to get the polynomial terms in r in the form you showed. You may wish to omit it. Ie,

collect(rhs(b1), [BesselI, BesselK, G[r], B[r]], factor);

BesselI(1, m*r)*_C2+BesselK(1, m*r)*_C1-(1/384)*(N-1)*(-m^4*r^4+3*h^2*m^4-24*m^2*r^2-192)*r*(Nb-Nt)*G[r]/(m^4*(N-2))+(1/16)*(N-1)*Nt*(2*h^2*m^2-m^2*r^2-8)*r*B[r]/(m^2*Nb*(N-2))-(1/2)*(N-1)*Dp*r/(N-2)

(1)

collect(rhs(b1), [BesselI, BesselK, G[r], B[r]], u->simplify(u,size));

BesselI(1, m*r)*_C2+BesselK(1, m*r)*_C1-(1/128)*(N-1)*(-64+(-(1/3)*r^4+h^2)*m^4-8*m^2*r^2)*r*(Nb-Nt)*G[r]/(m^4*(N-2))+(1/16)*(N-1)*Nt*(2*h^2*m^2-m^2*r^2-8)*r*B[r]/(m^2*Nb*(N-2))-(N-1)*Dp*r/(2*N-4)

(2)

 

Is this what you want?

restart;

with(plots):

R := 5; alpha := (1/9)*Pi; beta := (1/3)*Pi; n := 100; dt := 2*Pi/n;

5

(1/9)*Pi

(1/3)*Pi

100

(1/50)*Pi

C1 := plot([R*cos(t), R*sin(t), t = 0 .. 2*Pi], color = blue):

local O := [0, 0];

M := [R*cos(beta), R*sin(beta)];

[5/2, (5/2)*3^(1/2)]

OM := plot([O, M]):

A := [R*cos(alpha), R*sin(alpha)];

[5*cos((1/9)*Pi), 5*sin((1/9)*Pi)]

B := [R*cos(alpha+Pi), R*sin(alpha+Pi)];

[-5*cos((1/9)*Pi), -5*sin((1/9)*Pi)]

AB := plot([A, B]):

P := [R*cos(t0*dt), R*sin(t0*dt)];

[5*cos((1/50)*t0*Pi), 5*sin((1/50)*t0*Pi)]

Q := [R*cos(dt*t0+Pi), R*sin(dt*t0+Pi)];

[-5*cos((1/50)*t0*Pi), -5*sin((1/50)*t0*Pi)]

tp := textplot([[A[1]+.3, A[2], "A"], [B[1]-.3, B[2], "B"], [M[1]+.3, M[2]+.3, "M"]]):

 

all := [seq(display(AB, OM, C1, tp,
                    plot([M, P]), plot([M, Q]),
                    plot([P, Q], color = green)),
            t0 = 0 .. n)]:

 

display(all, insequence = true);


OManim.mw

The behavior in this attachment may be what the author of Statistics:-DataSummary had intended.

QuestionDataWeights_acc.mw

A more robust way to discern the position of the weights optional parameter could be devised. I only checked against Maple 2019.0, Maple 2018.0 and Maple 2018.2, but chances are reasonable that weights corresponds to either PARAM(4) or PARAM(5) of the InertForm in earlier versions that have DataSummary.

I use the ToInert/subsop/FromInert to repair what seems like a typo in the source of DataSummary. On the source line of that assigns to w,ww the first argument passed to PreProcessData is edited from X to weights. See showstat(Statistics:-DataSummary) for clarification.

I'll submit a bug report.

The use of nested loops to assign distinct deep copies of a Record to the Array entries seems somewhat unnatural and stilted to me.

The Array command accepts an initializer procedure as an optional argument, which can be used to populate each entry. The distinct Records (and Vectors in their fields) can be created on the fly.

This gets rid of the need for nested loops (as many as there are dimensions) to assign all the entries their Records. It gets rid of the need for copy or deep copying of a Record. The values of the indices (i,j,k) could be passed to the constructing procedure, in the case that some Record field were desired to be set up depending on their values.

The Records could also share some Vector(s), if desired. The copy approach gets in the way of that.

restart;

create_nodeprop := () -> Record(

                    'inode'      = 0 ,
                    'nnb'        = 6 ,
                    'nind'       = Vector( 1 .. 6 ,  0  ) ,
                    'U'          = 100. ,
                    'Ul'         = 100. ,
                    'Uprev'      = 100. ,
                    'Ueq'        = 100. ,
                    'Ueql'       = 100. ,
                    'blin'       = 0.   ,
                    'q'          = 0.   ,
                    'rx'         = 2.   ,
                    'ry'         = 2.   ,
                    'rz'         = 2.   ,
                    'rfin'       = Vector( 1 .. 6 , -1. ) ,
                    'conduct'    = Vector( 1 .. 6 , -1. ) ,
                    'c'          = 2.   ,
                    'tau'        = 0.   ,
                    'k'          = 0.   ,
                    'sumconduct' = 0.   ,
                    'sharedB'    = `if`(B::rtable,ArrayTools:-Alias(B),NULL)
                  ) :

nx  := 10  :    # X direction dimension.
ny  :=  4  :    # Y direction dimension.
nz  :=  1  :    # Z direction dimension.

 

B := Vector([1,2,3,4]):

 

node := Array( 1 .. ny , 1 .. nx , 1 .. nz,
               (i,j,k) -> create_nodeprop() ) :

ic := 0 :
for k from 1 to nz do
    for j from 1 to ny do
        for i from 1 to nx do

            ic := ic + 1 :
            
            node[  j , i , k ]:-inode := ic :
            
        od : # for i from 1 to nx do
    od : # for j from 1 to ny do
od : # for k from 1 to nz do

for k from 1 to nz do
    for j from 1 to ny do
        for i from 1 to nx do

            ic := 0 :

#          +X
            if ( i + 1 <= nx ) then
                 ic := ic + 1 :
                 node[ j , i , k ]:-nind[ic] := node[ j , i + 1 , k ]:-inode :
            end if : # if ( i + 1 <= nx )

#          -X
            if ( 1 <= i - 1 ) then
                 ic := ic + 1 :
                 node[ j , i , k ]:-nind[ic] := node[ j , i - 1 , k ]:-inode :
            end if : # if ( 1 <= i - 1 )

#          +Y
            if ( j + 1 <= ny ) then
                 ic := ic + 1 :
                 node[ j , i , k ]:-nind[ic] := node[ j + 1 , i , k ]:-inode :
            end if : # if ( j + 1 <= ny )

#          -Y
            if ( 1 <= j - 1 ) then
                 ic := ic + 1 :
                 node[ j , i , k ]:-nind[ic] := node[ j - 1 , i , k ]:-inode :
            end if : # if ( 1 <= j - 1 )

#          +Z
            if ( k + 1 <= nz ) then
                 ic := ic + 1 :
                 node[ j , i , k ]:-nind[ic] := node[ j , i , k + 1 ]:-inode :
            end if : # if ( k + 1 <= nz )

#          -Z
            if ( 1 <= k - 1 ) then
                 ic := ic + 1 :
                 node[ j , i , k ]:-nind[ic] := node[ j , i , k - 1 ]:-inode :
            end if : # if ( 1 <= k - 1 )

            node[ j , i , k ]:-nnb := ic : # Store number of neighbors.

        od : # for i from 1 to nx do
     od : # for j from 1 to ny do
od : # for k from 1 to nz do

 

node[2,3,1]:-nind

Vector(6, {(1) = 14, (2) = 12, (3) = 23, (4) = 3, (5) = 0, (6) = 0})

 

node[2,3,1]:-nind   should be [ 14 , 12 , 23 , 3 , 0 , 0 ]

 

 

node[3,6,1]:-nind

Vector(6, {(1) = 27, (2) = 25, (3) = 36, (4) = 16, (5) = 0, (6) = 0})

 

node[ 3,6,1 ]:-nind   should be [ 27 , 25 , 36 ,16 , 0 , 0 ]

 

``

node[ 1 , 1 , 1 ]:-inode;

1

node[ 2 , 2 , 1 ]:-inode;

12

 

We've already seen that the Records get their own distinct Vectors.

But they can also share certain designated Vectors.

 

node[1,2,1]:-sharedB[3] := 27:

node[4,3,1]:-sharedB;

Vector(4, {(1) = 1, (2) = 2, (3) = 27, (4) = 4})

 

We can even unassign the name global name `B`, for re-use.

 

B := 'B';

B

node[4,3,1]:-sharedB[2] := 64:

node[1,2,1]:-sharedB;

Vector(4, {(1) = 1, (2) = 64, (3) = 27, (4) = 4})

 

``

Download Why_initializer.mw

Are you looking for something like this?

GAMMA(3/2);

                      1/2
                    Pi
                    -----
                      2

convert((x)!, GAMMA);

                 GAMMA(x + 1)                                                  
                                                                                                              
convert((1/2)!, GAMMA);

                      1/2
                    Pi
                    -----
                      2
First 156 157 158 159 160 161 162 Last Page 158 of 336