Carl Love

Carl Love

28025 Reputation

25 Badges

12 years, 311 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

After 3 years, here is a solution using LPSolve for 10 deliveries with more than 1 truckload. (VV's combinatorial solution is for 8 deliveries. Brian's ILP solution is only of the 1-load TSP (i.e., travelling salesman problem = 1 load) subproblem.) The truck is required to return to the warehouse when it's finished each round of deliveries. The worksheet below specifically shows the solution for 2 loads (truck capacity K = 100). To increase the number of loads, set the value of lower. 
 

restart
:

#Vector of required deliveries. Index 1 is warehouse, so that entry is 0:
d:= <0, 67, 15, 2, 6, 12, 17, 2, 5, 13, 28>
:
#Number of destinations, including warehouse:
n:= numelems(d):
interface(rtablesize= n):
N:= {$1..n}: #index set, for convenience
C:= {$2..n}: #customer indices

#Truck capacity:
K:= 100
:

#min Truckloads required:
Vs:= ceil(add(d)/K);
V:= {$1..Vs}:

2

#Symmetric Matrix of costs (or distances) to drive between given locations:
Cost:= Matrix(
    n$2, shape= symmetric, scan= triangular,
    [0, op]~([
        [376, 246, 297, 599, 318, 132, 633, 478, 346, 585],
             [326, 240, 638, 457, 184, 746, 278, 314, 321],
                  [ 13, 619, 346,  45, 621, 323, 212, 376],
                       [645, 317, 128, 648, 379, 297, 337],
                            [545, 653, 482, 869, 742, 883],
                                 [440, 573, 340, 459, 638],
                                      [501, 365, 318, 406],
                                          [1045, 850, 1087],
                                                [243,  24],
                                                      [239]
    ])
);

Matrix(11, 11, {(1, 1) = 0, (1, 2) = 376, (1, 3) = 246, (1, 4) = 297, (1, 5) = 599, (1, 6) = 318, (1, 7) = 132, (1, 8) = 633, (1, 9) = 478, (1, 10) = 346, (1, 11) = 585, (2, 1) = 376, (2, 2) = 0, (2, 3) = 326, (2, 4) = 240, (2, 5) = 638, (2, 6) = 457, (2, 7) = 184, (2, 8) = 746, (2, 9) = 278, (2, 10) = 314, (2, 11) = 321, (3, 1) = 246, (3, 2) = 326, (3, 3) = 0, (3, 4) = 13, (3, 5) = 619, (3, 6) = 346, (3, 7) = 45, (3, 8) = 621, (3, 9) = 323, (3, 10) = 212, (3, 11) = 376, (4, 1) = 297, (4, 2) = 240, (4, 3) = 13, (4, 4) = 0, (4, 5) = 645, (4, 6) = 317, (4, 7) = 128, (4, 8) = 648, (4, 9) = 379, (4, 10) = 297, (4, 11) = 337, (5, 1) = 599, (5, 2) = 638, (5, 3) = 619, (5, 4) = 645, (5, 5) = 0, (5, 6) = 545, (5, 7) = 653, (5, 8) = 482, (5, 9) = 869, (5, 10) = 742, (5, 11) = 883, (6, 1) = 318, (6, 2) = 457, (6, 3) = 346, (6, 4) = 317, (6, 5) = 545, (6, 6) = 0, (6, 7) = 440, (6, 8) = 573, (6, 9) = 340, (6, 10) = 459, (6, 11) = 638, (7, 1) = 132, (7, 2) = 184, (7, 3) = 45, (7, 4) = 128, (7, 5) = 653, (7, 6) = 440, (7, 7) = 0, (7, 8) = 501, (7, 9) = 365, (7, 10) = 318, (7, 11) = 406, (8, 1) = 633, (8, 2) = 746, (8, 3) = 621, (8, 4) = 648, (8, 5) = 482, (8, 6) = 573, (8, 7) = 501, (8, 8) = 0, (8, 9) = 1045, (8, 10) = 850, (8, 11) = 1087, (9, 1) = 478, (9, 2) = 278, (9, 3) = 323, (9, 4) = 379, (9, 5) = 869, (9, 6) = 340, (9, 7) = 365, (9, 8) = 1045, (9, 9) = 0, (9, 10) = 243, (9, 11) = 24, (10, 1) = 346, (10, 2) = 314, (10, 3) = 212, (10, 4) = 297, (10, 5) = 742, (10, 6) = 459, (10, 7) = 318, (10, 8) = 850, (10, 9) = 243, (10, 10) = 0, (10, 11) = 239, (11, 1) = 585, (11, 2) = 321, (11, 3) = 376, (11, 4) = 337, (11, 5) = 883, (11, 6) = 638, (11, 7) = 406, (11, 8) = 1087, (11, 9) = 24, (11, 10) = 239, (11, 11) = 0})

#Check and ENFORCE triangle inequality:
for T in combinat:-choose(n,3) do
    if (S:= Cost[T[1],T[2]] + Cost[T[2],T[3]]) < Cost[T[1],T[3]] then
        printf("Triangle inequality violation at %a corrected.\n", T);
        Cost[T[1],T[3]]:= S #Enforcement
    fi
