Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 175 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Some speed-ups...

tomleslie 13876
July 24 2018
2 0

on Maple 2018.1, the worksheet supply

  1. fails to evaluate the sum assigned to 'P' with the error "too many levels of recursion"
  2. with the above failure, executing the complete worksheet takes about 125 secs (on my machine)

If you convert the sum() and int() commands to the inert forms 'Sum()' and 'Int()', with evalf() commands where necessary, then everything in the worksheet executes correctly, and takes about 25secs on my machine (see the attached).

You state that you are using this in a 1000-step loop. Whether or not the approached suggested in the attached, will always work for your loop depends a lot which variables you are looping over and the values they are taking, so I'd try it over a loop with a few iterations first.

Again depending, on which variables you are looping over, it is always advisable to shift as many calculations "outside" the loop, and then just evaluate the calculations within the loop. For example, I doubt if it is strictly necessary to repeatedly compute the symbolic differentials within the loop


 

NULL

> 

restart
> 

tstart := time()
> 

`λλ`[1] := .1111111; `λλ`[2] := 2222222; `αα` := 1.51222222
> 

NULL
> 

P := Sum(Sum((t+1)*`αα`^2*(1-`αα`)^(t+T)/(t+1+`λλ`[1]*(T+1)/`λλ`[2]), t = 0 .. infinity), T = 0 .. infinity); evalf(%)

Sum(Sum(2.286816043*(t+1)*(-.51222222)^(t+T)/(t+1.000000050+0.5000000000e-7*T), t = 0 .. infinity), T = 0 .. infinity)

  

.9999999598

 (1)
> 

r[1] := exp(-C*(Int(lambda[1]*alpha^2*exp(lambda[1]*Z)/((exp(lambda[1]*Z)-1+alpha)^2*(exp(lambda[2]*Z)-1+alpha)), Z = 0 .. infinity)))

exp(-C*(Int(lambda[1]*alpha^2*exp(lambda[1]*Z)/((exp(lambda[1]*Z)-1+alpha)^2*(exp(lambda[2]*Z)-1+alpha)), Z = 0 .. infinity)))

 (2)
> 

``
> 

R[1] := diff(r[1], lambda[1]); R[1, 1] := diff(r[1], lambda[1], lambda[1]); R[1, 2] := diff(r[1], lambda[1], lambda[2]); R[1, 3] := diff(r[1], lambda[1], alpha); R[2] := diff(r[1], lambda[2]); R[2, 1] := diff(r[1], lambda[2], lambda[1]); R[2, 2] := diff(r[1], lambda[2], lambda[2]); R[2, 3] := diff(r[1], lambda[2], alpha); R[3] := diff(r[1], alpha); R[3, 1] := diff(r[1], alpha, lambda[1]); R[3, 2] := diff(r[1], alpha, lambda[2]); R[3, 3] := diff(r[1], alpha, alpha)

 

> 

lambda[1] := 1.117480; lambda[2] := .124861; alpha := 7.629478
> 

NULL
> 

C := -2
> 

r_r[1] := evalf(r[1]); R_R[1] := evalf(R[1]); R_R[1, 1] := evalf(R[1, 1]); R_R[1, 2] := evalf(R[1, 2]); R_R[1, 3] := evalf(R[1, 3]); R_R[2] := evalf(R[2]); R_R[2, 1] := evalf(R[2, 1]); R_R[2, 2] := evalf(R[2, 2]); R_R[2, 3] := evalf(R[2, 3]); R_R[3] := evalf(R[3]); R_R[3, 1] := evalf(R[3, 1]); R_R[3, 2] := evalf(R[3, 2]); R_R[3, 3] := evalf(R[3, 3])

6.833764322

  

.5388679374

  

-1.033616758

  

4.934912437

  

-0.3143256678e-1

  

-4.822756148

  

4.934912437

  

-5.541440351

  

.2813149492

  

0.4011667570e-1

  

-0.3143256678e-1

  

.2813149492

  

-0.7207815682e-2

 (3)
> 

tstop := time()-tstart

24.165

 (4)

``


 

Download sumInt.mw

Your maple installation is OK (probably)...

tomleslie 13876
July 24 2018
1 1

OS = 64-bit Win7

Your worksheet executes with no problems in

Maple 18.02
Maple 2015.2
Maple 2016.2
Maple 2017.3

But Maple 2018.1 "hangs" ( it might finish eventually? ). I have no idea why. I'll try a few experiements and let you know

I have to admit that I was sceptical...

tomleslie 13876
July 23 2018
0 1

On 64-bit Win7, Maple 2016 returns

> 

restart;
kernelopts(version);
z[x]:=z;
z[x]:=z(x);

`Maple 2016.2, X86 64 WINDOWS, Jan 13 2017, Build ID 1194701`

  

z

  

(table( [( x ) = ] ))(x)

 (1)
> 

 

Download op2016.mw

On same OS, Maple 2018 returns

> 

restart;
kernelopts(version);
z[x]:=z;
z[x]:=z(x);

`Maple 2018.1, X86 64 WINDOWS, Jun 8 2018, Build ID 1321769`

  

z

  

Error, recursive assignment

  
> 

 

Download op2018.mw

I wish I could show the output for the same code in Maple 2017.1, but I get an mserver.exe crash every time I try to execute it!

