At about 6-8 sec per objective evaluation, allowing NLPSolve to use its default methods which utilize derivatives, some sort of result returned after about 30 min. It took longer, hard coding the default quadrature method to _Gquad rather than _CCquad.

h := (c0, c1, c2, R, r) -> sqrt(R^2-r^2)*(c0+c1*r^2/R^2+c2*r^4/R^4):

K := (c0, c1, c2, R, r) -> ((diff(h(c0, c1, c2, R, r), r))*
   (1+(diff(h(c0, c1, c2, R, r), r))^2)+
   r*(diff(h(c0, c1, c2, R, r), r, r)))/
   (r*(1+(diff(h(c0, c1, c2, R, r), r))^2)^(3/2)):

f := (c0, c1, c2, R, r) -> r*sqrt(1+(diff(h(c0, c1, c2, R, r), r))^2)*
   K(c0, c1, c2, R, r)^2:

F:=proc(c0,c1,c2,R)
local res;
if not type(c0,numeric) then
  return 'procname'(args);
else
  res:=evalf(Int(f(c0,c1,c2,R,x),x=0..R,method=_CCquad));
  if not type(res,numeric) then
    res:=evalf(Int(f(c0,c1,c2,R,x),
                   x=0..R,epsilon=0.00001));
    if not type(res,numeric) then
      res:=1000.0;
    end if;
  end if;
end if;
res;
end proc:

E := proc (c0, c1, c2) F(c0, c1, c2, 3.91) end proc;

Optimization:-NLPSolve(E(c0, c1, c2),
       variables=[c0,c1,c2],
       initialpoint=[c0=.1035805,c1=1.001279,c2=-.561381]);

That returned,

[4.00000000007146994, [c0 = 1.00000053433322145,

                                 -5                              -5
    c1 = -0.328526474190508903 10  , c2 = 0.514721350878141067 10  ]]

Note that it might only be a local optimum.

And the smoothness of f might have been compromised in the above code, with that "fake" 1000.0 return value for individual failing objective value computations. Also, some of the individual quadrature results might only be accurate to the fallback tolerance of epsilon=0.00001 so that too may be affecting the result's accuracy (but not always).

Starting from initial point c0=1.0, c1=0.0, c2=0.0 took ten mins to get,

[4.00000000005466028, [c0 = 1.00000000001311972,

                                -8                              -6
    c1 = 0.165449448974536908 10  , c2 = 0.154581310067891980 10  ]]

Without raising Digits and adjusting the fallback epsilon (ie. smaller accuracy tolerance) it wasn't clear to me whether 1.0,0.0,0.0 was the actual extreme point. I didn't think about it analytically.

acer

Are you trying to plot the partial sums? If so, then then there are several way to do it. Here's one (not especially efficient) way.

# A
T:=Sum(1/(3*n^2+7*n+1)^n,n=1..N):
L:=evalf(subs(N=infinity,T));
plots:-display(`if`(L<infinity,plot(L,0..10),NULL),plots:-pointplot([seq([N,T],N=1..10)]));

# E
T:=Sum(2/(2*n+1),n=1..N):
L:=evalf(subs(N=infinity,T));
plots:-display(`if`(L<infinity,plot(L,0..10),NULL),plots:-pointplot([seq([N,T],N=1..10)]));

# G
T:=Sum(1/n^4,n=1..N):
L:=evalf(subs(N=infinity,T));
plots:-display(`if`(L<infinity,plot(L,0..10),NULL),plots:-pointplot([seq([N,T],N=1..10)]));

And so on.

acer

Are you saying that you've exported from Maple to LaTeX source, and now you need to tell your latex engine where to find the maple2e.sty file?

On Windows, it is located in the ETC folder, under the Maple installation. For example,

C:\Program Files\Maple 13\ETC\maple2e.sty

acer

It's not completely clear what you have in mind generally, as your only example is very simple.

Using the techniques laid out here, you can get an effect like this,

> a := 1:
> b := 2:

> c := a+b = f(a+b);

                              c := 1 + 2 = 3

