Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 175 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

And neither is "correct"...

tomleslie 13876
December 10 2018
1 1

Your expression 'f' has no solution for omega>0. Check the semi-logarithmic graph using

plots:-semilogplot(lhs(f), omega=1e-6..1e9);

Note that this "asymptotes" to 1 for large omega. Technically this is the "result" returned by DirectSearch:SolveEquations(), which reports a "residual" of 1, together with almost any value of omega which is "large enough". Whenever a DirectSearch command returns a large residual (the firt entry in the output list), then you should be suspicious

The expression does have a solution for small negative omega, which you can display with

plot( lhs(f), omega=-1e-05..-1e-06);

DirectSearch:-SolveEquations() will find this if you appropriately define the range, with

DirectSearch:-SolveEquations([f], {omega = -1 .. 0});

 

It is nearly always better...

tomleslie 13876
December 10 2018
2 0

to use add() rather than sum() whenever there are a finite (not too large?) number of terms. It also avoids some issues you had with premature evaluation of the summand.

The attached seems to work

> 

restart:
J:= proc(a,b,c)
        local u,v,p,w;
        p:=(a+b+c):
        if   type(p, even)
        then u:=b+c-a:
             v:=c+a-b:
             w:=a+b-c:
             return (-1)^p*sqrt(binomial(2*u,u)*binomial(2*v,v)*binomial(2*w,w)/((2*p+1)*binomial(2*p,p)))
        else return 0
        end if
    end proc:
    aa:=2:bb:=4:
    add( J(aa,bb,cc)*J(aa,bb,cc)*sqrt(2*(2*aa+1)*(2*bb+1)*(2*cc+1))*psi(1,cc),
         cc=abs(bb-aa)..(aa+bb)
        );

(2450/7293)*2^(1/2)*psi(1, 2)+(324/4199)*10^(1/2)*psi(1, 4)+(297/7429)*130^(1/2)*psi(1, 6)

 (1)
> 

 


 

Download getSum.mw

Use sum()...

tomleslie 13876
December 10 2018
1 0

Maple is reasonably good at finding closed-form expressions for sums. So for you example case

S:=simplify(sum(k^2, k=1..n));

will return the desired formula.

If you plan on obtaining values of this sum for a variety of upper indexes, then it is probably easier to write is as a function, as in

SS:=n->sum(k^2, k=1..n);

Then SS(4) will return 30 etc

 

Complicated way to do something simple??...

tomleslie 13876
December 09 2018
1 3

Defining the relevant vector and performing the "outer product" to get the matrix you want is trivial - takes one line for each.

All you then need to do is come up with a more convenient way to generate the individual entries which you have written as psi[n,m](t). I suspect the latter would more simply done if you didn't use 'indices' but considered these as arguments - psi(n, m, t). Every necessary entry could then be generated by one simple function definition like

psi:=( p, q, t )-> piecewise( t<=f(p,q),  0
                                              t<  g(p,q), h(p,q,t)
                                                               0
                                           )

See the attached for this partial solution. (The matrices look neater in a worksheet than they do on this site - honest)

If you want the "complete" solution then provide a readable definition for all functions psi[n,m](t) for any n,m

> 

  restart;
  with(LinearAlgebra):
#
# This works for any k, M so the next two values
# are just for example purposes
#
  k:=2;
  M:=3;
#
# Define the vector and generate the outer product
#
  Psi__t:=Vector[row]([seq( seq( psi[p,q](t), q=0..M-1), p=1..k)]);
  prd:=OuterProductMatrix(Psi__t, Transpose(Psi__t));

k := 2

  

M := 3

  

Vector[row](6, {(1) = psi[1, 0](t), (2) = psi[1, 1](t), (3) = psi[1, 2](t), (4) = psi[2, 0](t), (5) = psi[2, 1](t), (6) = psi[2, 2](t)})

  

Matrix(%id = 18446744074396598510)

 (1)
> 

#
# It would probably be easier to change the 'indices'
# on psi to arguments - ie rather than having
# psi[i,j](t), use psi(i,j,t) because then the actual
# functions could be (probably) be produced as a
# one-liner like
#
# psi:=( p, q, t )-> piecewise( t<=f(p,q),  0
#                               t< g(p,q), h(p,q,t)
#                               0
#                             )
#
# where f(), g() , and h() are 'known' expressions
# which depend on (p,q)
#
# Notice that with this usage, defining the relevant
# vector and outer product doesn't change (much)
#
  Psi__t:=Vector[row]([seq( seq( psi(p,q, t), q=0..M-1), p=1..k)]);
  prd:=OuterProductMatrix(Psi__t, Transpose(Psi__t));

Vector[row](6, {(1) = psi(1, 0, t), (2) = psi(1, 1, t), (3) = psi(1, 2, t), (4) = psi(2, 0, t), (5) = psi(2, 1, t), (6) = psi(2, 2, t)})

  

Matrix(%id = 18446744074381243326)

 (2)
> 

 

> 

 


Download the_answer.mw

Some random(?) plots...

tomleslie 13876
December 09 2018
1 0

It is still not bvious (to me at least) what you what to plot against what please try to be soecific.

The attatched shows a couple of possibilities

  1. V(t) versus t in the range t=0..1, for various values of U[0]
  2. V(t=1) versus U[0]

If neither of the above the above is what you require then I'm pretty sure I can produce what you want - provided you can specify it, clearly and accurately

> 

  restart;
#
# Define equations
#
eqs:= diff(r[1](t),t$2)+(.3754075596+.3045306134*U[0])*(diff(r[1](t),t))+740.0274933*r[1](t)-0.1181114312e-2*V(t) = -0.5234785518e-2*U[0]^2*q(t),
      2.030000000*10^(-8)*(diff(V(t), t))+4.065040650*10^(-7)*V(t)+0.1181114312e-2*(diff(r[1](t), t)) = 0,
      diff(q(t), t, t)-4.997988315*U[0]*(diff(q(t), t))+277.5543022*U[0]^2*q(t) = 0;
#
# Define initial conditions
#
  ics := V(0) = 0, q(0) = 0.1e-3, (D(q))(0) = 0, r[1](0) = 0, (D(r[1]))(0) = 0;