The code itself doesn't make much sense, so I wouldn't have expected anything very sensible from any of Maple2016.2, Maple 2017.1 or Maple 2018.1 - but the outright crach in Maple 2017.1 did surprise me!!

Congratulations?! This is the shortest code sequence I have ever seen which causes Maple to reproducibly crash

I agree with vv...

tomleslie 13876
July 23 2018
0 0

in that technically in the construction in your second example

object:=Object||n

the name 'Object' will be treated as a 'gllobal', whilst 'n' is local. In general trems the attributes of the valeu returned through use of the || construct reflect those of the first argument, so I'm guessing that anything constructed by this means will actually return a global - which isn't what you want.

Best workaround I can come up with for this isuue is

restart:

f := proc()
              local Object1, Object2, Object3, Choice, A:
              Choice := proc(a)
                                        local l, n, object:
                                        l := [seq(args, k=1..nargs)]:
                                        n := ListTools:-Search(true, l):
                                       object := [Object1, Object2, Object3][n]:
                               end proc:

             Object1 := {x1, y1, z1}:
             Object2 := {x2, y2}:
             Object3 := {x3}:

          # here I give A some values, in the true application A is the sequence
          # returned by a Maplet through a Shutdown command
          # More precisely A correspond to “N” buttons of the same group … then
          # only one is “true”.

             A := false, true, false:
             Choice(A);
     end proc;

  f();

 

Some thoughts...

tomleslie 13876
July 23 2018
1 2

You have three equations in the four unknowns p,v,w,t.

It is possible to solve/plot this system - ie obtain explicit expressions for p(t), v(t) and w(t) and plot each with respect to 't', with a few caveats

  1. There isn't a unique "solution", there are five possible solutions.
  2. For the first solution set, v(t) can only be evaluated to reals over the range t=-1..1. Outside this range, it is complex: p(t) and w(t) appear to be real everywhere
  3. For the second solution set, w(t) can only be evaluated to reals over the range t=-1..1. Outside this range, it is complex: p(t) and v(t) appear to be real everywhere
  4. For the third solution set, v(t) and w(t) are complex ( everywhere? ): p(t) appears to be real everywhere
  5. For the fourth solution set v(t) and w(t) are complex ( everywhere ?): p(t) appears to be real everywhere
  6. For the fifth solution set 'v' is arbitrary,  w=f(v), so is also arbitrary, and p=f(v,t) so is also arbitrary

Check the attached. ( Normally I would display the contents here, but having used Array(plot...), this never seems to render here, so you will have to download to examine)

Download plotProb.mw

Three issues...

tomleslie 13876
July 21 2018
0 1
  1. It is not clear whether you actually want to write your own version of Dijkstra's algorithm (eg for educational purposes), or whether using MAple's in-built version is acceptable. I'm assuming the latter for the moment
  2. It is not clear whether your test graph is directed or undirected. The attached does both: I have assumed that for the case of a directed graph, the direction is from lower to the higher vertex number. If this is incorrect it would be pretty trivial to change the appropriate adjacency matrix.
  3. I think mmcdara may have misintepreted the tuple [14, 16, 1] which (s)he is interpreting as an 'edge' for vertex 14 to vertex 16, with weight 1. I interpret this according to the comment in your worksheet, that your graph has 14 vertexes, and 16 edges, with 'source node'=1

The attached shows the possibilities

> 

  restart;
  interface(rtablesize=15):
  with(GraphTheory):
#
# OP's graph shown below for reference. Not obvious
# whether this is intended to be a directed or an
# undirected graph, so work it out both ways
#
> 

#
# Define the adjacency matrix for the graph on the
# basis the all edges are directed, with the implied
# direction being from lower->higher vertex number.
#
# The adjacency matrix for such a graph will be upper
# triangular. For the case of an undirected graph
# (where the adjacency will be symmetric), see later
#
  A:= Matrix( 14,
              14,
              [ [ 0, 2, 0, 0, 0, 0, 3, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 7, 0, 0, 0, 0, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 3, 0, 0, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 6, 0, 0, 0, 3 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 4, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 3, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 7, 0 ,0, 0, 0, 3],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 3 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,3, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 2, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0, 5, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0, 0, 3],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0, 0, 0]
              ]
           ):
  G1:= Graph(A);
#
# Run Dijkstras Algorithm to give the cost and vertex list
# of every arc which starts at vertex 1. Notie that the
# output of this command i a list of lists, where each sublist
# comprises two elements
#
# 1  the vertex sequence for the desired arc: note thst the
#    final vertex in these sublists goes from 1 to 14
#    sequentially, so all possible end vertexes are covered
# 2  the cost (ie sum of weights) for the arc
#
  DijkstrasAlgorithm(G1, 1);
#
# Do the same thing starting from vertex 6. Note that because
# the graph is directed as described above, no vertexes
# numbered lower than 6 are attainable, so show with no vertex
# sequence and infinite cost
#
  DijkstrasAlgorithm(G1, 6);
#
# Specifying both the start and end-point is done with the
# following, where the second will fail because of the
# directed nature of the graph
#
  DijkstrasAlgorithm(G1, 6, 10);
  DijkstrasAlgorithm(G1, 10,6);

GRAPHLN(directed, weighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14], Array(%id = 18446744074331554990), `GRAPHLN/table/1`, Matrix(%id = 18446744074331554150))

  