Using the undocumented (and undebuggable) kernelopts(Evaluator) one could add a little more polish to that, I think. (It's a useful little gem that I only ever saw used internally to the "Warnings package" on the Applications Centre.) The same output as you showed could likely be produced just from c:=a+b as input.

Here's an Evaluator that uses an uneval trick.

> restart:

> kernelopts(Evaluator=proc(x::uneval) x=eval(x); end proc):

> a:=1:
> b:=2:

> a+b;
                                  a + b = 3

So, it might work by combining the two methodologies together. Or maybe ditch the Evaluator suggestion, and just refine the literal expressions code from the above link.

acer

Congratulations on the acute observation, as well as your circumspection.

You are quite right in that `op` can be used in this way, to get the dimensions of a Matrix.

And it's sensible to ask whether it can be relied upon to continue to work. As a general rule, I would usually advise using the command (such as LinearAlgebra:-Dimension), over picking apart the structure.

But let's observe two things. First, the `op` command has worked in this way for every release from Maple 6 through to Maple 13. And second, `op` is used in this way in many places within the Maple Library, including LinearAlgebra routines.

I would thus suggest that your choice might depend on the performance cost versus legibility. Will your code call this operation many times, so that the extra Library level function calls would be a measurable hit?

acer

Using names for the numeric values forestalls (numeric) evaluation. But it is (all) the arithmetic operations which the OP wanted to delay. So, how about rebinding those operators with inert versions. And then, just to get it to display nicely, create print extensions for those inert operators which themselves utilize names-for-numerics.

> restart:

> p:=module() export `+`, `^`, abs;
> option package;
>   `+`:=proc() &+(args); end proc:
>   `^`:=proc() &^(args); end proc:
>   abs:=proc() &abs(args); end proc:
> end module:

> `print/&+`:=proc(a,b) :-`+`(`if`(type(a,numeric),convert(a,name),a),
                             
> `if`(type(b,numeric),convert(b,name),b)); end proc:

> `print/&^`:=proc(a,b) :-`^`(`if`(type(a,numeric),convert(a,name),a),
                            
> `if`(type(b,numeric) and b<-1,convert(b,name),b)); end proc:

> `print/&abs`:=proc(a) :-`abs`(`if`(type(a,numeric),convert(a,name),a)); end proc:

> f:=proc(X)
>    subsindets(X,specfunc(numeric,
>       {`&^`,`&+`,`&abs`}),
>       t->`if`(op(0,t)=`&^`,:-`^`(op(t)),
>          `if`(op(0,t)=`&+`,:-`+`(op(t)),:-abs(op(t)))));
> end proc:

> with(p):

> expr:=12/(1-abs(2^2-3^2))+abs(-7)/abs(6-2^2)^2;
                                  12             | -7 |
                    expr := --------------- + ------------
                                   2    2             2  2
                            1 - | 2  - 3  |   | -6 + 2  |
 
> f(expr);
                               12              7
                         -------------- + -----------
                         1 - | 4 + -9 |             2
                                          | 6 + -4 |
 
> f(%);
                                  12         7
                              ---------- + ------
                              1 - | -5 |        2
                                           | 2 |
 
> f(%);
                                   12      7
                                 ------ + ----
                                 1 + -5     2
                                           2
 
> f(%);
                                   12
                                  ---- + 7/4
                                   -4
 
> f(%);
                                   -3 + 7/4
 
> f(%);
                                     -5/4

acer

If you know some keyword(s) that you remember appearing in the thread then you could try google itself (using `mapleprimes` & keywords). But I don't  advise spending time trying mapleprimes own search facility -- it is heavily broken and does not seem to have been recording new posts for many months now. (See here.)

Or you could manually search all the forums.

acer

Maybe you could use the Statistics package, with high working precision.

> restart:with(Statistics):

> X:=RandomVariable('Normal'(0,1)):

> evalf[50](2*CDF(X,1)-1);
             0.68268949213708589717046509126407584495582593345319
 
> evalf[50](Quantile(X,(%+1)/2));
             0.99999999999999999999999999999999999999999999999997

> evalf[50](Quantile(X,CDF(X,1)));
             0.99999999999999999999999999999999999999999999999996

Did I do that right?

acer

Use Array instead of array, with the first letter capitalized.

A:=Array(1..3):
for i from 1 to 3 do
  A[i]:=f(i);
end do:
A;
                           [f(1) f(2) f(3)]

B:=Array(1..3,i->f(i));
                         B:=[f(1) f(2) f(3)]

See the help-page ?Array

acer

Have you tried using `map` to apply `diff` to each entry of the Vector?

acer

> restart:
> H:=Array(1..2^16,datatype=float[8],storage=rectangular,order=C_order);
                              [ 1..65536 1-D Array   ]
                         H := [ Data Type: float[8]  ]
                              [ Storage: rectangular ]
                              [ Order: C_order       ]
 
> rtable_dims(H);
                                  1 .. 65536
 
> rtable_options(H);
 datatype = float[8], subtype = Array, storage = rectangular, order = C_order
 
> rtable_options(H,order,storage);
                             C_order, rectangular

This is one of a few such helpful routines which can be used together to copy an Array/Matrix/Vector.

> P:=Array(1..3,[11,13,17],datatype=float[8],storage=rectangular,order=C_order);
                             P := [11., 13., 17.]
 
> rtable(rtable_dims(P),rtable_elems(P),rtable_options(P));
                                [11., 13., 17.]

> interface(verboseproc=3):
> showstat(copy,2);
 
copy := proc(A)
local X, r, e;
       ...
   2     rtable(rtable_indfns(A),rtable_dims(A),A,rtable_options(A),readonly = false)
       ...
end proc

acer

Note that below I use an Array (capitalized).

> X:=Array(3): # originally with 3 elements only

> for i from 1 to 5 do
>  X(i):=i^2;
> end do:

> X; # using Maple 13, X can be "grown" to be of length 5
                        [1 4 9 16 25]

In much older versions of Maple, I would not have been able to grow X in size, and would have been limited to its original 3 entry spots. And I would have had to use the syntax X[i]:=... instead of X(i):=...

I could also have used a table structure, which is yet another type of mutable structure with an open-ended number of entries.

Note also that X was not created using X:=[...] as that would have created a list rather than an Array. See here for a bit more about comparison between Maple's data structures.

acer

I may be in the minority here, but I prefer using the Matrix() constructor with its listlist argument for the data. For example, I find that I only very rarely get the syntax wrong for something like this,

Matrix([[1,3],[-3,1]]);

You might also consider guessing what the common syntax mistakes are with using the angle brackets, and perhaps even convert/accept typical mistakes. For example, you could convert, <<1,3>,<3,-1>> and give partial or full marks.

acer

The Tolerances package uses evalr.

> R1:=INTERVAL(76.0 .. 84.0):
> R2:=INTERVAL(114.0 .. 126.0):
> Req:=R1*R2/(R1+R2):

> evalr(Req);
                     INTERVAL(41.25714285 .. 55.70526317)
 
> evalr(1/(expand(1/Req)));
                     INTERVAL(45.59999998 .. 50.40000003)

Now, that second result was just lucky. Terms simplified and repeated instances of R1 and R2 vanished before evalr got its hands on them.

Are you trying to bound the behaviour of the system, given bounds on the capabilities of the components? If so, then I can imagine that a Tolerances package that functions well could be useful.

If not, then are you trying to characterize behaviour given physical measurements of the resistances (with of course have an associated measurement error)? In that case, is the ScientificErrorAnalysis package of any use? Note that for that package, the error does not represent an interval, but rather a distribution.

> restart:

> R1:=ScientificErrorAnalysis:-Quantity(80.1,4.0)*Unit(ohm):
> R2:=ScientificErrorAnalysis:-Quantity(120.1,6.0)*Unit(ohm):

> Req:=combine(R1*R2/(R1+R2),errors);
           Req := Quantity(48.05199800 [Omega], 1.730531812 [Omega])

acer

map(trunc,a);

map(round,a);

trunc~(a);

round~(a);

acer