#
# Solve ode system with U[0] as parameter
#
  sol:=dsolve( {eqs, ics}, numeric, parameters=[U[0]]);

diff(diff(r[1](t), t), t)+(.3754075596+.3045306134*U[0])*(diff(r[1](t), t))+740.0274933*r[1](t)-0.1181114312e-2*V(t) = -0.5234785518e-2*U[0]^2*q(t), 0.2030000000e-7*(diff(V(t), t))+0.4065040650e-6*V(t)+0.1181114312e-2*(diff(r[1](t), t)) = 0, diff(diff(q(t), t), t)-4.997988315*U[0]*(diff(q(t), t))+277.5543022*U[0]^2*q(t) = 0

  

V(0) = 0, q(0) = 0.1e-3, (D(q))(0) = 0, r[1](0) = 0, (D(r[1]))(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 := [U[0] = `U[0]`]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 5, (2) = 5, (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) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .0, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..6, {(1) = 0., (2) = 0.1e-3, (3) = 0., (4) = 0., (5) = 0., (6) = Float(undefined)})), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..5, {(1) = .1, (2) = .1, (3) = .1, (4) = .1, (5) = .1}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0}, datatype = float[8], order = C_order), Array(1..5, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0}, datatype = float[8], order = C_order), Array(1..5, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0}, datatype = integer[8]), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..10, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = V(t), Y[2] = q(t), Y[3] = diff(q(t),t), Y[4] = r[1](t), Y[5] = diff(r[1](t),t)]`; YP[1] := -20.024830788177339901*Y[1]-58182.971034482758621*Y[5]; YP[3] := 4.997988315*Y[6]*Y[3]-277.5543022*Y[6]^2*Y[2]; YP[5] := -0.5234785518e-2*Y[6]^2*Y[2]-(.3754075596+.3045306134*Y[6])*Y[5]-740.0274933*Y[4]+0.1181114312e-2*Y[1]; YP[2] := Y[3]; YP[4] := Y[5]; 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] = q(t), Y[3] = diff(q(t),t), Y[4] = r[1](t), Y[5] = diff(r[1](t),t)]`; YP[1] := -20.024830788177339901*Y[1]-58182.971034482758621*Y[5]; YP[3] := 4.997988315*Y[6]*Y[3]-277.5543022*Y[6]^2*Y[2]; YP[5] := -0.5234785518e-2*Y[6]^2*Y[2]-(.3754075596+.3045306134*Y[6])*Y[5]-740.0274933*Y[4]+0.1181114312e-2*Y[1]; YP[2] := Y[3]; YP[4] := Y[5]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..6, {(1) = 0., (2) = 0., (3) = 0.1e-3, (4) = 0., (5) = 0., (6) = 0.}); _vmap := array( 1 .. 5, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4), ( 5 ) = (5) ] ); _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 _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, V(t), q(t), diff(q(t), t), r[1](t), diff(r[1](t), t)], (4) = [U[0] = `U[0]`]}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (1)
> 

#
# Plot solutions of V(t) versus t for the parameter
# U[0]=0..1 in steps of 0.2.
#
  cols:=[ red, orange, yellow, green, blue, black]:
  for j from 0 to 5 by 1 do
      sol(parameters=[U[0]=j/5]);
      plt[j]:= plot( 'rhs(sol(t)[2])', t=0..1, color=cols[j+1])
  od:
  plots:-display( [seq(plt[j], j=0..5)],
                  legend=[seq(typeset( "U[0]=",j/5), j=0..5)],
                  legendstyle=[font=[times, bold, 20]],
                  title="V(t) versus t",
                  titlefont=[times, bold, 20],
                  axes=boxed
                );

 
> 

#
# Plot the value of V(t=1) as a function of
# U[0]
#
  pts:=NULL:
  for j from 0 to 100 by 1 do
       sol(parameters=[U[0]=j/100]);
       pts:=pts, [j/100,rhs(sol(1)[2])]
  od:
  plots:-pointplot( [pts],
                    title="V(t=1) versus U[0]",
                    titlefont=[times, bold, 20],
                    labels=[U[0], V(1) ],
                    labelfont=[times, bold, 16],
                    style =line);;

 
> 

 

Download odePlots.mw

check...

tomleslie 13876
December 06 2018
1 1

kernelopts(datalimit)

This *ought* to return infinity, indicating (I think) that no software limit is being applied

Assuming that this command returns infinity, then it would seem that you might be hitting a hardware limit so the obvious debug stages would be

  1. Why is your code burning so much memory - can it be optimised to use less? Throwing more resources (eg cpu cores, RAM) at a problem should never be your first response!
  2. In an earlier thread you were interested in running parallel computations - are you doing so now? If you are, consider running just one of the parallel tasks, and check how much memory is being consumed. Then apply (1) above to reduce the memory requirement. In a Grid-based (ie separate mservers) parallel implementation the amount of memory required is probably (more or less) NumberOfTasks*MemoryPerTask. If the MemoryPerTask is high, then you might have to reduce the NumberOfTasks

One way...

tomleslie 13876
December 06 2018
2 0

The most awkward part is to combine/convert all the exp() terms into the simplest combination of cosh()/sinh(). There may be neater ways to do this, but the attached *works*

