Your un is created as a list. Don't use a list for this sort of task. Use a mutable data structure (for which element replacement is both possible and appropriate) such as a Vector, Array, or table instead. For example,

un := Vector(6):

See the section on data structures in the Maple Portal, or the Introductory Programming Guide, for details on those objects.

ps. That Portal page could explain "replacement" of a list element more fully (by describing subsop on a list and what that entails.)

acer

Don't use vector, matrix, and array. They are deprecated. Use Vector, Matrix and Array.

CurveFitting:-PolynomialInterpolation doesn't accept lowercase vector and matrix. The printing difference you mentioned is due to the fact that lowercase vector and matrix have last_name_eval.

Maple is a case-sensitive language.

ps. It is highly misleading to have the help-page vector(deprecated) be titled "Overview of Vectors". It's no wonder that new users find all this confusing. Fortunately the matrix(deprecated) and array(deprecated) help-pages haven't been treated to such a misguided mix-up with the capitalization in their titles.

acer

The purpose of that op() call is to extract the integrand of the inert Int call.

The general idea of the code seems to be to extract the integrand expression, get an alternate form, make a procedure that computes that, optimize that procedure, instantiate that at a dummy (r) to get a new integrand expression, recreate a new Int using the new integrand, and then finally to do numeric quadrature.

 F := Int(expr,t=a..b);
                                       b
                                      /
                                     |
                               F :=  |   expr dt
                                     |
                                    /
                                      a

> op(F)[1];
                                     expr

> IntegrationTools:-GetIntegrand(F); # another way
                                     expr

Why do you insert uneval quotes about commands, which prevent them from doing their intended job?

acer

It's not clear at what precision you intend the calculations to be done. That Generate call specifies 20 digits for the random values, but the environment variable Digits is not set above the default value (10).

If you would be satisfied with double precision then (leaving Digits alone) you could create Array A with the datatype=float[8] option. That could allow Statistics:-Mean (etc) to work with it using fast, memory efficient compiled routines.

For your loop of 10 million iterations I could imagine that quite a bit of Maple's time would be spent in creating and garbage collecting all those lists (t, and sorted t). They are small, but there are 20 million of them produced.

It might be useful if the ArrayTools package got a new routine to sort Vectors efficiently (in the memory sense, acting either inplace in the orginal or on a supplied container Vector for the result). More below of getting such efficient functionality yourself.

It's been mentioned before on this site that an improved Statistics:-Sample might allow one to re-use a float[8] Vector by populating it inplace. (That's how some good non-Maple random number generators work.) In the absence of that, you might consider the memory hit of generating all n*10^7 random values up front in a single float[8] Vector. Then you could simply act on clumps of 6 entries at a time. The code would keep track of which clump number it was working on and simply index into the large Vector appropriately. Such a large Vector would take about 500MB to allocate, and less than 10 seconds to create.

> with(Statistics):
> Dist := RandomVariable(Uniform(0,1)):
> st,ba:=time(),kernelopts(bytesalloc):
> T := Sample(Dist,6*10^7):
> time()-st,kernelopts(bytesalloc)-ba;
                               9.260, 480065524

This can also give an idea of how fast Maple could compute the mean of such a hardware float[8] Vector,

> st,ba:=time(),kernelopts(bytesalloc):
> Mean(T);
                                 0.4999735416

> time()-st,kernelopts(bytesalloc)-ba;
                                   0.148, 0

If 500MB allocation at once is too much then you could consider acting on a a fewer, repeated number of data Vectors -- say ten float[8] Vectors each with n*10^6 entries. That way only about 90MB would need to be allocated at once, with only six large Vectors that could be garbage-collected quickly. The idea is to play off total memory allocation against the cost of garbage collection of all those t's and sorted-t's.

If you feel adventurous, you might write a procedure which accepts a Vector and sorts it inplace. That could then be hit with Maple's Compiler:-Compile to provide a version which runs very fast and leanly. Since your n=6 then the complexity of the sorting algorithm wouldn't be relatively important. You could either create a single workspace float[8] Vector named sorted_t and repeatedly use ArrayTools:-Copy to get the current clump of the large T data Vector into it, or even better you could write the sorting procedure to accept full T as well as the "clump number" and then index into T appropriately in order to sort just the current clump of n entries inplace.

The procedure to be compiled could also do that max(...) work and populate float[8] Array A. In other words, write a single procedure to do it all, the sorting, the max&abs work, and the assignment into A. For fun, you could precompute one (or two) float[8] Vectors containing the i/n (or (i-1)/n) for quick re-use.

acer