[[[1], 0], [[1, 2], 2], [[1, 2, 3], 9], [[1, 2, 3, 4], 10], [[1, 2, 3, 4, 5], 13], [[1, 2, 3, 4, 5, 6], 19], [[1, 7], 3], [[1, 7, 8], 6], [[1, 7, 8, 9], 13], [[1, 2, 3, 4, 5, 10], 16], [[1, 2, 3, 4, 5, 10, 11], 19], [[1, 2, 3, 4, 5, 10, 11, 12], 21], [[1, 2, 3, 4, 5, 10, 11, 12, 13], 26], [[1, 7, 8, 14], 9]]

  

[[[], infinity], [[], infinity], [[], infinity], [[], infinity], [[], infinity], [[6], 0], [[6, 7], 4], [[6, 7, 8], 7], [[6, 7, 8, 9], 14], [[6, 7, 8, 9, 10], 17], [[6, 7, 8, 9, 10, 11], 20], [[6, 7, 8, 9, 10, 11, 12], 22], [[6, 7, 8, 9, 10, 11, 12, 13], 27], [[6, 7, 8, 14], 10]]

  

[[6, 7, 8, 9, 10], 17]

  

[[], infinity]

 (1)
> 

#
# Adjust the adjeacency matrix so that it is symmetric, and
# hence defines a weighted, undirected graph
#
  B:= Matrix( 14,
              14,
              [ [ 0, 2, 0, 0, 0, 0, 3, 0, 0, 0 ,0, 0, 0, 0],
                [ 2, 0, 7, 0, 0, 0, 0, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 7, 0, 1, 0, 0, 0, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 1, 0, 3, 0, 0, 0, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 0, 3, 0, 6, 0, 0, 0, 3 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 6, 0, 4, 0, 0, 0 ,0, 0, 0, 0],
                [ 3, 0, 0, 0, 0, 4, 0, 3, 0, 0 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 3, 0, 7, 0 ,0, 0, 0, 3],
                [ 0, 0, 0, 0, 0, 0, 0, 7, 0, 3 ,0, 0, 0, 0],
                [ 0, 0, 0, 0, 3, 0, 0, 0, 3, 0 ,3, 0, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 3 ,0, 2, 0, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,2, 0, 5, 0],
                [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 5, 0, 3],
                [ 0, 0, 0, 0, 0, 0, 0, 3, 0, 0 ,0, 0, 3, 0]
              ]
           ):\
  G2:=Graph(B);
#
# Same four case as tested before, although since the graph
# is now undirected, all arcs are now possible (and 6->10
# has the same arc and cost as 10->6
#
# All arcs+costs from vertex 1
#
  DijkstrasAlgorithm(G2, 1);
#
# All arcs+costs from vertex 6
#
  DijkstrasAlgorithm(G2, 6);
#
# arc+cost from vertex 6 to vertex 10
#
  DijkstrasAlgorithm(G2, 6, 10);
#
# arc+cost from vertex 10 to vertex 6. Should be same cost
# as above and 'reversed' vertex list - which it is
#
  DijkstrasAlgorithm(G2, 10, 6);

GRAPHLN(undirected, weighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14], Array(%id = 18446744074368532838), `GRAPHLN/table/2`, Matrix(%id = 18446744074368540294))

  

[[[1], 0], [[1, 2], 2], [[1, 2, 3], 9], [[1, 2, 3, 4], 10], [[1, 7, 6, 5], 13], [[1, 7, 6], 7], [[1, 7], 3], [[1, 7, 8], 6], [[1, 7, 8, 9], 13], [[1, 7, 6, 5, 10], 16], [[1, 7, 6, 5, 10, 11], 19], [[1, 7, 8, 14, 13, 12], 17], [[1, 7, 8, 14, 13], 12], [[1, 7, 8, 14], 9]]

  

[[[6, 7, 1], 7], [[6, 7, 1, 2], 9], [[6, 5, 4, 3], 10], [[6, 5, 4], 9], [[6, 5], 6], [[6], 0], [[6, 7], 4], [[6, 7, 8], 7], [[6, 5, 10, 9], 12], [[6, 5, 10], 9], [[6, 5, 10, 11], 12], [[6, 5, 10, 11, 12], 14], [[6, 7, 8, 14, 13], 13], [[6, 7, 8, 14], 10]]

  

[[6, 5, 10], 9]

  

[[10, 5, 6], 9]

 (2)
> 

 

> 

 


 

Download Dijkstra.mw

Interesting and not very reproducible!!...

tomleslie 13876
July 20 2018
0 0

I think Acer's comments are apposite - and should be considered carefully.

  1. Other than the actual timing procedure GTS2Timer2(), didn't examine your code at all, so I have no idea where time is being consumed, whether code will always terminate, whatever.
  2. Rather obviously, any results are going to be machine-dependent, so you (probably) won't be able to reproduce the output shown in the attached
  3. A little experimentation has shown me that results are not even reproducible on my own machine! About the only conclusion I could come to, is that it is always(?) the same three execution groups which cause problems. However the reported times, vary a lot
  4. I modified your try/catch statement to 'catch' the 'Time Expired' error, and just print a message when, this error occurs
  5. The most interesting case is probably the final one in the attached, where the 'Time Expired' error occurs on iterations 1, 3, and 5. Note however that for this case, the time taken on iterations 2 and 4 is actually longer than that taken on iteration 1, which did produce the error
  6. The situation described in (5) above is probably(?) an example of the limitations which Acer has described. For example, if (somewhere) a recursive function call is being used, then I doubt that the 'time' will be polled between recursions, but only when the recursion terminates (assuming it does!!). When the recursion terminates, and an answer has been obtained, there wouldn't really be any point in returning the 'Time Expired' error, rather the answer obtained. And if the recursion never termnates, then you could be waiting for a very loonnngggg time, which would be unaffected by the timelimit() setting