(It wouldn't surprise me if the necessary simplifications were slightly different in other Maple versions - the attached uses MAple 2018.2)

> 

  restart;
  kernelopts(version);
  with(inttrans):
  f:= piecewise
      ( t<0, 0,
        t<1, cosh(t),
        0
      );
  ans1:= simplify
         ( convert
           ( fouriersin
             ( convert(f, exp),
               t,
               s
             ),
             trig
           )
         );
  ans2:= simplify
         ( convert
           ( fouriercos
             ( convert(f, exp),
               t,
               s
             ),
             trig
           )
         );

`Maple 2018.2, X86 64 WINDOWS, Nov 16 2018, Build ID 1362973`

  

piecewise(t < 0, 0, t < 1, cosh(t), 0)

  

-2^(1/2)*(s*cos(s)*cosh(1)-sin(s)*sinh(1)-s)/(Pi^(1/2)*(s^2+1))

  

2^(1/2)*(s*sin(s)*cosh(1)+cos(s)*sinh(1))/(Pi^(1/2)*(s^2+1))

 (1)
> 

 

Download sincostrans.mw

Any size matrix is allowed...

tomleslie 13876
December 05 2018
1 0

By default, Maple restricts the maximum size of displayed matrices to 10*10 - just to stop unwary users from displaying really big matrices. You can change this display default by resetting the interface variable 'rtablesize'. So setting rtalbesize=20, the following will allow a 15*15 matrix to be displayed

with(LinearAlgebra):
interface(rtablesize=20):
RandomMatrix(15,15);

Note that this limitation only affects the display, not the ability to define or manipulate large matrices. After all if you were dealing with 10000*10000, matrices, I doubt if you would want them displayed -unless you have a really big monitor ;-)

I'm betting Maple is correct...

tomleslie 13876
December 05 2018
1 1

You state

In calculating the real part I see two arguments

Err no - when I chck the contents of the function being plotted with indets(.., 'name') it returns one variable, ie theta. Check the indets(zz, 'name') command in the attached.

You also state

the plot is not smooth

Now this can arise for two reasons

  1. There are too few points being plotted. Increasing the number of plotted points by adding the numpoints=10000 option to the plot command will determine if this is a plotting issue - It isn't
  2. The function you want to plot contains discontinuities. Now these can be difficult to spot, but I always get suspicious when I see terms incolving the two-argument 'arctan()' function. One of the terms in your expression is sin(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037)). In the attached, I have added a plot of this term on its own, over the range of theta which you are using. It produces a really nice discontinuity. There are several such two -argument arctan() functions in the expression your are plotting, each of whihc has (at least) the potential to introduce discontinuities, so I'm not really surrised that yor final graph has  a few "steps" in it

Check the attached


 

NULL

 

 
> 

##Toya complex variable method

> 

NULL

> 

restart;

> 

stress_c:=-(1+1/nu_c)*k*p2*zeta_c/2;

-(1/2)*(1+1/nu_c)*k*p2*zeta_c

(1.1)
> 

p2:=(c0_c-d_1c/k)*(z-a*(cos(alpha)+2*lambda*sin(alpha)))+(1-k)/k*a*(N_infty-T_infty)*exp(2*I*phi_c+2*lambda*(alpha-Pi))*((a*(cos(alpha)-2*lambda*sin(alpha)))/z-a^2/z^2)

(c0_c-d_1c/k)*(z-a*(cos(alpha)+2*lambda*sin(alpha)))+(1-k)*a*(N_infty-T_infty)*exp((2*I)*phi_c+2*lambda*(alpha-Pi))*(a*(cos(alpha)-2*lambda*sin(alpha))/z-a^2/z^2)/k

 (1.2)
> 

NULL

> 

z := exp(I*theta)

exp(I*theta)

(1.3)
> 

NULL

> 

k := beta_c/(1+nu_c)

beta_c/(1+nu_c)

(1.4)
> 

nu_c := (kappa2*mu+mu2)/(kappa*mu2+mu)

(kappa2*mu+mu2)/(kappa*mu2+mu)

(1.5)
> 

d_1c := (N_infty+T_infty)*(1/2)

(1/2)*N_infty+(1/2)*T_infty

(1.6)
> 

lambda := -evalf(ln(nu_c)/(2*Pi))

-.1591549430*ln((kappa2*mu+mu2)/(kappa*mu2+mu))

(1.7)
> 

NULL

> 

beta_c := mu*(1+kappa2)/(kappa*mu2+mu)

mu*(1+kappa2)/(kappa*mu2+mu)

(1.8)
> 

zeta_c := ((z-a*exp(I*alpha))/(z-a*exp(-I*alpha)))^(I*lambda)/((z-a*exp(I*alpha))^.5*(z-a*exp(-I*alpha))^.5)

((exp(I*theta)-a*exp(I*alpha))/(exp(I*theta)-a*exp(-I*alpha)))^(-(.1591549430*I)*ln((kappa2*mu+mu2)/(kappa*mu2+mu)))/((exp(I*theta)-a*exp(I*alpha))^.5*(exp(I*theta)-a*exp(-I*alpha))^.5)

(1.9)
> 

NULL

> 

c0_c := G_c+I*H_c

G_c+I*H_c

(1.10)
> 

G_c:=(0.5*(T_infty+N_infty)*(1-(cos(alpha)+2*lambda*sin(alpha))*exp(2*lambda*(evalf(Pi)-alpha)))-0.5*(1-k)*(1+4*lambda^2)*(N_infty-T_infty)*(sin(alpha))^2*cos(2*phi_c))/(2-k-k*(cos(alpha)+2*lambda*sin(alpha))*exp(evalf(2*lambda*(Pi-alpha))));

(.5*(N_infty+T_infty)*(1-(cos(alpha)-.3183098860*ln((kappa2*mu+mu2)/(kappa*mu2+mu))*sin(alpha))*exp(-.3183098860*ln((kappa2*mu+mu2)/(kappa*mu2+mu))*(3.141592654-alpha)))-.5*(1-mu*(1+kappa2)/((kappa*mu2+mu)*(1+(kappa2*mu+mu2)/(kappa*mu2+mu))))*(.1013211835*ln((kappa2*mu+mu2)/(kappa*mu2+mu))^2+1)*(N_infty-T_infty)*sin(alpha)^2*cos(2*phi_c))/(2-mu*(1+kappa2)/((kappa*mu2+mu)*(1+(kappa2*mu+mu2)/(kappa*mu2+mu)))-mu*(1+kappa2)*(cos(alpha)-.3183098860*ln((kappa2*mu+mu2)/(kappa*mu2+mu))*sin(alpha))*exp(-.3183098860*ln((kappa2*mu+mu2)/(kappa*mu2+mu))*(3.141592654-1.*alpha))/((kappa*mu2+mu)*(1+(kappa2*mu+mu2)/(kappa*mu2+mu))))

 (1.11)