od:  

Triangle inequality violation at [1, 3, 4] corrected.
Triangle inequality violation at [1, 7, 11] corrected.
Triangle inequality violation at [1, 9, 11] corrected.
Triangle inequality violation at [2, 7, 8] corrected.

Triangle inequality violation at [2, 9, 11] corrected.
Triangle inequality violation at [3, 4, 6] corrected.
Triangle inequality violation at [3, 4, 11] corrected.
Triangle inequality violation at [3, 7, 8] corrected.
Triangle inequality violation at [3, 9, 11] corrected.
Triangle inequality violation at [4, 7, 8] corrected.
Triangle inequality violation at [6, 9, 11] corrected.
Triangle inequality violation at [7, 9, 11] corrected.
Triangle inequality violation at [8, 9, 11] corrected.

#View corrected matrix:
Cost;

Matrix(11, 11, {(1, 1) = 0, (1, 2) = 376, (1, 3) = 246, (1, 4) = 259, (1, 5) = 599, (1, 6) = 318, (1, 7) = 132, (1, 8) = 633, (1, 9) = 478, (1, 10) = 346, (1, 11) = 502, (2, 1) = 376, (2, 2) = 0, (2, 3) = 326, (2, 4) = 240, (2, 5) = 638, (2, 6) = 457, (2, 7) = 184, (2, 8) = 685, (2, 9) = 278, (2, 10) = 314, (2, 11) = 302, (3, 1) = 246, (3, 2) = 326, (3, 3) = 0, (3, 4) = 13, (3, 5) = 619, (3, 6) = 330, (3, 7) = 45, (3, 8) = 546, (3, 9) = 323, (3, 10) = 212, (3, 11) = 347, (4, 1) = 259, (4, 2) = 240, (4, 3) = 13, (4, 4) = 0, (4, 5) = 645, (4, 6) = 317, (4, 7) = 128, (4, 8) = 629, (4, 9) = 379, (4, 10) = 297, (4, 11) = 337, (5, 1) = 599, (5, 2) = 638, (5, 3) = 619, (5, 4) = 645, (5, 5) = 0, (5, 6) = 545, (5, 7) = 653, (5, 8) = 482, (5, 9) = 869, (5, 10) = 742, (5, 11) = 883, (6, 1) = 318, (6, 2) = 457, (6, 3) = 330, (6, 4) = 317, (6, 5) = 545, (6, 6) = 0, (6, 7) = 440, (6, 8) = 573, (6, 9) = 340, (6, 10) = 459, (6, 11) = 364, (7, 1) = 132, (7, 2) = 184, (7, 3) = 45, (7, 4) = 128, (7, 5) = 653, (7, 6) = 440, (7, 7) = 0, (7, 8) = 501, (7, 9) = 365, (7, 10) = 318, (7, 11) = 389, (8, 1) = 633, (8, 2) = 685, (8, 3) = 546, (8, 4) = 629, (8, 5) = 482, (8, 6) = 573, (8, 7) = 501, (8, 8) = 0, (8, 9) = 1045, (8, 10) = 850, (8, 11) = 1069, (9, 1) = 478, (9, 2) = 278, (9, 3) = 323, (9, 4) = 379, (9, 5) = 869, (9, 6) = 340, (9, 7) = 365, (9, 8) = 1045, (9, 9) = 0, (9, 10) = 243, (9, 11) = 24, (10, 1) = 346, (10, 2) = 314, (10, 3) = 212, (10, 4) = 297, (10, 5) = 742, (10, 6) = 459, (10, 7) = 318, (10, 8) = 850, (10, 9) = 243, (10, 10) = 0, (10, 11) = 239, (11, 1) = 502, (11, 2) = 302, (11, 3) = 347, (11, 4) = 337, (11, 5) = 883, (11, 6) = 364, (11, 7) = 389, (11, 8) = 1069, (11, 9) = 24, (11, 10) = 239, (11, 11) = 0})

#
# We define the decision variables thus:
#    x[i,j,v] = 1, if truckload v goes from destination i to j;
#             = 0, otherwise.
#

#Objective:   (Eq 2.2)
#   minimize:
Obj:= add(add(add(Cost[i,j]*x[i,j,v], i= N), j= N), v= V)
:

###############
# Constraints #
###############
Con:= table()
:
# Each truckload does not exceed the capacity: (Eq 2.4)
#
Con[capacity]:=
    seq(add(d[j]*add(x[i,j,v], i= N minus {j}), j= N) <= K, v= V)
:

# Each delivery is made exactly once:  (Eq 2.3)
#
Con[delivery]:=
    seq(add(add(x[i,j,v], j= N minus {i}), v= V) = 1, i= C)
:

# Trucks start at warehouse:  (Eq 2.5)
#
Con[start]:=
    seq(add(x[1,j,v], j= C) = 1, v= V)
:

# Trucks return to warehouse: (Carl's idea)
#
Con[`return`]:=
    seq(add(x[i,1,v], i= C) = 1, v= V)
:

# Continuity: (Eq 2.6)
#
Con[continuity]:=
    seq(
        seq(
            add(x[i,k,v], i= N minus {k})
                = add(x[k,j,v], j= N minus {k}),
            k= C
        ),
        v= V
    )
:

# Avoid subtours:  (VV's idea) (u's are artificial variables.)
#
Con[antisubtour]:=
    seq(
        seq(
            seq(
                u[i,v]-u[j,v]+n*x[i,j,v] <= n-1,
                i= N minus {1,j}
            ),
            j= C
        ),
        v= V
    )
:

BV:= binaryvariables
    = {seq(seq(seq(x[i,j,v], i= N minus {j}), j= N), v= V)}
:

infolevel[Optimization]:= 5:
Sol:= CodeTools:-Usage(
    Optimization:-LPSolve(
        Obj, [entries(Con, nolist)],
        BV, assume= nonnegative, feasibilitytolerance= 1e-5
    )
):

LPSolve: calling ILP solver

LPSolve: number of problem variables 240
LPSolve: number of general linear constraints 216
LPSolve: feasibility tolerance set to 0.1e-4
LPSolve: integer tolerance set to 0.10e-4
LPSolve: infinite bound set to 0.10e21
LPSolve: iteration limit set to 2280
LPSolve: depth limit set to 360
LPSolve: node limit set to 0

memory used=5.86MiB, alloc change=2.76MiB, cpu time=29.52s, real time=29.68s, gc time=0ns

Stops:= {lhs~(select(e-> rhs(e)::1, Sol[2]))[]}:

Tour:= table():
for X in Stops do
    Tour[op(3,X)][op(1,X)]:= op(2,X)
od:

for v in V do
    A:= Array(1..1, [1]);
    for i do A,= (r:=  Tour[v][A[i]]) until r=1;
    Tour[v]:= seq(A)
od:

Routes:= [entries(Tour)];

[[1, 2, 10, 1], [1, 6, 9, 11, 4, 3, 7, 8, 5, 1]]

COST = Sol[1];

COST = 3695

 


 

Download VehicleRoutingILP.mw

Do you mean something like this?

unapply([for k to 3 do x^k od], x);

This uses syntax new to Maple 2019. The example above would be better done with seq, but there are more-elaborate examples where a do-loop (or for-loop) works and seq doesn't. There are also cases where seq is the less-efficient option.

I would've given the same answer as Joe and Acer if I'd gotten to it in time. I'll add that the same level of type checking can also be applied to the return value of a procedure:

P:= proc()::posint; 0 end proc:  

P(); #error when assertlevel is set to 2

But mostly I want to draw your attention to the distinction between types and properties because I think that you may be referring to the latter. Types are inherent to data structures (i.e., anything that you can assign to a variable): Given any structure and any type ttype(x,t) should return true or false; nebulous or ambiguous cases should be avoided when coding t.

Properties, on the other hand, can be given to symbolic variables (i.e., names) with assume and related commands. Some "property-aware" commands will respect the properties, but most commands will ignore them.

Maple's type system is highly developed, but the property system is nascent and has been so AFAIK since created by Gaston Gonnet, et al, many years ago. My guess is that making significant improvements to the system would require redesign of Maple from the ground up. I believe that the strongly typed "free" CAS named Axiom has something like this implemented. Maybe it'd be possible to get Maple and Axiom to interact?

Use:

cat("", V, ".xls")

The plot3d command needs more quotes. This works:

plot3d(['f(x,y,'g(x, y)')'], x=0..1,y=0..1)

Use

f~([$0..5])

If you put immediately before an operator, it makes it inert. So, you can do a %/ b * c, and it'll display like you want. The inertness will remain until you use the value command.

Unfortunately, you do need to use some special syntax to get things to display as you want.

If an appears in the first column of input, it's a shell escape, equivalent to a system command. See ?system.

@digerdiga You wrote:

  • In a module the structure is somewhat different compared to a procedure.

Yes, modules have two types of local variables: locals and exports. The only difference between them is that exports can be accessed from outside the module (no distinction is made between read access and write access). By far the most important thing about modules is that both of these types of variables retain their values between uses of the module, up until a restart. In this way, these variables are like globals, but (1) they don't corrupt the namespace, and (2) there can simultaneously exist in memory multiple copies of essentially the same module with different values of the variables. None of this applies to the variables that are local to the procedures inside a module; those act just like any procedure local.

  • but the following two lines I just asked about are not directly run when calling Rand(), or?

The line ModuleApply:= proc ... end proc (as well as all other assignment statements in the localexport, and "body" sections of the module) are executed at the time that the module is defined, not when the module is used.

The ModuleApply procedure itself (as opposed to its assignment) is directly run every time that you call Rand(). The ModuleLoad procedure is only run if and when the module is loaded from a library, or it's directly invoked via ModuleLoad().

  • I mean in order for the procedure ModuleApply to run it needs the seed, but this seed is only given later in the export section.

