V:=Vector[row](8,(i)->i);

                        V := [1, 2, 3, 4, 5, 6, 7, 8]

with(ArrayTools):

m1:=Reshape(V,[2,4]);

                                   [1  3  5  7]
                             m1 := [          ]
                                   [2  4  6  8]

m2:=Alias(V,[2,4]);

                                   [1  3  5  7]
                             m2 := [          ]
                                   [2  4  6  8]

V[3]:=17:

m1[1,2]; # V and m1 do not share data

                                      3

m2[1,2]; # the difference, V and m2 share data

                                     17

whattype(m1), whattype(m2);

                               Matrix, Matrix

acer

Using the GlobalOptimization package (under 64bit Maple 14.01 on Windows 7),

restart:
with(GlobalOptimization):

GlobalSolve(
   (3*a/sqrt(1-u^2)+b/sqrt(1-t^2))/c,
   {a^2+2*b*c*u >= b^2+c^2, a*t+b*u <= c,
   b^2*(t^2-u^2)/(t^2-1)+c^2 <= 2*b*c*u},
   t=0..1, u=0..1, a=0..100, b=0..100, c=0..100);

   [4.00000000000000711, [a = 35.3564579217052, b = 35.3563778458042, 

     c = 50.0015242889233, t = 0.707107562133025, u = 0.707105960663738]]

That was for the default 'multistart' method. Unfortunately it's not clear how it chooses initial points. As the following userinfo messages indicates it appears to have already obtained a good result prior to the userinfo printouts. These messages seem to indicate that it has found the objective value very close to 4.0000 during some initial local search phase, before any global search phase has begun. Notice that the default feasibility tolerance used by GlobalSolve is 1.0e-4 which is really quite coarse a possible violation. But I was also able to get the same result from `GlobalSolve` using feasibility=1e-11, but no finer. The regular `Optimization:-Minimize` command can also find the same objective value if supplied with that same coarse feasibility tolerance (but doesn't get nearly as close to 4.00000000000000 for feasibilitytolerance=1e-5 unless the iterationlimit is supplied as high as 400).

restart:
with(GlobalOptimization):
infolevel[GlobalOptimization]:=100:
GlobalSolve(
   (3*a/sqrt(1-u^2)+b/sqrt(1-t^2))/c,
   {a^2+2*b*c*u >= b^2+c^2, a*t+b*u <= c,
   b^2*(t^2-u^2)/(t^2-1)+c^2 <= 2*b*c*u},
   t=0..1, u=0..1, a=0..100, b=0..100, c=0..100,
   method=singlestart);