> 

H_c:=0.5*(1-k)*(1+4*lambda^2)*(-T_infty+N_infty)*(sin(alpha))^2*sin(2*phi_c)/(k*(1+(cos(alpha)+2*lambda*sin(alpha))*exp(2*lambda*(evalf(Pi)-alpha))));

.5*(1-mu*(1+kappa2)/((kappa*mu2+mu)*(1+(kappa2*mu+mu2)/(kappa*mu2+mu))))*(.1013211835*ln((kappa2*mu+mu2)/(kappa*mu2+mu))^2+1)*(N_infty-T_infty)*sin(alpha)^2*sin(2*phi_c)*(kappa*mu2+mu)*(1+(kappa2*mu+mu2)/(kappa*mu2+mu))/(mu*(1+kappa2)*(1+(cos(alpha)-.3183098860*ln((kappa2*mu+mu2)/(kappa*mu2+mu))*sin(alpha))*exp(-.3183098860*ln((kappa2*mu+mu2)/(kappa*mu2+mu))*(3.141592654-alpha))))

 (1.12)
> 

##Input

> 

alpha:=evalf(Pi/6)

.5235987758

(1.13)
> 

phi_c:=alpha;

.5235987758

(1.14)
> 

N_infty:=0;

0

(1.15)
> 

T_infty:=1;

1

(1.16)
> 

a:=1;nu2:=22/100;kappa2:=3-4*nu2;nu:=35/100;kappa:=3-4*nu;mu:=239/100;mu2:=442/10;

1

 

11/50

 

53/25

 

7/20

 

8/5

 

239/100

 

221/5

(1.17)
> 

NULL

> 

stress_c

-(9321/123167)*(((.5586916801-.5*(.8660254037-.1591549431*ln(123167/182775))*exp(-.8333333329*ln(123167/182775))+0.5946710490e-2*ln(123167/182775)^2)/(22817/11767-(717/11767)*(.8660254037-.1591549431*ln(123167/182775))*exp(-.8333333329*ln(123167/182775)))-(1.668336947*I)*(.1013211835*ln(123167/182775)^2+1)/(1+(.8660254037-.1591549431*ln(123167/182775))*exp(-.8333333329*ln(123167/182775)))-11767/1434)*(exp(I*theta)-.8660254037+.1591549431*ln(123167/182775))-(11050/717)*exp(1.047197552*I+.8333333328*ln(123167/182775))*((.8660254037+.1591549431*ln(123167/182775))/exp(I*theta)-1/(exp(I*theta))^2))*((exp(I*theta)+(-.8660254037-.5000000002*I))/(exp(I*theta)+(-.8660254037+.5000000002*I)))^(-(.1591549430*I)*ln(123167/182775))/((exp(I*theta)+(-.8660254037-.5000000002*I))^.5*(exp(I*theta)+(-.8660254037+.5000000002*I))^.5)

 (1.18)
> 

assume((1/6)*Pi < theta, theta < 2*Pi-(1/6)*Pi)

> 

zz := simplify(evalc(Re(stress_c))); indets(zz, 'name')

-0.1184144190e-10*((((-0.5262862456e-1*cos(theta)^7-0.3541068761e-1*cos(theta)^6+.4653205519*cos(theta)^5+(84.91408010*sin(theta)-.1768534236)*cos(theta)^4-.9987575018*cos(theta)^3+(-84.78230062*sin(theta)+0.1188037308e12)*cos(theta)^2+(-0.6859136587e11*sin(theta)-0.5228707886e11)*cos(theta)-0.5515218277e11-0.2348357835e11*sin(theta))*cos(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037))+(-9.444917506*cos(theta)^7+0.3941955295e-1*cos(theta)^6+84.08911190*cos(theta)^5+(-81.24793278*sin(theta)+88.26451743)*cos(theta)^4-84.11584179*cos(theta)^3+(.1320908588*sin(theta)-0.6859136597e11)*cos(theta)^2+(-0.1188037310e12*sin(theta)+0.7857647929e11)*cos(theta)-0.1310336552e11+0.4313662459e11*sin(theta))*sin(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037)))*sin(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037))+(1.000000000*cos(theta)^7-0.3941955295e-1*cos(theta)^6+.3730125636*cos(theta)^5+(67.69143114*sin(theta)-88.26451777)*cos(theta)^4+0.5657578745e-1*cos(theta)^3+(16.75774422*sin(theta)+0.6859136597e11)*cos(theta)^2+(0.1188037311e12*sin(theta)-0.7857647933e11)*cos(theta)+0.1310336557e11-0.4313662455e11*sin(theta))*cos(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037))*cos(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037))+(-0.5262862456e-1*cos(theta)^7-0.3541067072e-1*cos(theta)^6+.7986537095*cos(theta)^5+(84.91408001*sin(theta)-.1768532293)*cos(theta)^4-.1320910201*cos(theta)^3+(-84.78230062*sin(theta)+0.1188037308e12)*cos(theta)^2+(-0.6859136586e11*sin(theta)-0.5228707882e11)*cos(theta)-0.5515218310e11-0.2348357852e11*sin(theta))*cos(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037))*sin(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037)))*cos(0.314104002e-1*ln(-1292820323.*cos(theta)+746410161.*sin(theta)+1492820323.)-0.314104002e-1*ln(-1292820322.*cos(theta)-746410161.4*sin(theta)+1492820322.))+(((-.9999999974*cos(theta)^7+33.41981464*cos(theta)^5+(-64.35809780*sin(theta)+88.26451768)*cos(theta)^4-25.00141918*cos(theta)^3+(-16.75774414*sin(theta)-0.6859136597e11)*cos(theta)^2+(-0.1188037310e12*sin(theta)+0.7857647933e11)*cos(theta)-0.1310336578e11+0.4313662444e11*sin(theta))*cos(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037))+(0.5262859922e-1*cos(theta)^7+0.2077356812e-2*cos(theta)^6-.7986538447*cos(theta)^5+(-84.50969801*sin(theta)+.1965889728)*cos(theta)^4+.1320909323*cos(theta)^3+(84.50685874*sin(theta)-0.1188037310e12)*cos(theta)^2+(0.6859136587e11*sin(theta)+0.5228707882e11)*cos(theta)+0.5515218258e11+0.2348357843e11*sin(theta))*sin(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037)))*sin(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037))+(-0.3541068677e-1*cos(theta)^6+.7985498379*cos(theta)^5+(-.1768532969+42.48739162*sin(theta))*cos(theta)^4-.1320910201*cos(theta)^3+(0.1188037310e12-42.61560461*sin(theta))*cos(theta)^2+(-0.5228707881e11-0.6859136587e11*sin(theta))*cos(theta)-0.5515218290e11-0.2348357825e11*sin(theta))*cos(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037))*cos(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037))+(33.40655361*cos(theta)^5+(-81.24793279*sin(theta)+87.91408004)*cos(theta)^4-25.33343679*cos(theta)^3+(.1320911054*sin(theta)-0.6859136596e11)*cos(theta)^2+(-0.1188037310e12*sin(theta)+0.7857647929e11)*cos(theta)-0.1310336533e11+0.4313662486e11*sin(theta))*cos(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037))*sin(.5*arctan(sin(theta)-.5000000002, cos(theta)-.8660254037)))*sin(0.314104002e-1*ln(-1292820323.*cos(theta)+746410161.*sin(theta)+1492820323.)-0.314104002e-1*ln(-1292820322.*cos(theta)-746410161.4*sin(theta)+1492820322.)))/((-sin(theta)+2.-1.732050807*cos(theta))^(1/4)*(sin(theta)+2.-1.732050807*cos(theta))^(1/4))

  