So with all of the above caveats, the best I could come up with, was the attached

Set up the model

 

``

> 

F := [k[a1]*C[T]*(R-x[1]-x[2])-k[d1]*x[1], k[a2]*C[T]*(R-x[1]-x[2])-k[d2]*x[2]]; H := 1000*(x[1]+x[2])

[k[a1]*C[T]*(R-x[1]-x[2])-k[d1]*x[1], k[a2]*C[T]*(R-x[1]-x[2])-k[d2]*x[2]]

 

1000*x[1]+1000*x[2]

(1.1)

Set up the Programs

 

 

> 

with(LinearAlgebra); with(Groebner); lieDer := proc (H, F) local N, V, vars; N := nops(F); vars := [seq(x[t], t = 1 .. N)]; V := map(proc (a, b) options operator, arrow; diff(b, a) end proc, vars, H); DotProduct(Vector(F), Vector(V), conjugate = false) end proc; listLieDer := proc (H, F, N) local L, i, tmp; L := [y[0]-H]; tmp := H; for i to N do tmp := lieDer(tmp, F); L := [op(L), y[i]-tmp] end do end proc; TypeTools:-AddType(hname, proc (n) local r; r := n; if not n::name then return false end if; while not r::symbol do r := op(0, r) end do; evalb(substring(r, -1 .. -1) = ':-h') end proc)

GTS := proc (H, F, Na, Nd) local Ha, Hd, hatsubs, Hdhat, Hahat, Sols; Ha := listLieDer(H, F, Na); Ha := subs([x[1] = 0, x[2] = 0, y = ya], Ha); Hd := listLieDer(H, subs(C[T] = 0, F), Nd); Hd := subs(y = yd, Hd); hatsubs := [k = kh, C = Ch, R = Rh, x = xh]; Hdhat := subs(hatsubs, Hd); Hahat := subs(hatsubs, Ha); Sols := solve([op(`~`[`-`](Ha, Hahat)), op(`~`[`-`](Hd, Hdhat)), x[1] <> 0, Ch[T] <> 0, C[T] <> 0]) end proc; GTS2 := proc (H, F, Na, Nd) local Sols, Sols2, Sols3; Sols := [GTS(H, F, Na, Nd)]; numelems(Sols); Sols2 := remove(proc (s) options operator, arrow; `in`(true, map(proc (e) options operator, arrow; `assuming`([is(rhs(e) <= 0)], [positive]) end proc, s)) end proc, Sols); numelems(Sols2); Sols3 := select(proc (s) options operator, arrow; andmap(hastype, remove(evalb, s), hname) end proc, Sols2); numelems(Sols3); Sols3 end proc

GTS2timer2 := proc (H, F, Na, Nd) local st; st := time(); try timelimit(20, GTS2(H, F, Na, Nd)) catch "time expired": printf("Time Expired - I'm baling out on iteration i=%d\n", Nd) finally time()-st end try end proc

proc (H, F, Na, Nd) local Sols, Sols2, Sols3; Sols := [GTS(H, F, Na, Nd)]; numelems(Sols); Sols2 := remove(proc (s) options operator, arrow; `in`(true, map(proc (e) options operator, arrow; `assuming`([is(rhs(e) <= 0)], [positive]) end proc, s)) end proc, Sols); numelems(Sols2); Sols3 := select(proc (s) options operator, arrow; andmap(hastype, remove(evalb, s), hname) end proc, Sols2); numelems(Sols3); Sols3 end proc

  

proc (H, F, Na, Nd) local st; st := time(); try timelimit(20, GTS2(H, F, Na, Nd)) catch "time expired": printf("Time Expired - I'm baling out on iteration i=%d ", Nd) finally time()-st end try end proc

 (2.1)

Can we Graph the time it takes as Na Nd varies

 
> 

Na := 1; Nd := 5; v1 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

1

 

5

 

[.702, .718, .639, .687, 4.368]

(3.1)
> 

Na := 2; Nd := 5; v2 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

2

 

5

 

[3.386, 6.162, 5.506, 5.944, 11.123]

(3.2)
> 

Na := 3; Nd := 5; v3 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

3

 

5

 

Time Expired - I'm baling out on iteration i=5

 

[3.105, 8.127, 5.772, 5.554, 46.847]

(3.3)
> 

Na := 4; Nd := 5; v4 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

4

 

5

 

[7.301, 23.572, 5.491, 5.678, 5.975]

(3.4)
> 

Na := 5; Nd := 5; v5 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

5

 

5

 

[9.469, 6.303, 6.989, 7.441, 7.270]

(3.5)
> 

Na := 6; Nd := 5; v6 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

6

 

5

 

[16.115, 8.003, 10.420, 10.983, 14.430]

(3.6)
> 

Na := 7; Nd := 5; v7 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

7

 

5

 

Time Expired - I'm baling out on iteration i=1

 

[22.807, 16.178, 20.498, 20.452, 16.879]

(3.7)
> 