You are mistaking the "body" of the module as being part of the export section. Note the semicolon (I always left-justify this semicolon) that ends the export section. Everything after that semicolon up until end module is the module body. These body statements are executed when (and only when) the module's code is evaluated (like when it's typed into a worksheet). In this case, the only body statement is seed:= 1. Thus, seed has a value before you call Rand().

  • So when calling Rand() does the module read everything till the end and then afterwards calls the ModuleApply?

No, that reading (and assigning) is only done once, when the module's code is evaluated.

  • Do you have an example of how ModuleLoad could be read from a Library? Not sure I understand you right.

Whether the ModuleLoad itself is read from a library is irrelevant. The key thing is that when the module itself is read (or "loaded") from a library (aka "repository"), its ModuleLoad procedure (if it has one) gets executed.

Any Maple data structure that can be assigned to a name can be saved in a disk file. There are three types of such files: plaintext, ".m", and libraries. Most large collections of code (including the vast majority of ordinary Maple commands) are stored in libraries.

The purpose of my ModuleLoad is to make sure that seed has a value when it's needed if the code is loaded from a library. The only way that it would get into a library would be if someone explicitly put it there. So, the ModuleLoad in this example does nothing. It's only there so that the code will work if someone decides to put it into a library.

The three types of modules: There are three types of modules: records, just plain modules, and objects. Records are vastly simplified modules that are only used as fancy containers, similar to tables. They only have exports, not locals or bodies. (See ?Record). On the other hand, objects are vastly more complicated than plain modules. While everything that I wrote above is also true of objects, this post couldn't be considered a primer on them because it barely scratches the surface.

The parameter sequence of a procedure can be extracted as its first operand, op(1, ...). This includes type declarations and default values. If those appear, they can be stripped off, also with op(1, ...) (just a coincidence that the same op command does all 3 things). So, the following procedure returns as a list just the parameter names, regardless of how many or whether the procedure is given in proc form or arrow form:

Params:= (f::procedure)-> op~(1, op~(1, [op(1, eval(f))])):

Warning: This ignores the possibility that the parameter names have been assigned values as globals.
Edit: That assigned-global problem is solved in a Reply below.

A complete implementation of your Iter follows in the first Reply.