{theta}

(1.19)
> 

plot(zz, theta = (1/6)*Pi .. 2*Pi-(1/6)*Pi, numpoints = 10000)

 
> 

#
# Plot of one term extracted from OP's expression
# Note how well its discontinuity aligns with the
# "discontinuity" in the above
#
# OP's expression contains a few other two-element
# arctan functions, each of whihc has at least the
# capability to introduce similar discontinuities
#
  plot( sin(.5*arctan(sin(theta)+.5000000002, cos(theta)-.8660254037)),
        theta=Pi/6..2*Pi-Pi/6
      );

 
> 

 


 

Download arctanIssue.mw

After fixing the parenthesis problem...

tomleslie 13876
December 04 2018
1 3

mentioned by PrebenAlshom, I managed to get a solution for this system using the DirectSearch() add-on package, whihc is available (free) from the Maple Application Centre at

https://www.maplesoft.com/applications/Category.aspx?cid=1600

The residuals aren't great, but only you will know if they are acceptable

The fsolve() command might have worked - eventually, but I ran out of patience. Even the DirectSearch:-SolveEquations took about 5minutes on my (admittedly again) hardware.

See the attached, but note that you will not be able to re-execute this worksheet unless you have the DirectSearch() package installed

> 

Loading LinearAlgebra

> 

restart; lambda[0] := 0.8e-3

0.8e-3

 (1)
> 

lambda[1] := 0.8e-3

0.8e-3

 (2)
> 

lambda[2] := 0.8e-3

0.8e-3

 (3)
> 

lambda[3] := 0.8e-3; n := 1

0.8e-3

  

1

 (4)
> 

e25 := p[s] = (n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3]))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n

p[s] = (t[0]/(1-t[0])+t[1]/(1-t[1])+t[2]/(1-t[2])+t[3]/(1-t[3]))*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])

 (5)
> 

e26 := p[b] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n

p[b] = 1-(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])

 (6)
> 

e21 := r[0] = lambda[0]*(20*p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+20*r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5)+p[0, b]*(11258.065*p[s]/p[b]+469.725)*(p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5))+(469.725*(-p[0, p]*r[0]+p[0, p]+r[0]))*(-4*p[0, c]^5+p[0, c]^4+p[0, c]^3+p[0, c]^2+p[0, c])+11727.79-11727.79*p[0, c]^5*(-p[0, p]*r[0]+p[0, p]+r[0]))

r[0] = 0.160e-1*p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+0.160e-1*r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5)+0.8e-3*p[0, b]*(11258.065*p[s]/p[b]+469.725)*(p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5))+.3757800*(-p[0, p]*r[0]+p[0, p]+r[0])*(-4*p[0, c]^5+p[0, c]^4+p[0, c]^3+p[0, c]^2+p[0, c])+9.382232-9.382232*p[0, c]^5*(-p[0, p]*r[0]+p[0, p]+r[0])

 (7)
> 

e22 := r[1] = lambda[1]*(20*p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+20*r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5)+p[1, b]*(11258.065*p[s]/p[b]+469.725)*(p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5))+(469.725*(-p[1, p]*r[1]+p[1, p]+r[1]))*(-4*p[1, c]^5+p[1, c]^4+p[1, c]^3+p[1, c]^2+p[1, c])+11727.79-11727.79*p[1, c]^5*(-p[1, p]*r[1]+p[1, p]+r[1]))

r[1] = 0.160e-1*p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+0.160e-1*r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5)+0.8e-3*p[1, b]*(11258.065*p[s]/p[b]+469.725)*(p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5))+.3757800*(-p[1, p]*r[1]+p[1, p]+r[1])*(-4*p[1, c]^5+p[1, c]^4+p[1, c]^3+p[1, c]^2+p[1, c])+9.382232-9.382232*p[1, c]^5*(-p[1, p]*r[1]+p[1, p]+r[1])

 (8)
> 

e23 := r[1] = lambda[2]*(20*p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+20*r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5)+p[2, b]*(11258.065*p[s]/p[b]+469.725)*(p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5))+(469.725*(-p[2, p]*r[2]+p[2, p]+r[2]))*(-4*p[2, c]^5+p[2, c]^4+p[2, c]^3+p[2, c]^2+p[2, c])+11727.79-11727.79*p[2, c]^5*(-p[2, p]*r[2]+p[2, p]+r[2]))