Na := 8; Nd := 5; v1 := [seq(GTS2timer2(H, F, Na, i), i = 1 .. Nd)]

8

 

5

 

Time Expired - I'm baling out on iteration i=1

Time Expired - I'm baling out on iteration i=3

Time Expired - I'm baling out on iteration i=5

 

[22.932, 35.552, 39.250, 34.258, 41.746]

(3.8)
> 

with(plots); M := Matrix([v1, v2, v3, v4, v5, v6, v7]); matrixplot(M); matrixplot(`~`[log10](M)); m2 := Matrix([[1, 2, 3], [4, 5, 6]])

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

  

Matrix(7, 5, {(1, 1) = 22.932, (1, 2) = 35.552, (1, 3) = 39.250, (1, 4) = 34.258, (1, 5) = 41.746, (2, 1) = 3.386, (2, 2) = 6.162, (2, 3) = 5.506, (2, 4) = 5.944, (2, 5) = 11.123, (3, 1) = 3.105, (3, 2) = 8.127, (3, 3) = 5.772, (3, 4) = 5.554, (3, 5) = 46.847, (4, 1) = 7.301, (4, 2) = 23.572, (4, 3) = 5.491, (4, 4) = 5.678, (4, 5) = 5.975, (5, 1) = 9.469, (5, 2) = 6.303, (5, 3) = 6.989, (5, 4) = 7.441, (5, 5) = 7.270, (6, 1) = 16.115, (6, 2) = 8.003, (6, 3) = 10.420, (6, 4) = 10.983, (6, 5) = 14.430, (7, 1) = 22.807, (7, 2) = 16.178, (7, 3) = 20.498, (7, 4) = 20.452, (7, 5) = 16.879})

 

 

 

Matrix(%id = 18446744074608213398)

(3.9)

 

``

``

``

Download timeProb.mw

Not sure but I don't think so...

tomleslie 13876
July 19 2018
0 0

Best I can think of is to set up aliases, as in
 

> 

restart;
with(VectorCalculus):
SetCoordinates(cartesian[x, y, z]);
alias( vx=vx(x,y,z),
       vy=vy(x,y,z),
       vz=vz(x,y,z)
     ):
v := VectorField([vx, vy, vz]);

Jacobian(v);

cartesian[x, y, z]

  

vx, vy, vz

  

Vector(3, {(1) = vx(x, y, z), (2) = vy(x, y, z), (3) = vz(x, y, z)})

  

Matrix(%id = 18446744074372878198)

 (1)
> 

 


 

Download VCalias.mw

In the explict code...

tomleslie 13876
July 19 2018
1 1

you are using two differnt 'timer' functions. One of these is called GTS2timer() and the other is called GTS2timer2(). However, only the latter is defined (as far as I can tell). So when you call the former everything crashes

Changing all the 'timer' calls to GTS2timer2() (my best guess) doesn't help much, because every(?) test relies on executing the procedure GTS2(), which fails every time because it relies on the undefined type name 'hname'

I have no idea what the type name 'hname' should refer to: so I can't even guess a fix for this. Until code executes, not much one can do about debugging timing

Well I think I understand the problem...

tomleslie 13876
July 19 2018
0 0

because something very similar was asked, and answered before at

https://www.mapleprimes.com/questions/225167-Comparison-Of-Coefficient-Of-Exponential-Function

However the OP seems reluctant to read  advice previously given.

The attached will work, but is seriously "quick and dirty".

> 

restart
> 