If the objective function is passed as a procedure (as suggested in the last paragraph of Acer's Reply), then you can print whatever you want from inside that procedure. Doing this on Kitonum's example, we see that the actual number of iterations is much higher than are probably expected:

infolevel[ShowIter]:= 1:
Optimization:-Minimize(
    proc(x,y)
    local r:= exp((x-1)^2+(y-2)^2);
        userinfo(1, ShowIter, NoName, "x=", x, "y=", y, "val"= r);
        r
    end proc
);
x= HFloat(0.0) y= HFloat(0.0) val = HFloat(148.4131591025766)
x= 0. y= HFloat(0.0) val = 148.41315910257660342111558
x= 0.50000000000000000000000000e-8 y= HFloat(0.0) val = 148.41315761844502352633642
x= -0.50000000000000000000000000e-8 y= HFloat(0.0) val = 148.41316058670820557786861
x= HFloat(0.0) y= 0. val = 148.41315910257660342111558
x= HFloat(0.0) y= 0.50000000000000000000000000e-8 val = 148.41315613431345476254404
x= HFloat(0.0) y= -0.50000000000000000000000000e-8 val = 148.41316207083981886560872
x= HFloat(296.826318205153) y= HFloat(593.652636410306) val = HFloat(HFloat(infinity))
x= 296.8263182051530 y= HFloat(593.652636410306) val = 0.33278359824732988270238406e190033
x= 296.82631821015300000000000 y= HFloat(593.652636410306) val = 0.33278458271025232621688077e190033
x= 296.82631820015300000000000 y= HFloat(593.652636410306) val = 0.33278261378731975193693085e190033
x= HFloat(296.826318205153) y= 593.6526364103059 val = 0.33278359822409659209765980e190033
x= HFloat(296.826318205153) y= 593.65263641530590000000000 val = 0.33278556715285363807149492e190033
x= HFloat(296.826318205153) y= 593.65263640530590000000000 val = 0.33278162930698874720167114e190033
x= HFloat(148.4131591025765) y= HFloat(296.826318205153) val = HFloat(HFloat(infinity))
x= 148.4131591025765 y= HFloat(296.826318205153) val = 0.30483482837354069869543186e47188
x= 148.41315910757650000000000 y= HFloat(296.826318205153) val = 0.30483527774052246962191513e47188
x= 148.41315909757650000000000 y= HFloat(296.826318205153) val = 0.30483437900722136858926819e47188
x= HFloat(148.4131591025765) y= 296.8263182051530 val = 0.30483482837731537856171930e47188
x= HFloat(148.4131591025765) y= 296.82631821015300000000000 val = 0.30483572711194135047156588e47188
x= HFloat(148.4131591025765) y= 296.82631820015300000000000 val = 0.30483392964533912419576772e47188
x= HFloat(74.20657955128824) y= HFloat(148.4131591025765) val = HFloat(HFloat(infinity))
x= 74.20657955128824 y= HFloat(148.4131591025765) val = 0.23018676697883155390839898e11638
x= 74.206579556288240000000000 y= HFloat(148.4131591025765) val = 0.23018693549075192545878372e11638
x= 74.206579546288240000000000 y= HFloat(148.4131591025765) val = 0.23018659846703455563500392e11638
x= HFloat(74.20657955128824) y= 148.4131591025765 val = 0.23018676697987969766940068e11638
x= HFloat(74.20657955128824) y= 148.41315910757650000000000 val = 0.23018710400384379840854317e11638
x= HFloat(74.20657955128824) y= 148.41315909757650000000000 val = 0.23018642995640905551112493e11638
x= HFloat(37.10328977564412) y= HFloat(74.20657955128824) val = HFloat(HFloat(infinity))
x= 37.10328977564412 y= HFloat(74.20657955128824) val = 0.25139367418419102481734851e2831
x= 37.103289780644120000000000 y= HFloat(74.20657955128824) val = 0.25139376494559408338451085e2831
x= 37.103289770644120000000000 y= HFloat(74.20657955128824) val = 0.25139358342282074666630316e2831
x= HFloat(37.10328977564412) y= 74.20657955128824 val = 0.25139367418410752434162267e2831
x= HFloat(37.10328977564412) y= 74.206579556288240000000000 val = 0.25139385570694640298907617e2831
x= HFloat(37.10328977564412) y= 74.206579546288240000000000 val = 0.25139349266139972964958592e2831
x= HFloat(18.55164488782206) y= HFloat(37.10328977564412) val = HFloat(HFloat(infinity))
x= 18.55164488782206 y= HFloat(37.10328977564412) val = 0.87964587087268480177328466e669
x= 18.551644892822060000000000 y= HFloat(37.10328977564412) val = 0.87964602526501789892447776e669
x= 18.551644882822060000000000 y= HFloat(37.10328977564412) val = 0.87964571648037884699604237e669
x= HFloat(18.55164488782206) y= 37.10328977564412 val = 0.87964587087261069345987837e669
x= HFloat(18.55164488782206) y= 37.103289780644120000000000 val = 0.87964617965730396414151683e669
x= HFloat(18.55164488782206) y= 37.103289770644120000000000 val = 0.87964556208802586032716252e669
x= HFloat(9.27582244391103) y= HFloat(18.55164488782206) val = HFloat(5.278235480625041e148)
x= 9.275822443911030 y= HFloat(18.55164488782206) val = 0.52782354806250521273629132e149
x= 9.2758224489110300000000000 y= HFloat(18.55164488782206) val = 0.52782359174424668824786704e149
x= 9.2758224389110300000000000 y= HFloat(18.55164488782206) val = 0.52782350438076737863910726e149
x= HFloat(9.27582244391103) y= 18.55164488782206 val = 0.52782354806249464175529183e149
x= HFloat(9.27582244391103) y= 18.551644892822060000000000 val = 0.52782363542598119460461830e149
x= HFloat(9.27582244391103) y= 18.551644882822060000000000 val = 0.52782346069902257538999975e149
x= HFloat(0.009881287988267002) y= HFloat(0.019762575976534004) val = HFloat(134.51494651573813)
x= 0.9881287988267002e-2 y= HFloat(0.019762575976534004) val = 134.51494651573822415377608
x= 0.98812929882670020000000000e-2 y= HFloat(0.019762575976534004) val = 134.51494518388057820523680
x= 0.98812829882670020000000000e-2 y= HFloat(0.019762575976534004) val = 134.51494784759589001503456
x= HFloat(0.009881287988267002) y= 0.1976257597653400e-1 val = 134.51494651573822623946518
x= HFloat(0.009881287988267002) y= 0.19762580976534000000000000e-1 val = 134.51494385202294416648466
x= HFloat(0.009881287988267002) y= 0.19762570976534000000000000e-1 val = 134.51494917945356778608051
x= HFloat(0.09630756341701092) y= HFloat(0.19261512683402177) val = HFloat(59.34097821750855)
x= 0.9630756341701092e-1 y= HFloat(0.19261512683402177) val = 59.340978217508519414268186
x= 0.96307568417010920000000000e-1 y= HFloat(0.19261512683402177) val = 59.340977681248591374879303
x= 0.96307558417010920000000000e-1 y= HFloat(0.19261512683402177) val = 59.340978753768455266846425
x= HFloat(0.09630756341701092) y= .1926151268340218 val = 59.340978217508511627773979
x= HFloat(0.09630756341701092) y= .19261513183402180000000000 val = 59.340977144988658911612276
x= HFloat(0.09630756341701092) y= .19261512183402180000000000 val = 59.340979290028386695546375
x= HFloat(0.15456181325214174) y= HFloat(0.3091236265042795) val = HFloat(35.65283397381811)
x= .1545618132521417 y= HFloat(0.3091236265042795) val = 35.652833973818131700218521
x= .15456181825214170000000000 y= HFloat(0.3091236265042795) val = 35.652833672395460793238144
x= .15456180825214170000000000 y= HFloat(0.3091236265042795) val = 35.652834275240806938182978
x= HFloat(0.15456181325214174) y= .3091236265042795 val = 35.652833973818128083146444
x= HFloat(0.15456181325214174) y= .30912363150427950000000000 val = 35.652833370972787926205873
x= HFloat(0.15456181325214174) y= .30912362150427950000000000 val = 35.652834576663480216098241
x= HFloat(0.22933338820309085) y= HFloat(0.4586667764061795) val = HFloat(19.48480892911679)
x= .2293333882030908 y= HFloat(0.4586667764061795) val = 19.484808929116789468211163
x= .22933339320309080000000000 y= HFloat(0.4586667764061795) val = 19.484808778953873744833124
x= .22933338320309080000000000 y= HFloat(0.4586667764061795) val = 19.484809079279707323085112
x= HFloat(0.22933338820309085) y= .4586667764061795 val = 19.484808929116786615115743
x= HFloat(0.22933338820309085) y= .45866678140617950000000000 val = 19.484808628790955838494530
x= HFloat(0.22933338820309085) y= .45866677140617950000000000 val = 19.484809229442622994999254
x= HFloat(0.30356276727546316) y= HFloat(0.6071255345509279) val = HFloat(11.303632063533023)
x= .3035627672754632 y= HFloat(0.6071255345509279) val = 11.303632063533021787011651
x= .30356277227546320000000000 y= HFloat(0.6071255345509279) val = 11.303631984810320003096722
x= .30356276227546320000000000 y= HFloat(0.6071255345509279) val = 11.303632142255724684362392
x= HFloat(0.30356276727546316) y= .6071255345509279 val = 11.303632063533021771267111
x= HFloat(0.30356276727546316) y= .60712553955092790000000000 val = 11.303631906087618469100828
x= HFloat(0.30356276727546316) y= .60712552955092790000000000 val = 11.303632220978427831631836
x= HFloat(0.3853589567012047) y= HFloat(0.7707179134024081) val = HFloat(6.612210797744427)
x= .3853589567012047 y= HFloat(0.7707179134024081) val = 6.6122107977444285321483808
x= .38535896170120470000000000 y= HFloat(0.7707179134024081) val = 6.6122107571030673899809423
x= .38535895170120470000000000 y= HFloat(0.7707179134024081) val = 6.6122108383857902547248470
x= HFloat(0.3853589567012047) y= .7707179134024081 val = 6.6122107977444293287190656
x= HFloat(0.3853589567012047) y= .77071791840240810000000000 val = 6.6122107164617071288773094
x= HFloat(0.3853589567012047) y= .77071790840240810000000000 val = 6.6122108790271528583653134
x= HFloat(0.4726539000765073) y= HFloat(0.9453078001530171) val = HFloat(4.016735650512733)
x= .4726539000765073 y= HFloat(0.9453078001530171) val = 4.0167356505127340463505823
x= .47265390507650730000000000 y= HFloat(0.9453078001530171) val = 4.0167356293306354054050219
x= .47265389507650730000000000 y= HFloat(0.9453078001530171) val = 4.0167356716948329998358971
x= HFloat(0.4726539000765073) y= .9453078001530171 val = 4.0167356505127337667468781
x= HFloat(0.4726539000765073) y= .94530780515301710000000000 val = 4.0167356081485364961404407
x= HFloat(0.4726539000765073) y= .94530779515301710000000000 val = 4.0167356928769316850019857
x= HFloat(0.5676775476700361) y= HFloat(1.135355095340071) val = HFloat(2.545974575281497)
x= .5676775476700361 y= HFloat(1.135355095340071) val = 2.5459745752814969298092981
x= .56767755267003610000000000 y= HFloat(1.135355095340071) val = 2.5459745642746772976967882
x= .56767754267003610000000000 y= HFloat(1.135355095340071) val = 2.5459745862883167368054897
x= HFloat(0.5676775476700361) y= 1.135355095340071 val = 2.5459745752814972225907026
x= HFloat(0.5676775476700361) y= 1.1353551003400710000000000 val = 2.5459745532678579423012381
x= HFloat(0.5676775476700361) y= 1.1353550903400710000000000 val = 2.5459745972951368205187076
x= HFloat(0.6704666849677644) y= HFloat(1.3409333699355328) val = HFloat(1.7210955376872163)
x= .6704666849677644 y= HFloat(1.3409333699355328) val = 1.7210955376872164962165431
x= .67046668996776440000000000 y= HFloat(1.3409333699355328) val = 1.7210955320156333683762451
x= .67046667996776440000000000 y= HFloat(1.3409333699355328) val = 1.7210955433587997288013740
x= HFloat(0.6704666849677644) y= 1.340933369935533 val = 1.7210955376872160379526221
x= HFloat(0.6704666849677644) y= 1.3409333749355330000000000 val = 1.7210955263440497579344690
x= HFloat(0.6704666849677644) y= 1.3409333649355330000000000 val = 1.7210955490303824787845764
x= HFloat(0.7797359062867295) y= HFloat(1.559471812573459) val = HFloat(1.2745349346035728)
x= .7797359062867295 y= HFloat(1.559471812573459) val = 1.2745349346035727507930728
x= .77973591128672950000000000 y= HFloat(1.559471812573459) val = 1.2745349317962299629846458
x= .77973590128672950000000000 y= HFloat(1.559471812573459) val = 1.2745349374109156085118147
x= HFloat(0.7797359062867295) y= 1.559471812573459 val = 1.2745349346035727687600668
x= HFloat(0.7797359062867295) y= 1.5594718175734590000000000 val = 1.2745349289888871674634077
x= HFloat(0.7797359062867295) y= 1.5594718075734590000000000 val = 1.2745349402182584585177456
x= HFloat(0.8868345174144543) y= HFloat(1.7736690348289055) val = HFloat(1.06612665525306)
x= .8868345174144543 y= HFloat(1.7736690348289055) val = 1.0661266552530599880945141
x= .88683452241445430000000000 y= HFloat(1.7736690348289055) val = 1.0661266540465726410400816
x= .88683451241445430000000000 y= HFloat(1.7736690348289055) val = 1.0661266564595473898206066
x= HFloat(0.8868345174144543) y= 1.773669034828905 val = 1.0661266552530602214291723
x= HFloat(0.8868345174144543) y= 1.7736690398289050000000000 val = 1.0661266528400855020324293
x= HFloat(0.8868345174144543) y= 1.7736690298289050000000000 val = 1.0661266576660349995935570
x= HFloat(0.9675495640300781) y= HFloat(1.935099128060163) val = HFloat(1.005279039255056)
x= .9675495640300781 y= HFloat(1.935099128060163) val = 1.0052790392550560115855033
x= .96754956903007810000000000 y= HFloat(1.935099128060163) val = 1.0052790389288386058179008
x= .96754955903007810000000000 y= HFloat(1.935099128060163) val = 1.0052790395812734677229168
x= HFloat(0.9675495640300781) y= 1.935099128060163 val = 1.0052790392550560181098519
x= HFloat(0.9675495640300781) y= 1.9350991330601630000000000 val = 1.0052790386026211815485982
x= HFloat(0.9675495640300781) y= 1.9350991230601630000000000 val = 1.0052790399074909053584935
x= HFloat(0.9974615874297784) y= HFloat(1.9949231748595513) val = HFloat(1.0000322182108787)
x= .9974615874297784 y= HFloat(1.9949231748595513) val = 1.0000322182108787599366068
x= .99746159242977840000000000 y= HFloat(1.9949231748595513) val = 1.0000322181854938414044036
x= .99746158242977840000000000 y= HFloat(1.9949231748595513) val = 1.0000322182362637284710653
x= HFloat(0.9974615874297784) y= 1.994923174859551 val = 1.0000322182108787626883347
x= HFloat(0.9974615874297784) y= 1.9949231798595510000000000 val = 1.0000322181601089006237092
x= HFloat(0.9974615874297784) y= 1.9949231698595510000000000 val = 1.0000322182616486747571486
x= HFloat(0.9999856333913124) y= HFloat(1.9999712667826288) val = HFloat(1.0000000010319972)
x= .9999856333913124 y= HFloat(1.9999712667826288) val = 1.0000000010319972262146130
x= .99998563839131240000000000 y= HFloat(1.9999712667826288) val = 1.0000000010318535851275888
x= .99998562839131240000000000 y= HFloat(1.9999712667826288) val = 1.0000000010321409173016373
x= HFloat(0.9999856333913124) y= 1.999971266782629 val = 1.0000000010319972262039243
x= HFloat(0.9999856333913124) y= 1.9999712717826290000000000 val = 1.0000000010317099190299178
x= HFloat(0.9999856333913124) y= 1.9999712617826290000000000 val = 1.0000000010322845833779309
x= HFloat(0.9999999995345299) y= HFloat(1.9999999990690573) val = HFloat(1.0)
x= .9999999995345299 y= HFloat(1.9999999990690573) val = 1.0000000000000000010833166
x= 1.0000000045345299000000000 y= HFloat(1.9999999990690573) val = 1.0000000000000000214286156
x= .99999999453452990000000000 y= HFloat(1.9999999990690573) val = 1.0000000000000000307380176
x= HFloat(0.9999999995345299) y= 1.999999999069057 val = 1.0000000000000000010833173
x= HFloat(0.9999999995345299) y= 2.0000000040690570000000000 val = 1.0000000000000000167738873
x= HFloat(0.9999999995345299) y= 1.9999999940690570000000000 val = 1.0000000000000000353927473
x= HFloat(0.9999999999999993) y= HFloat(2.0000000000000004) val = HFloat(1.0)
x= .9999999999999993 y= HFloat(2.0000000000000004) val = 1.
x= 1.0000000049999993000000000 y= HFloat(2.0000000000000004) val = 1.0000000000000000249999930
x= .99999999499999930000000000 y= HFloat(2.0000000000000004) val = 1.0000000000000000250000070
x= HFloat(0.9999999999999993) y= 2.000000000000000 val = 1.
x= HFloat(0.9999999999999993) y= 2.0000000050000000000000000 val = 1.0000000000000000250000000
x= HFloat(0.9999999999999993) y= 1.9999999950000000000000000 val = 1.0000000000000000250000000
                   [    [0.999999999999999]]
                   [1., [                 ]]
                   [    [2.00000000000000 ]]

There is a possibility that using procedure-form input increases the number of iterations because if it can't "see" the function symbolically, then it can't get symbolic partial derivatives (or gradient). Finding the zeros of the gradient of the function in this example is trivial once the gradient is known symbolically.

I saw that you just asked a Question about the derivative of Eq. 6. That Question got deleted, I guess because it's too similar to this one.

Anyway, in that Question, you took the derivative of both sides of Eq6 to get Eq6_2, and you asked for a plot of its solution. This operation is redundant or superfluous because taking the derivative of both sides of an ODE doesn't fundamentally change its solution (provided that an appropriate initial condition is also added). I suspect that what you really want is a plot of the derivative of f(eta), which is a different thing than the solution of the derivative of an equation that contains f(eta). To plot the derivative of f(eta), all that you need to do is change f(eta) to diff(f(eta), eta) in the odeplot command. No other change is needed. Doing that, and also removing the plot title, I get

Here is a complete example of a pseudo-random number generator (PRNG), just intended as an educational example so that you can see how they work. It's surprisingly short and simple. It generates integers between 0 and 100, inclusive. The only difference between this and an actual PNRG is that the actual one would use much larger integers and probably something a little bit more complicated than a linear modular congruence. The way that the Randomize in this module works is almost identical to Maple's own randomize. Note that floats are generated from integers. That's the way it actually works in Maple.

Rand:= module()
local 
    A:= 47, B:= 53, M:= 101, #Note that these are prime.
    seed,
    ModuleApply:= proc() seed:= irem(A*seed+B, M) end proc,
    ModuleLoad:= proc() seed:= 1 end proc
;
export 
    Randomize:= proc(s::posint:= 0) seed:= `if`(s=0, Value(Time()), s) end proc,
    Float:= ()-> evalf(Rand()/M)
;
    seed:= 1
end module
:    
'Rand()' $ 9;
                100, 6, 32, 42, 7, 79, 29, 2, 46
Seed1:= Rand:-Randomize();
                     Seed1 := 1577955051735
'Rand()' $ 9;
               34, 35, 82, 69, 64, 31, 96, 20, 84
'Rand()' $ 9;
               62, 38, 21, 30, 49, 33, 89, 95, 74
Rand:-Randomize(Seed1): #Restart sequence from Seed1.
'Rand()' $ 9;
               34, 35, 82, 69, 64, 31, 96, 20, 84

 

These 3 boundary conditions inside the loop are nonsense:

D(F[k])(last) = 0, H[k](last) = 0, T[k](last) = 0

They're nonsense because last is the number of loop iterations; it's not an endpoint of a real interval.

Yes, compiling on each node helps. Here is an example. The part about wasting time is just so that the task is nontrivial and thus I can measure the efficiency of the multiprocessing and know for sure that multiple processors were used.
 

restart:

Sq:= proc(x::integer[4])::integer[4];
local i::integer[4]:= 0, k::integer[4];
    #This loop is only to waste time (approx 4 secs):
    for k to 2\000\000\000 do i:= i+1 od;
    x^2
end proc:

Sq_comp:= Compiler:-Compile(Sq):

#Measure time for one simple run:
rt:= CodeTools:-Usage(Sq_comp(1), output= realtime);

memory used=1.30KiB, alloc change=0 bytes, cpu time=4.34s, real time=4.34s, gc time=0ns

4.342

nc:= kernelopts(numcpus);

8

m:= 3: #can be any small posint

S:= [$1..m*nc]: #input data for distributed processing

#
#Here I show that this naive technique generates an error:
st:= time[real]():
R:= Grid:-Map(Sq_comp, S):
rt2:= time[real]() - st:

Error, (in Grid:-Map) invalid input: op expects 1 or 2 arguments, but received 0

#######################
# Try a different way #
#######################

#Save definitions so that I can restart:
currentdir("C:/Users/carlj/desktop"):
save Sq, rt, nc, m, S, "GridTest.mpl"; #not Sq_comp

#
restart:

#
currentdir("C:/Users/carlj/desktop"):
read "GridTest.mpl":

#
#Correct that error by compiling on each node:
Grid:-Set(Sq);
Grid:-Set(Sq_comp= 'Compiler:-Compile'(Sq));

#
#Attempt distributed processing again:
st:= time[real]():
R:= Grid:-Map(Sq_comp, S):
rt2:= time[real]() - st:

#
#Check accuracy:
evalb(R = S^~2);

true

#
#Check efficiency:
E:= m*rt/rt2;

.6957961647

 


 

Download GridComp.mw

There may also be a way to do it without compiling on each node. I don't know.

Like this:

plots:-textplot(
   [seq(
       [X[k], Y[k], sprintf("%3d", k), color= COLOR(HSV, .85*(k-1)/(N-1), 1, .8)],
       k= 1..N
   )],
   font= ["times", "bold", 9], axes= none, size= [700$2]
);

First 118 119 120 121 122 123 124 Last Page 120 of 395