r[1] = 0.160e-1*p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+0.160e-1*r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5)+0.8e-3*p[2, b]*(11258.065*p[s]/p[b]+469.725)*(p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5))+.3757800*(-p[2, p]*r[2]+p[2, p]+r[2])*(-4*p[2, c]^5+p[2, c]^4+p[2, c]^3+p[2, c]^2+p[2, c])+9.382232-9.382232*p[2, c]^5*(-p[2, p]*r[2]+p[2, p]+r[2])

 (9)
> 

e24 := r[3] = lambda[3]*(20*p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+20*r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5)+p[3, b]*(11258.065*p[s]/p[b]+469.725)*(p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5))+(469.725*(-p[3, p]*r[3]+p[3, p]+r[3]))*(-4*p[3, c]^5+p[3, c]^4+p[3, c]^3+p[3, c]^2+p[3, c])+11727.79-11727.79*p[3, c]^5*(-p[3, p]*r[3]+p[3, p]+r[3]))

r[3] = 0.160e-1*p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+0.160e-1*r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5)+0.8e-3*p[3, b]*(11258.065*p[s]/p[b]+469.725)*(p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5))+.3757800*(-p[3, p]*r[3]+p[3, p]+r[3])*(-4*p[3, c]^5+p[3, c]^4+p[3, c]^3+p[3, c]^2+p[3, c])+9.382232-9.382232*p[3, c]^5*(-p[3, p]*r[3]+p[3, p]+r[3])

 (10)
> 

e1 := t[0] = (-p[0, c]^5+1)/((1-p[0, c])*(r[0]+(1-r[0])/p[0, p]+(-p[0, c]^5+1)/(1-p[0, c])+(7*(p[0, b]*(1-r[0])+r[0]))/(2-2*p[0, b])+(15*p[0, c]^4+15*p[0, c]^3+15*p[0, c]^2+15*p[0, c])/(2-2*p[0, b])))

t[0] = (-p[0, c]^5+1)/((1-p[0, c])*(r[0]+(1-r[0])/p[0, p]+(-p[0, c]^5+1)/(1-p[0, c])+7*(p[0, b]*(1-r[0])+r[0])/(2-2*p[0, b])+(15*p[0, c]^4+15*p[0, c]^3+15*p[0, c]^2+15*p[0, c])/(2-2*p[0, b])))

 (11)
> 

e2 := t[1] = (-p[1, c]^5+1)/((1-p[1, c])*(r[1]+(1-r[1])/p[1, p]+(-p[1, c]^5+1)/(1-p[1, c])+(15*(p[1, b]*(1-r[1])+r[1]))/(2-2*p[1, b])+(31*p[1, c]^4+31*p[1, c]^3+31*p[1, c]^2+31*p[1, c])/(2-2*p[1, b])))

t[1] = (-p[1, c]^5+1)/((1-p[1, c])*(r[1]+(1-r[1])/p[1, p]+(-p[1, c]^5+1)/(1-p[1, c])+15*(p[1, b]*(1-r[1])+r[1])/(2-2*p[1, b])+(31*p[1, c]^4+31*p[1, c]^3+31*p[1, c]^2+31*p[1, c])/(2-2*p[1, b])))

 (12)
> 

e3 := t[2] = (-p[2, c]^5+1)/((1-p[2, c])*(r[2]+(1-r[2])/p[2, p]+(-p[2, c]^5+1)/(1-p[2, c])+(31*(p[2, b]*(1-r[2])+r[2]))/(2-2*p[2, b])+(511*p[2, c]^4+255*p[2, c]^3+127*p[2, c]^2+63*p[2, c])/(2-2*p[2, b])))

t[2] = (-p[2, c]^5+1)/((1-p[2, c])*(r[2]+(1-r[2])/p[2, p]+(-p[2, c]^5+1)/(1-p[2, c])+31*(p[2, b]*(1-r[2])+r[2])/(2-2*p[2, b])+(511*p[2, c]^4+255*p[2, c]^3+127*p[2, c]^2+63*p[2, c])/(2-2*p[2, b])))

 (13)
> 

e4 := t[3] = (-p[3, c]^5+1)/((1-p[3, c])*(r[3]+(1-r[3])/p[3, p]+(-p[3, c]^5+1)/(1-p[3, c])+(31*(p[3, b]*(1-r[3])+r[3]))/(2-2*p[3, b])+(511*p[3, c]^4+255*p[3, c]^3+127*p[3, c]^2+63*p[3, c])/(2-2*p[3, b])))

t[3] = (-p[3, c]^5+1)/((1-p[3, c])*(r[3]+(1-r[3])/p[3, p]+(-p[3, c]^5+1)/(1-p[3, c])+31*(p[3, b]*(1-r[3])+r[3])/(2-2*p[3, b])+(511*p[3, c]^4+255*p[3, c]^3+127*p[3, c]^2+63*p[3, c])/(2-2*p[3, b])))

 (14)
> 

e5 := p[0, p] = 1-exp(-lambda[0]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725))

p[0, p] = 1-exp(.3597800*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-9.0064520*(t[0]/(1-t[0])+t[1]/(1-t[1])+t[2]/(1-t[2])+t[3]/(1-t[3]))*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-.3757800)

 (15)
> 

e6 := p[1, p] = 1-exp(-lambda[1]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725))

p[1, p] = 1-exp(.3597800*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-9.0064520*(t[0]/(1-t[0])+t[1]/(1-t[1])+t[2]/(1-t[2])+t[3]/(1-t[3]))*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-.3757800)

 (16)
> 

e7 := p[2, p] = 1-exp(-lambda[2]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725))

p[2, p] = 1-exp(.3597800*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-9.0064520*(t[0]/(1-t[0])+t[1]/(1-t[1])+t[2]/(1-t[2])+t[3]/(1-t[3]))*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-.3757800)

 (17)
> 

e8 := p[3, p] = 1-exp(-lambda[3]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725))

p[3, p] = 1-exp(.3597800*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-9.0064520*(t[0]/(1-t[0])+t[1]/(1-t[1])+t[2]/(1-t[2])+t[3]/(1-t[3]))*(1-t[0])*(1-t[1])*(1-t[2])*(1-t[3])-.3757800)

 (18)
> 

e9 := p[0, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0])