GlobalSolve: calling NLP solver
SolveGeneral: calling global optimization solver
SolveGeneral: number of problem variables 5
SolveGeneral: number of nonlinear inequality constraints 3
SolveGeneral: number of nonlinear equality constraints 0
SolveGeneral: method singlestart
SolveGeneral: merit function evaluation limit 8000
SolveGeneral: non-improving merit function evaluation limit 1600
SolveGeneral: constraint penalty multiplier 100.0
SolveGeneral: target merit function value -0.10e11
SolveGeneral: local search target objective function value -0.10e11
SolveGeneral: local search feasibility tolerance 0.10e-5
SolveGeneral: local search optimality tolerance 0.10e-5
SolveGeneral: time limit in seconds 800
SolveGeneral: trying evalhf mode
Extcall: Lower bounds, nominal values, and upper bounds of decision variables:
Extcall: 1   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 2   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 3   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 4   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: 5   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: --- Initial function evaluation at nominal solution ---
Extcall: Objective function value 4.61880215351700585
Extcall: 1   5.000000e+001
Extcall: 2   5.000000e+001
Extcall: 3   5.000000e+001
Extcall: 4   5.000000e-001
Extcall: 5   5.000000e-001
Extcall: Constraint function values, 1, 0.
Extcall: Constraint function values, 2, 0.
Extcall: Constraint function values, 3, 0.
Extcall: The relative objective function improvement isless than .100000000000000002e-7 in 3 subsequent major iterations.
Extcall: The local search phase is terminated
Extcall: Number of function evaluations: 259
Extcall:  Current optimum estimate: 4.00000000000000533
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: 1   -9.094947e-013
Extcall: 2   0.000000e+000
Extcall: 3   4.547474e-013
Extcall: --- Automatic Global Optimization Procedure ---
Extcall: --- Global Adaptive Random Search ---
SolveGeneral: trying evalf mode
Extcall: Lower bounds, nominal values, and upper bounds of decision variables:
Extcall: 1   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 2   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 3   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 4   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: 5   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: --- Initial function evaluation at nominal solution ---
Extcall: Objective function value 4.61880215351700585
Extcall: 1   5.000000e+001
Extcall: 2   5.000000e+001
Extcall: 3   5.000000e+001
Extcall: 4   5.000000e-001
Extcall: 5   5.000000e-001
Extcall: Constraint function values, 1, 0.
Extcall: Constraint function values, 2, 0.
Extcall: Constraint function values, 3, 0.
Extcall: The relative objective function improvement isless than .100000000000000002e-7 in 3 subsequent major iterations.
Extcall: The local search phase is terminated
Extcall: Number of function evaluations: 259
Extcall:  Current optimum estimate: 4.00000000000000711
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: 1   -1.364242e-012
Extcall: 2   0.000000e+000
Extcall: 3   1.818989e-012
Extcall: --- Automatic Global Optimization Procedure ---
Extcall: --- Global Adaptive Random Search ---
Extcall: The relative objective function improvement isless than .100000000000000002e-7 in 3 subsequent major iterations.
Extcall: The local search phase is terminated
Extcall: Best solution found in GARS + LS search mode
Extcall: Number of function evaluations: 8398
Extcall:  Current optimum estimate: 4.00000000000000711
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: Constraint function values, 1, -.136424205265939236e-11
Extcall: Constraint function values, 2, 0.
Extcall: Constraint function values, 3, .181898940354585648e-11
Extcall: --- Summary of Optimization Results ---
Extcall: Total number of function evaluations: 8399
Extcall: Estimated optimum (including penalties) 4.00000000000000711
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: All stated constraints are met within tolerance .99999999999999995e-6
Extcall: Sum of constraint violations: .181898940354585648e-11
Extcall: Solution time (seconds), including user and system I/O operations, 0.
SolveGeneral: total number of function evaluations 8399
SolveGeneral: runtime in external solver 0.
SolveGeneral: maximum constraint infeasibility 0.181898940354585648e-11
SolveGeneral: cycling or stall detected in solver

       [4.00000000000000711, [a = 35.3564579217052, b = 35.3563778458042, 
                              c = 50.0015242889233, t = 0.707107562133025,
                              u = 0.707105960663738]]

acer

Does it help at all, to set u=t=sqrt(2)/2 ?  Just curious.

What about if you search on a=0.0..0.0001, b=0.0..0.0001, c=0.0..0.001 ?

acer

@Subjectivity Here's an attempt.

Download noiseweight.mw

What hapens if you "declare" the parameters of procedure `iang1`, using double-colon?

I mean, something like,

iang1 := proc (m1::float[8], vF1::float[8], omega::float[8], q::float[8]) ....

This would help the Compiler know that these were not going to be integer, with whatever other effects would ensure for that in the code (casts, etc).

acer

You want S.E.T to be approximately diagonal, right? (...by defn of Smith Form)

If so, then you'll need to increase the working precision (ie. `Digits`). But as the working precision gets higher, the rank of floating-point Matrix E may change. (The small entries of E can become significant, as the working precision increases.)

restart:
interface(rtablesize=100):

E := ImportMatrix(cat(kernelopts(homedir),
                      "/Downloads/model15_E.txt"),
                  delimiter=" "):