This is called 2-argument eval, or evalat ("eval at", because it evaluates at some values). Notice that the instantiation can be obtained without having to actually assign to alpha. That can be useful because it leaves the original object as is, while also leaving name alpha unassigned and thus free immediately for continued symbolic use (no need to unassign it).

> eval(Result,alpha=Pi/2);
                                 [0    0    1]
                                 [           ]
                                 [1    0    0]
                                 [           ]
                                 [0    0    1]

I should mention that lowercase matrix and linalg are deprecated in favour of capitalized Matrix and LinearAlgebra.

acer

> s := series(-1/(x-1),x,3);
                                          2      3
                            s := 1 + x + x  + O(x )
 
> convert(s,polynom);
                                           2
                                  1 + x + x

acer

Look at the conditioning of your Matrix A.

> with(LinearAlgebra):
> S:=<<7476.45,303.62>|<303.62,26.19>>:
> I2:=IdentityMatrix(2):
> charPoly:=Determinant(x * I2 - S):

> EV:=solve(charPoly=0,x):
> ConditionNumber( EV[1]*I2 - S );
                                              11
                               0.7167270145 10
 
> Digits:=60:
> EV:=solve(charPoly=0,x):
> ConditionNumber( EV[1]*I2 - S ): Digits:=10: evalf(%%);
                                              61
                               0.1947557443 10
 
> Digits:=600:
> EV:=solve(charPoly=0,x):
> ConditionNumber( EV[1]*I2 - S ): Digits:=10: evalf(%%);
                                             602
                              0.1324339061 10

This is why singular values and singular vectors are used to compute the rank and the nullspace of floating-point Matrices. Numerically appropriate algorithms get used.

Trying to compute the RREF of a float Matrix is dubious. (Even a Matrix with two rows which are "obviously" simple scalar multiples of each other can "reduce" to something which appears to be the identity Matrix at the current working precision.) Also dubious is attempting to calculate the nullspace (kernel) using a method that is not appropriate for the floating-point domain.

> Digits:=10:
> ReducedRowEchelonForm( EV[1]*I2 - S );
                      [                              -14]
                      [1.    -0.299550326690223283 10   ]
                      [                                 ]
                      [0.                1.             ]
 
> infolevel[LinearAlgebra]:=1:
> NullSpace( EV[1]*I2 - S );
MatrixScalarMultiply:   "calling external function"
MatrixScalarMultiply:   "NAG"   hw_f06edf
MatrixScalarMultiply:   "calling external function"
MatrixScalarMultiply:   "NAG"   hw_f06edf
SingularValues:   "calling external function"
SingularValues:   "NAG"   hw_f08kef
SingularValues:   "NAG"   hw_f08kff
SingularValues:   "NAG"   hw_f08mef
                            [0.999173371900192842 ]
                           {[                     ]}
                            [0.0406518497192806250]
 
> Rank( EV[1]*I2 - S );
MatrixScalarMultiply:   "calling external function"
MatrixScalarMultiply:   "NAG"   hw_f06edf
MatrixScalarMultiply:   "calling external function"
MatrixScalarMultiply:   "NAG"   hw_f06edf
SingularValues:   "calling external function"
SingularValues:   "NAG"   hw_f08kef
SingularValues:   "NAG"   hw_f08mef
                                       1

ps. All else aside, LinearSolve(A) might be replaced by LinearSolve(A,Vector(2)) if you want to compute the nullspace of A rather than have A get interpreted as an augmented system.

acer

What are the odds that every one of 10 positive integers in a group (row) of 10, all taken from the pool of 1..50, is distinct?

Isn't  it (49/50)*...*(41/50)=0.38=38%? So the odds of at least one duplicate is 62%, no?

acer

I find this example interesting because the OP starts with an operator, and so presumably wishes to end up with an operator as well.

Thomas' reply shows how to get that, using unapply. Great. (He might also have started with the function application r(t) to get the Vector over which he maps diff.)

But why shouldn't D work directly on r?

There are less-well-known routines codegen[GRADIENT], codegen[JACOBIAN], and codegen[HESSIAN] for similar purposes.

Was there not a codegen[DIFF] many releases ago, or am I misremembering? Is D supposed to do that job (now)? If so, should it work here?

acer

Don't use square-brackets or curly braces as delimiters. They are for other things (creating lists and sets, respectively).

You have an expression, not an equation. Perhaps you implicitly consider it to be equal to zero.

The_equation := -(u+delta)*(1-c)*(c-c2h-t)+(1-u)/(1-q^2)/(4-q^2)^2*((2-q^2)*(theta-c2l-t)
                -q*(1-c))^2-(1-u)*(1-q^2)*(1/(4-q^2)/(1-q^2)*(2-q^2)*(theta-c2h-t)-q*(1-c))^2;