p[0, c] = 1-(1-t[1])*(1-t[2])*(1-t[3])

 (19)
> 

e10 := p[1, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1])

p[1, c] = 1-(1-t[0])*(1-t[2])*(1-t[3])

 (20)
> 

e11 := p[2, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2])

p[2, c] = 1-(1-t[0])*(1-t[1])*(1-t[3])

 (21)
> 

e12 := p[3, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3])

p[3, c] = 1-(1-t[0])*(1-t[1])*(1-t[2])

 (22)
> 

e13 := p[0, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0]))^1

p[0, b] = 1-(1-t[1])*(1-t[2])*(1-t[3])

 (23)
> 

e14 := p[1, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1]))^1

p[1, b] = 1-(1-t[0])*(1-t[2])*(1-t[3])

 (24)
> 

e15 := p[2, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2]))^2

p[2, b] = 1-(1-t[0])^2*(1-t[1])^2*(1-t[3])^2

 (25)
> 

e16 := p[3, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3]))^6

p[3, b] = 1-(1-t[0])^6*(1-t[1])^6*(1-t[2])^6

 (26)
> 

e17 := p[0, s] = n*t[0]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0])

p[0, s] = t[0]*(1-t[1])*(1-t[2])*(1-t[3])

 (27)
> 

e18 := p[1, s] = n*t[1]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1])

p[1, s] = t[1]*(1-t[0])*(1-t[2])*(1-t[3])

 (28)
> 

e19 := p[2, s] = n*t[2]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2])

p[2, s] = t[2]*(1-t[0])*(1-t[1])*(1-t[3])

 (29)
> 

e20 := p[3, s] = n*t[3]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3])

p[3, s] = t[3]*(1-t[0])*(1-t[1])*(1-t[2])

 (30)
> 

NULL
> 

ans := DirectSearch:-SolveEquations([e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21, e22, e23, e24, e25, e26], evaluationlimit = 1000000)

ans := [0.320885135297404e-1, Vector(26, {(1) = HFloat(-0.08831144493706443), (2) = HFloat(-0.0693418665433532), (3) = HFloat(-0.0021464255952535795), (4) = HFloat(0.02525739122791197), (5) = HFloat(0.0031580748336259568), (6) = HFloat(-4.6403354143986775e-4), (7) = HFloat(0.006983895861944145), (8) = HFloat(0.003132676742497509), (9) = HFloat(-0.022270327700032677), (10) = HFloat(-0.013231467797989893), (11) = HFloat(0.04808558148212372), (12) = HFloat(0.005073701553590659), (13) = HFloat(0.018261555070266167), (14) = HFloat(0.04625740108087373), (15) = HFloat(-0.06042360717583484), (16) = HFloat(-0.00758733871746442), (17) = HFloat(0.022151692326873005), (18) = HFloat(-0.01805271745523712), (19) = HFloat(0.004837188626218878), (20) = HFloat(0.053288827972262566), (21) = HFloat(0.0014045469167574128), (22) = HFloat(-0.0016129299471571634), (23) = HFloat(0.0013595219275934608), (24) = HFloat(-0.001089044420167725), (25) = HFloat(-0.05420818842959208), (26) = HFloat(-0.05434463383178145)}), [p[b] = 0.620386094631374e-4, p[s] = -0.392427958021412e-4, p[0, b] = 0.756730548801514e-1, p[0, c] = 0.351411721098526e-1, p[0, p] = .410681795321087, p[0, s] = 0.191468649582325e-1, p[1, b] = .104150213695914, p[1, c] = 0.446613448170502e-1, p[1, p] = .407059686946021, p[1, s] = -0.215388576290325e-1, p[2, b] = 0.211148176770432e-1, p[2, c] = 0.897215717514149e-1, p[2, p] = .414507616349405, p[2, s] = 0.176078707981724e-1, p[3, b] = 0.308856640998852e-1, p[3, c] = 0.115911429905629e-1, p[3, p] = .410656397229959, p[3, s] = .101178058976535, r[0] = 3.93313839960829, r[1] = 1.78084258789638, r[2] = 4.51730734679446, r[3] = 2.82725885070482, t[0] = -0.318784641233630e-2, t[1] = -0.370036469360981e-2, t[2] = 0.133255026715183e-1, t[3] = 0.482033938004299e-1], 496524]

 (31)
> 

 


 

Download eqSolve.mw

The volume and surface area...

tomleslie 13876
December 04 2018
0 0