with(LinearAlgebra):

for i from 10 to 17 by 1 do
   Digits:=i;
   S,SNF,T := evalf(SmithForm(E,output=['U','S','V']));
   printf("   %ld      %3.3g     %.3e   %ld   %ld\n",
          Digits, Determinant(E), Norm(S.E.T - SNF),
          Rank(SNF), Rank(E));
end do:

for i from 10 to 17 do
   Digits:=i;
   S,SNF,T := evalf(SmithForm(E,output=['U','S','V'])):
   Z:=S.E.T:
   printf("Digits=%ld\n",i);
   printf("%.3f\n\n\n",map(fnormal,Z));
end do:

acer

You've defined `r` as an expression, but you use it as if it were a parameter name in your call to the plot3d command.

After assigning to your constants, the indeterminate names left in `Top` and `Bottom` are `phi` and `v`.

Did you instead mean something like,

plot3d([Top, Bottom], phi = s .. 20, v = 0 .. 2*Pi);

What coordinate system are you trying to plot in?

What are these two lines supposed to accomplish? (nb. `u` doesn't even appear anywhere else.)

u := u1+(u2-u1)*JacobiSN((1/2)*phi*sqrt(s*(u3-u1)), k)^2; 

r := 1/u;

acer

Try it with `algsubs` instead of `subs`.

algsubs(y2*y3 = 0, u(y2, y3));

                             2     2
                           y2  + y3 

acer

When you write of testing inequality, do you mean that you want it to return 'false' if at least one pair of entries violates the condition? (That's as opposed to all pairs having to violate it.)

If so then it can be faster to have a process which quits upon the first elementwise violation and which does not necessarily have to walk all entries. (I'm sure Joe Riel can beat whatever I come up with, soundly, for speed.)

I don't think that there are commands `orzip` and `andzip`.

Suppose that we wish to test whether A ≥ B, and by which we have taken the meaning to be that this is true only if A[i,j]≥B[i,j] for all i,j.

restart:
randomize(4567890):
with(LinearAlgebra):
A:=RandomMatrix(4,4,generator=10..100):
st:=time():
count:=0:
for i from 1 to 10000 do
   B:=RandomMatrix(4,4,generator=-80..20);
   if ormap(`<`,A-B,0) then
      # restriction was false, for at least one pair of entries
      count:=count+1;
  end if;
end do:
printf("\n     %a iterations produced false in %.3f seconds\n\n",count,time()-st);

     887 iterations produced false in 8.767 seconds

restart:
randomize(4567890):
with(LinearAlgebra):
A:=RandomMatrix(4,4,generator=10..100):
st:=time():
count:=0:
for i from 1 to 10000 do
   B:=RandomMatrix(4,4,generator=-80..20);
   try
      zip(proc(a,b) if a < b then error "found"; end if; end proc, A, B);
   catch "found":
      # restriction was false, for at least one pair of entries
      count:=count+1;
   end try;
end do:
printf("\n     %a iterations produced false in %.3f seconds\n\n",count,time()-st);

     887 iterations produced false in 0.671 seconds

restart:
randomize(4567890):
with(LinearAlgebra):
A:=RandomMatrix(4,4,generator=10..100):
st:=time():
count:=0:
for i from 1 to 10000 do
   B:=RandomMatrix(4,4,generator=-80..20);
   if ormap(evalb,`<`~(A,B)) then
      # restriction was false, for at least one pair of entrie
      count:=count+1;
  end if;
end do:
printf("\n     %a iterations produced false in %.3f seconds\n\n",count,time()-st);

     887 iterations produced false in 0.405 seconds

acer

Are you saying that you don't like the choice of colors in the default palette of Maple 16?

If so, then you can change the color palette using the plots:-setcolors command, eg.

plots:-setcolors("Classic"):

See also the setcolors help-page.

restart:
plot([x,x^2,x^3,x^4,x^5,x^6,x^7,x^8,x^9],x=0..Pi);
plots:-setcolors("Classic"):
plot([x,x^2,x^3,x^4,x^5,x^6,x^7,x^8,x^9],x=0..Pi);
plots:-setcolors("oldplots"):
plot([x,x^2,x^3,x^4,x^5,x^6,x^7,x^8,x^9],x=0..Pi);
plots:-setcolors("default"):
plot([x,x^2,x^3,x^4,x^5,x^6,x^7,x^8,x^9],x=0..Pi);

acer

How about a different approach, since you're just summing terms.

G := sqrt(3)/4 + Sum( sqrt(3)/12*(2/3)^(2*n-2), n=1..N );

                      /  N                         \
                      |-----                       |
                      | \                 (2 n - 2)|
           1  (1/2)   |  )   1   (1/2) /2\         |
           - 3      + | /    -- 3      |-|         |
           4          |----- 12        \3/         |
                      \n = 1                       /

value(G);

                                        (N + 1)
                2  (1/2)   27  (1/2) /4\       
                - 3      - -- 3      |-|       
                5          80        \9/       

limit(%,N=infinity);

                            2  (1/2)
                            - 3     
                            5       

value(subs(N=infinity,G));

                            2  (1/2)
                            - 3     
                            5       

limit(G,N=infinity);

                            2  (1/2)
                            - 3     
                            5       

sqrt(3)/4 + sum( sqrt(3)/12*(2/3)^(2*n-2), n=1..infinity );

                            2  (1/2)
                            - 3     
                            5       

acer

Your worksheet specifically does great many calls to int(...,numeric=false). If these are replaced with calls to evalf(Int(...)), and if all the `assume` stuff is removed, then I see a reduction from about 2min to about 25sec. The numeric X results and the plots seem the same.

If you thought that those integrals might succeed generically (symbolically), then they'd then be better placed outside the loop. The assumptions might serve a purpose, in this case. The same goes for attempting a generic symbolic `solve` instead of using fsolve in the loop. I didn't try any of this.

Another modication is for the fsolve calls, inside the loop. The idea is to use the previous iteration's fsolve solution as a starting point for the current iterations's fsolve call, on the supposition that the roots may not have moved much (your plots are mostly smooth, and perhaps expectedly so). If this call to fsolve (with initial point supplied) fails to converge then a second call to fsolve is made, without any initial point supplied. This produces the same final results, but in about 13sec.

Apart from removing the `assume` stuff, I also removed the assignments to {hb, hbt, etc) and their unassignments. It's enough to use 2-argument `eval`, to get the results into X.

The integrands may still contain variables whose values must be obtained by scoping. I mean the refences in the integrands to assignments made even outside of your procedure. If that is the case then it would mean that fast evalhf mode cannot be used for the numeric integration solvers. It's possible that you could get more savings if the integrands were all evalhf'able. I haven't the time to confirm this now, sorry.

qn_mp_v1_modif2.mw

acer

> a := binomial(37, x-105)*.85^(142-x)*0.15^(x-105):

> Optimization:-Maximize(a, x=105..142);

       [0.18320882465623625, [x = 110.194002601565]]

acer

Notice the differences between surd, ^, and root.

(You could also try the RealDomain package. And you might enjoy this list by Robert Israel.)

Download surdystuff.mw

acer

The first argument to StringTools:-Remove should be a predicate in Maple 12. That's not the same as true or false as a result.

a := "Just try to remove a letter using Maple 12's StringTools[Remove]":

StringTools:-Remove(t->evalb(t="e"), a);

         "Just try to rmov a lttr using Mapl 12's StringTools[Rmov]"

StringTools:-SubstituteAll(a,"e","");

         "Just try to rmov a lttr using Mapl 12's StringTools[Rmov]"

StringTools:-Remove("e", a); # Maple 16

         "Just try to rmov a lttr using Mapl 12's StringTools[Rmov]"

acer