new_eq := eval(The_equation,[c2l=1/2,delta=1/2,q=1/2,theta=1/2,u=1/2]);

plots:-implicitplot3d(new_eq,t=0..10,c=0..1,c2h=0..1,axes=boxed);

Alternatively, you could try to isolate new_eq (using solve or isolate) for one of c, c2h, or t in terms of the other two and then use simply plot3d.

acer

Note the full colon in the end do:

for i from 1 to 3 do
   something:=19;
   print(plot(sin(i*x)));
   somethingelse:=13;
end do:

acer

You might consider using a Worksheet rather than a Document for prototyping longer  or more involved bits of code. (I mean, purely computational parts.) Then you can see it all explicitly. You might even consider using 1D input rather than typset 2D Math for any such algorithmic development work. Once working to your satisfaction you might then copy into code-edit or start-up regions of another Document. (And if you plan on re-using such code in multiple Documents, you might consider looking at saving some of your routines to Library archives.)

One thing to be wary of, in either Worksheet or Document, is executing sections or lines out of order. Just because the GUI allows one to go back to an earlier line, to tweak a formula, etc, doesn't mean that doing so is always prudent. Other lines which depend on the given line can get "out of sync" by such actions. (It's not so obvious, but this is another reason why using customized procedures is often much better than laying out lots of "top-level" code.)

acer

If only one feasible point is wanted, then the following can be quick,

> eqs:={20<=a,a<=60, 20<=b,b<=60, 20<=c, c<=60,
>       20<=d,d<=60, a+b+c+d=160, 4*a+6*b+8*c+12*d=1200}:

> Optimization:-Minimize(1,eqs,assume=integer);
                     [1, [a = 40, b = 40, c = 40, d = 40]]
 
> eval(eqs,%[2]);
                 {160 = 160, 1200 = 1200, 20 <= 40, 40 <= 60}

If you want to characterize all the solution space, and can wait, then the following might do,

SolveTools:-Inequality:-LinearMultivariateSystem(eqs,[a,b,c,d]);

acer

Note, when you say "simplify units in 2D math" it appears that you actually mean simplify units using the context-menus.

1) After executing the assignment to m1, and right-clicking, I did get the Units submenu in the pop-up context-menus. (Unfortunately in Maple 13 applying the Units->Simplify menu-item produces new 2D Math output 2*Unit(kg) instead of the properly typeset 2*[[kg]] as it did in Maple 12. A bug.)

2) For the equation F=m1*a the context-menu system is not providing Units->Simplify. I suspect that this is because it is only coded to detect an algebraic object (with units). It would indeed be an improvement, if it also detected an equation where either rhs or lhs had units attached.

3) If you have assigned to symbol m then you can indeed unassign by issuing either of these commands (as you tried, and where in fact that unassignment succeeded).

> m := 'm':

> unassign('m');

The 1D Math input command Unit('m') should work even when symbol m has been assigned a value. One combination that does not seem to work is using the [[m]] palette entry in a Document following assignment to the global name m.

But one can still get the nicely typeset [['m']] as 2D Math with the quoted name 'm' in the input by using command-completion. That will display as [[m]] in the output, without quotes. Suppose that I type in 45*Unit<Ctl-shift-spacebar> to get command completion suggestions in a pop-up. I choose "unit" which is the top suggestion. Then for the highlighted x I type in 'm'.

Note. The exact keystrokes to get command-completion suggestions may vary with platform. On Linux, I used Ctl-shift-spacebar.

You tried to use Units:-Natural at one point. I strongly advise against this. That subpackage's environment makes all these namespace issues much worse. Also, another reason it didn't work for you because your Document restarts after loading that package. The restart clears whatever packages are loaded.

4a) I don't understand what you mean about "a value included in the numerator", sorry.

4b) To separate an expression and its units, one can use the command convert,unit_free. It may be appropriate to first apply combine,units. For example,

> restart:
> F := 45*Unit(kg)*Unit(m)/Unit(s^2);
                                    45 [m]
                               F := ------ [kg]
                                     [ 2]
                                     [s ]
 
> F_value := convert(combine(F,units),unit_free,F_units);
                                 F_value := 45
 
> F_units;
                                      [N]

Note. You might consider loading the Units:-Standard package, which would bring about some automatic combining & simplification of units for arithmetic done at the top-level. Even if you choose not to load that package, there are commands to manipulate several aspects of the units in expression, which might serve even when the smaller subset of functionality in the right-click context-menus isn't working out for you.

acer

You'll need to compile and link as .dll (or .so on Linux, etc).

See define_external and related help.

acer