Eq1 := -alpha[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^2+alpha[1]*beta[1]^2*exp(-t*omega[1]+x*alpha[1]+y*beta[1])+alpha[2]*beta[2]^2*exp(-t*omega[2]+x*alpha[2]+y*beta[2])+alpha[1]*omega[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*a[1, 2]+alpha[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*omega[2]*a[1, 2]+omega[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*a[1, 2]+exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*alpha[2]^3+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[1]^3+exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*alpha[1]^3+alpha[2]*omega[2]*exp(-t*omega[2]+x*alpha[2]+y*beta[2])+alpha[1]*omega[1]*exp(-t*omega[1]+x*alpha[1]+y*beta[1])+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^3*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^3-beta[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*beta[2]^2-alpha[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]+3*alpha[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^2*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*omega[2]*a[1, 2]+alpha[1]*beta[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*a[1, 2]+alpha[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]^2*a[1, 2]+beta[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*beta[2]^2*a[1, 2]+3*alpha[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*a[1, 2]+exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*alpha[2]*beta[2]^2+exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*alpha[1]*beta[1]^2+exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*alpha[2]*omega[2]+exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*alpha[1]*omega[1]+2*exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*alpha[2]*beta[1]*beta[2]+2*exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*alpha[1]*alpha[2]^2+2*exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*alpha[1]^2*alpha[2]+alpha[2]^3*exp(-t*omega[2]+x*alpha[2]+y*beta[2])+alpha[1]^3*exp(-t*omega[1]+x*alpha[1]+y*beta[1])-omega[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]-alpha[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*omega[2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*omega[2]+2*exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*alpha[1]*beta[1]*beta[2]+2*alpha[1]*beta[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*a[1, 2]+alpha[1]^3*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*a[1, 2]+alpha[1]*omega[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])-alpha[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]^2+alpha[1]*beta[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])+2*beta[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*beta[2]*a[1, 2]; Eq2 := -exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[1]^2*beta[2]-exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^2*beta[1]+exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*beta[2]^3+exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*beta[1]^3+omega[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*omega[2]*a[1, 2]+alpha[1]^2*beta[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*a[1, 2]+exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*beta[2]*omega[2]+exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*beta[1]*omega[1]+beta[1]*omega[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*a[1, 2]+beta[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*omega[2]*a[1, 2]+beta[1]^3*exp(-t*omega[1]+x*alpha[1]+y*beta[1])-exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[1]^2*beta[2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^2*beta[2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[1]*omega[1]-exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[1]*beta[2]^2+2*alpha[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*alpha[2]*a[1, 2]+beta[2]^3*exp(-t*omega[2]+x*alpha[2]+y*beta[2])+3*beta[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*a[1, 2]+3*beta[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]^2*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[1]^2*beta[1]+beta[1]*omega[1]*exp(-t*omega[1]+x*alpha[1]+y*beta[1])+2*exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*alpha[1]*alpha[2]*beta[1]+2*exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*alpha[1]*alpha[2]*beta[2]+beta[2]*omega[2]*exp(-t*omega[2]+x*alpha[2]+y*beta[2])+2*exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*beta[1]*beta[2]^2+2*exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*beta[1]^2*beta[2]+exp(-2*t*omega[1]-t*omega[2]+2*x*alpha[1]+x*alpha[2]+2*y*beta[1]+y*beta[2])*a[1, 2]*alpha[2]^2*beta[2]+exp(-t*omega[1]-2*t*omega[2]+x*alpha[1]+2*x*alpha[2]+y*beta[1]+2*y*beta[2])*a[1, 2]*alpha[1]^2*beta[1]+alpha[1]^2*beta[1]*exp(-t*omega[1]+x*alpha[1]+y*beta[1])+alpha[2]^2*beta[2]*exp(-t*omega[2]+x*alpha[2]+y*beta[2])-exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[1]*omega[2]+2*alpha[1]*beta[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*omega[2]-exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*omega[1]+beta[1]^3*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]^3*a[1, 2]+alpha[1]^2*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]*a[1, 2]+beta[1]*exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^2*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*alpha[2]^2*beta[2]*a[1, 2]+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[1]^3+exp(-t*omega[1]-t*omega[2]+x*alpha[1]+x*alpha[2]+y*beta[1]+y*beta[2])*beta[2]^3
> 

s1:=-t*omega[1]+x*alpha[1]+y*beta[1]=eta[1]:
s2:=-t*omega[2]+x*alpha[2]+y*beta[2]=eta[2]:
s1a:=exp(eta[1])=p:
s1b:=exp(eta[2])=q:
getCoeffs:= proc( eq )
               #
               # Since the problem has been so badly posed
               # the first thing we have to do is perform
               # a sequence of substitutions.
               #
               # The first step is to substitute eta[1] and
               # and eta[2] for some (lengthy) expressions
               #
               # The next step is to substitute dummy
               # variables for exp(eta[1])^power and
               # exp(eta[2])^power, so that one produces
               # a simple multivariate polynomial
               #
               # Then extract the coefficients of the simple
               # multivariate polynomial
               #
                  local myCoeffs:= [ coeffs
                                     ( collect
                                       ( algsubs
                                         ( s1b,
                                           expand
                                           ( algsubs
                                             ( s1a,
                                               expand
                                               ( algsubs
                                                 ( s1,
                                                   algsubs
                                                   ( s2,
                                                     eq
                                                   )
                                                 )
                                               )
                                             )
                                           )
                                         ),
                                         [p, q].
                                         'distributed'
                                       ),
                                       [p, q],
                                       'tt'
                                     )
                                   ];
               #
               # Now 'reverse' all of the substitutions which
               # were performed above - Wouldn't life be so
               # much easier if the starting equations were
               # presented in an appropriate form in the first
               # place?
                  return algsubs~
                         ( isolate
                           ( s2,
                             eta[2]
                           ),
                           algsubs~
                           ( isolate
                             ( s1,
                               eta[1]
                             ),
                             algsubs~
                             ( isolate
                               ( s1b,
                                 q
                               ),
                               algsubs~
                               ( isolate
                                 ( s1a,
                                   p
                                 ),
                                 [ seq
                                   ( [ myCoeffs[j],
                                       tt[j]
                                     ],
                                     j=1..numelems([tt])
                                   )
                                 ]
                               )
                             )
                           )
                         );
            end proc:
#
# Get the 'coefficients' of the exponential terms
# for each of the supplied equations.
#
# Results are obtained as a list of lists. Each
# sublist consist of two terms - the second is the
# exponential^power term, and the first is its coefficient
#
  getCoeffs(Eq1);
  getCoeffs(Eq2);

[[2*a[1, 2]*alpha[1]*alpha[2]^2+a[1, 2]*alpha[2]^3+2*a[1, 2]*alpha[2]*beta[1]*beta[2]+a[1, 2]*alpha[2]*beta[2]^2+a[1, 2]*alpha[2]*omega[2], (exp(-t*omega[1]+x*alpha[1]+y*beta[1]))^2*exp(-t*omega[2]+x*alpha[2]+y*beta[2])], [a[1, 2]*alpha[1]^3+2*a[1, 2]*alpha[1]^2*alpha[2]+a[1, 2]*alpha[1]*beta[1]^2+2*a[1, 2]*alpha[1]*beta[1]*beta[2]+a[1, 2]*alpha[1]*omega[1], exp(-t*omega[1]+x*alpha[1]+y*beta[1])*(exp(-t*omega[2]+x*alpha[2]+y*beta[2]))^2], [a[1, 2]*alpha[1]^3+3*a[1, 2]*alpha[1]^2*alpha[2]+3*a[1, 2]*alpha[1]*alpha[2]^2+a[1, 2]*alpha[1]*beta[1]^2+2*a[1, 2]*alpha[1]*beta[1]*beta[2]+a[1, 2]*alpha[1]*beta[2]^2+a[1, 2]*alpha[2]^3+a[1, 2]*alpha[2]*beta[1]^2+2*a[1, 2]*alpha[2]*beta[1]*beta[2]+a[1, 2]*alpha[2]*beta[2]^2+a[1, 2]*alpha[1]*omega[1]+a[1, 2]*alpha[1]*omega[2]+a[1, 2]*alpha[2]*omega[1]+a[1, 2]*alpha[2]*omega[2]+alpha[1]^3-alpha[1]^2*alpha[2]-alpha[1]*alpha[2]^2+alpha[1]*beta[1]^2-alpha[1]*beta[2]^2+alpha[2]^3-alpha[2]*beta[1]^2+alpha[2]*beta[2]^2+alpha[1]*omega[1]-alpha[1]*omega[2]-alpha[2]*omega[1]+alpha[2]*omega[2], exp(-t*omega[1]+x*alpha[1]+y*beta[1])*exp(-t*omega[2]+x*alpha[2]+y*beta[2])], [alpha[1]^3+alpha[1]*beta[1]^2+alpha[1]*omega[1], exp(-t*omega[1]+x*alpha[1]+y*beta[1])], [alpha[2]^3+alpha[2]*beta[2]^2+alpha[2]*omega[2], exp(-t*omega[2]+x*alpha[2]+y*beta[2])]]

  

[[2*a[1, 2]*alpha[1]*alpha[2]*beta[2]+a[1, 2]*alpha[2]^2*beta[2]+2*a[1, 2]*beta[1]*beta[2]^2+a[1, 2]*beta[2]^3+a[1, 2]*beta[2]*omega[2], (exp(-t*omega[1]+x*alpha[1]+y*beta[1]))^2*exp(-t*omega[2]+x*alpha[2]+y*beta[2])], [a[1, 2]*alpha[1]^2*beta[1]+2*a[1, 2]*alpha[1]*alpha[2]*beta[1]+a[1, 2]*beta[1]^3+2*a[1, 2]*beta[1]^2*beta[2]+a[1, 2]*beta[1]*omega[1], exp(-t*omega[1]+x*alpha[1]+y*beta[1])*(exp(-t*omega[2]+x*alpha[2]+y*beta[2]))^2], [a[1, 2]*alpha[1]^2*beta[1]+a[1, 2]*alpha[1]^2*beta[2]+2*a[1, 2]*alpha[1]*alpha[2]*beta[1]+2*a[1, 2]*alpha[1]*alpha[2]*beta[2]+a[1, 2]*alpha[2]^2*beta[1]+a[1, 2]*alpha[2]^2*beta[2]+a[1, 2]*beta[1]^3+3*a[1, 2]*beta[1]^2*beta[2]+3*a[1, 2]*beta[1]*beta[2]^2+a[1, 2]*beta[2]^3+a[1, 2]*beta[1]*omega[1]+a[1, 2]*beta[1]*omega[2]+a[1, 2]*beta[2]*omega[1]+a[1, 2]*beta[2]*omega[2]+alpha[1]^2*beta[1]-alpha[1]^2*beta[2]-alpha[2]^2*beta[1]+alpha[2]^2*beta[2]+beta[1]^3-beta[1]^2*beta[2]-beta[1]*beta[2]^2+beta[2]^3+beta[1]*omega[1]-beta[1]*omega[2]-beta[2]*omega[1]+beta[2]*omega[2], exp(-t*omega[1]+x*alpha[1]+y*beta[1])*exp(-t*omega[2]+x*alpha[2]+y*beta[2])], [alpha[1]^2*beta[1]+beta[1]^3+beta[1]*omega[1], exp(-t*omega[1]+x*alpha[1]+y*beta[1])], [alpha[2]^2*beta[2]+beta[2]^3+beta[2]*omega[2], exp(-t*omega[2]+x*alpha[2]+y*beta[2])]]

 (1)
> 

 


 

Download expPow.mw

 

Stopping criterion...

tomleslie 13876
July 19 2018
0 1

as the loop is executed, the loop index is incremented. Once the index value is maore than the end-value it stops.

If it helps think of the following sequence

  1. At start of loop, is index greater than terminal value?
  2. If it isn't, execute loop
  3. At end of loop iteration increment index and repeat from (1) above,
  4. So what is the loop index value when the loop terminates? Like every other programming language on the planet, it is always one more than the loop termination value

No chance...

tomleslie 13876
July 19 2018
0 1

Because re-executing your worksheet, re-runs the code which you wrote. In order for plot changes you make interactively to become 'permanent', these changes would have to edit/overwrite/modify the code which you wrote.

The importntt point to realise is that many plot options accept lists as values, so you can easily assign different values to different curves. Consider the attached


impPlot.mw





 

 

Your immediate problem...

tomleslie 13876
July 18 2018
1 3

is fixed in the final execution group of the attached. Notes

  1. I moved the plot p1 a little further to the south, so that it didn't bump into the text along the x-axis
  2. As a general rule, if you want to rotate, scale and translate plot structires, you should use the commands plottools:-rotate(), plottools:-scale() and plottools:-rotate(). So much easier than any other way, although I haven't bothered to apply this in the attached.
  3. Rendering of plots is different/better in a Maple worksheet than the way they are displayed here
> 

restart:
> 

with(Statistics):


GUMBEL COPULA

> 

N    := 1000:
v__1 := Sample(Uniform(0,1), N):
v__2 := Sample(Uniform(0,1), N):


theta := 5:  


K__c := w -> w*(1-log(w)/theta):

W := fsolve~(K__c(w) =~ v__2, w=0..1):


# (u__1, u__2) is a random vector with uniform marginal and a Gumbel copula dependence
u__1 := exp~(v__1^~(1/theta) *~ log~(W)):
u__2 := exp~((1 -~ v__1)^~(1/theta) *~ log~(W)):


# Marginals
f__1 := Quantile~(Normal(0, 1), u__1, numeric):
f__2 := Quantile~(LogNormal(1, 1/2), u__2, numeric):
> 


#-----------------------------------------------------
# f__1 histogram rotated by Pi radians
scale := max(f__2)-min(f__2):

h := Histogram(f__1):
hpg := plottools:-getdata(h):


# rotation and scaling of h
HPG := map(u -> u[3].<<1,0>|<0,-scale>> -~ convert([seq([0, min(f__2)], k=1..4)], Matrix), [hpg]):

p__1 := PLOT(seq(POLYGONS(HPG[k]), k=1..numelems(HPG)), COLOR(RGB, 0.8, 0.8, 0.8)):

#-----------------------------------------------------
# f__2 histogram rotated by Pi/2 radians

scale := max(f__1)-min(f__1):

h := Histogram(f__2):
hpg := plottools:-getdata(h):

# rotation and scaling of h

HPG := map(u -> (u[3].<<1,0>|<0,-scale>>)[.., [2, 1]] +~ convert([seq([min(f__1), 0], k=1..4)], Matrix), [hpg]):

p__2 := PLOT(seq(POLYGONS(HPG[k]), k=1..numelems(HPG)), COLOR(RGB, 0.8, 0.8, 0.8)):

#-----------------------------------------------------

plots:-display(
   p__1,
   p__2,
   plot(f__1, f__2, style=point, axes=boxed, title=cat("Gumbel copula (theta=", theta, ") with Normal / Lognormal marginals"))
)

 
> 

plots:-display(
   p__1,
   p__2,
   plot(f__1, f__2, style=point, title=cat("Gumbel copula (theta=", theta, ") with Normal / Lognormal marginals"))
)

 
> 

p1:= plottools:-translate( p__1, 0, -1):
p2:=plots:-display(p__2):
p3:=plot(f__1, f__2, style=point, title=cat("Gumbel copula (theta=", theta, ") with Normal / Lognormal marginals")):
pl1:=plottools:-line([-6,0], [3,0]):
pl2:=plottools:-line([0,0], [0,15]):
pl3:=plots:-textplot([seq( [0, j, convert(j,string)], j=5..15, 5)] , align=left):
pl4:=plots:-textplot([seq( [j, 0, convert(j,string)], j=-5..3)], align=below):
plots:-display([p1, p2, p3, pl1, pl2, pl3,pl4], axes=none);

 
> 

 

Download fixPlot.mw

More complete...

tomleslie 13876
July 18 2018
1 0

is the following


 

> 

  restart:
  with(plots):
  implicitplot
  ( [ seq( (x^2+y^2)^2 - 8*a*x*(x^2 - 3*y^2) + 18*a^2*(x^2+y^2)=27*a^4,
           a=1..3
         )
    ],
    x=-10..10,
    y=-10..10,
    numpoints=10000,
    color=[red,green,blue],
    title="Deltoids from implicitplot",
    titlefont=[times, bold, 18]
  );
  plot( [ seq( [ a*(2*cos(t) + cos(2*t)),
               a*(2*sin(t) - sin(2*t)),
               t=0..2*Pi
             ],
             a=1..3
           )
        ],
        color=[red, green, blue],
        title="Deltoids from parametric plot",
        titlefont=[times, bold, 18]
      )

 

 
> 

 


 

Download plotDelts.mw

exp/ln...

tomleslie 13876
July 18 2018
0 3

From the codegen:-cost() help page (my emphasis)

Note that the cost used for computing powers is as follows. For an integral power, repeated multiplication is assumed. For a general power it is assumed to be computed using exp and ln.

So that would make sqrt(x) be treated as exp( (1/2) * ln(x) ) - you may (or may not!) regard this as two functions+one multiplication, depending on how you think the (1/2) term should be dealt with.

As a general comment, I'm not sure I would expect to get "seriously optimized" Matlab code from Maple. My preferred approach would be to wrtie the Matlab code myself. If I need any Maple functionality form within Matlab, then (assuming you have the Maple-Matlab link set up, which should have happened on installation), you can actually drive Maple from Matlab as if it were a symbolic toolbox in the latter.

First 128 129 130 131 132 133 134 Last Page 130 of 207