of a spherical cap can be computed as shown in the attached (whicjh agrees with result in your quoted reference https://en.m.wikipedia.org/wiki/Spherical_cap

> 

restart:
with(Student[Calculus1]):
assume(r>0, h>0, h<r):
SurfArea:= SurfaceOfRevolution(sqrt(r^2-(x-r)^2), x=0..h);
Volume:= VolumeOfRevolution(sqrt(r^2-(x-r)^2), x=0..h);

2*h*Pi*r

  

Pi*r*h^2-(1/3)*Pi*h^3

 (1)
> 

 

Download sphericalCap.mw

Be careful...

tomleslie 13876
December 03 2018
0 4

Threads() - doesn't always seem to provide what you might expect. I started playing around with this stuff after a post (now resolved) by vv

https://www.mapleprimes.com/questions/226018-Random-Results-Using-Threads

and an earlier post here

https://www.mapleprimes.com/questions/226016-Is-There-A-Way-To-Coerce-The-Maple-To

by the OP of this thread

I kept getting somewhat puzzling and ambiguous answers on whether using "Threads" would actually provide an improvement or not.

The attached worksheet shows five different ways to add the integers between 1 and 10^9.

The first is a "braindead" add(), just for reference. The next two are slight variations on vv's code above using Thread:-Seq(). The final two are actually taken from various Maple help/example pages

Compared with a simple add()

two of the "parallelized" methods give significant(~2X) "real time" improvements

one is about the same

one is about 6X worse

I'm not sure what other comments to make, althouugh there are a few other observations in the attached


 

> 

  restart;
> 

#
# NB I have four cores and two/threads per core
# I'm nearly certain that whilst Maple whilst may
# report this as 8 cpus, Maple will only allow
# me to have 4 simultaneos "Maple Threads", although
# when one of these "Threads" hits a processor core
# Intel's hyperthreading may allow two "processor
# threads" to work on one "Maple Thread" - or I might
# be completely wrong!!
#
  kernelopts(version);
  kernelopts(multithreaded);
  numCP:=kernelopts(numcpus);
  N:=10^9:
  f:= proc(a,b) local j;
           add(j, j=a..b)
      end proc

`Maple 2018.2, X86 64 WINDOWS, Nov 16 2018, Build ID 1362973`

  

true

  

8

  

proc (a, b) local j; add(j, j = a .. b) end proc

 (1)
> 

#
# try evaluating the same thing three different ways
#
# 1. without 'Threads'
# 2. braindead 'Threads' implementation
# 3. 'Threads' implementation with 'tasksize' control
#
# Note that (1) and (2) above take identical 'real time'
# (NB real time figures are verified by the entirely
# mechanical, "certified chronometer" on my left wrist)
#
# Limiting the task size does produces a 2X improvement,
# but the following shows that simple use of Threads:-Seq()
# is sometimes of no benefit at all in real time
#
  CodeTools:-Usage(add([seq( f(k*N/4+1,(k+1)*N/4), k=0..3)]));
  CodeTools:-Usage(add([Threads:-Seq(f(k*N/4+1,(k+1)*N/4), k=0..3)]));
  CodeTools:-Usage(add([Threads:-Seq[tasksize=1](f(k*N/4+1,(k+1)*N/4), k=0..3)]));

memory used=2.50GiB, alloc change=-3.01MiB, cpu time=53.79s, real time=53.85s, gc time=2.65s

  

500000000500000000

  

memory used=2.50GiB, alloc change=15.31MiB, cpu time=81.87s, real time=55.49s, gc time=28.27s

  

500000000500000000

  

memory used=2.50GiB, alloc change=0 bytes, cpu time=73.20s, real time=31.31s, gc time=6.11h

  

500000000500000000

 (2)
> 

#
# The above "threaded" implementations use a small(?)
# number of threads with "lengthy" adds. There is an
# example on the help page "Threads:-Task:-Start"
# which does essentially the same calculation, except
# that the add() operations are restricted to 1000
# entries, and it adds the four outputs.
#
# Hence if I understand this correctly there is a
# consequent increase in the number  of tasks
# being created to (N/1000), ie in this example to
# 10^9/1000 =10^6 tasks - gulp, can this be true!???
#
# I would have thought the task creation/management
# overhead would make this pretty slow, but in fact
# it gets quite close to the best time above??
#
  cont := proc( a, b )
                return a + b;
          end proc:
  task := proc( i, j )
                local k;
                 if   ( j-i < 1000 )
                 then return add( k, k=i..j );
                 else k := floor( (j-i)/2 )+i;
                      Threads:-Task:-Continue
                      ( cont,
                        Task=[ task, i, k ],
                        Task=[ task, k+1, j ]
                      );
                 end if;
          end proc:
  CodeTools:-Usage(Threads:-Task:-Start( task, 1, N ));

memory used=7.12GiB, alloc change=240.00MiB, cpu time=3.91m, real time=34.91s, gc time=7.04s

  

500000000500000000

 (3)
> 

#
# The last implementation which does more or less
# the same thing is taken from the final section
# of the Threads,Examples worksheet. This time
# rather than restricting the length of the add()
# function it permits the user to specify the number
# of threads, and the 'length' of the add() function
# is adjusted appropriately - so for test purposes
# I set the number of "threads" as 4
#
# Note that - in real time - this example takes 6X!!!
# as long as the *worst* of the above calculations.
# Couldn't obvipusly describe this code as a
# motivating example!!!
#
  MAdd := proc( s, e, threadnum )
               local n, i, j, out, blocks, section;
               n := e-s+1;
               out := table();
               section := floor( n/threadnum );
               blocks := [ s,
                           seq
                           ( section*i+s,
                             i=1..threadnum-1
                           ),
                           e+1
                         ];
               Threads:-Wait
               ( seq
                 ( Threads:-Create
                   ( add( j[i],
                          j[i]=blocks[i]..blocks[i+1]-1
                        ),
                     out[i]
                   ),
                   i = 1..threadnum
                 )
               );
              add( out[i], i=1..threadnum );
           end proc:
  CodeTools:-Usage(MAdd( 1, N, 4 ));

memory used=2.50GiB, alloc change=8.75MiB, cpu time=17.67m, real time=4.66m, gc time=25.38s

  

500000000500000000

 (4)
> 

 


 

Download parallel.mw

 

 

If I was looking for multiple optima...

tomleslie 13876
December 02 2018
0 0

then I would probably use the DirectSearch() package which is available (free) from the Maple Applications centre at

https://www.maplesoft.com/applications/view.aspx?SID=101333

The GlobalSearch() command from this package will find "all" local/global optima of a given function. You can specify the total number of solutions to return: the command options can be set up for integer-only options

So far as I can tell...

tomleslie 13876
December 02 2018
2 3

the Intel i7-6900K has 8 cores and 16 threads.

By default, Maple will perform a calculation using a single core - so about 1/8 of your total processor capability

If your calculation can be appropriately subdivided then you can use the Threads() package (see also the Task Programming Model)  to distribute the calculation across multiple cores which should give a significant speed-up

A numerical solution...

tomleslie 13876
December 01 2018
0 1

is simple to obtain. I just guessed(!) some boundary conditions - namely that the pendulum is stationary( ie D(theta)(0)=0 ) and held at the horizontal (ie theta=Pi/2). If you don't like these ICs, then insert your own in the attached

BTW the solution which you have obtained is not "loads of bullshit". When Maple cannot find a completely general analytic solution, then it returns a "reduced order" ode which may (or may not) be solvable with further processing or assumptions. Since you appear to be happy with a numeric solution, you don't have to worry about this.

For some reason this worksheet won't display inline, so you have to download - sorry

Download pendulum.mw

 

First 117 118 119 120 121 122 123 Last Page 119 of 207