The recolor procedure from Joe Riel's post above can be written also as

recolor:=(p,c)->eval(p,COLOUR=proc() 
    convert(color=c,PLOToptions)[] end);

In case if the plot structures contain COLOR instead of COLOUR, COLOR can be used instead of COLOUR, or both of them can be added in the eval, as [COLOUR=..., COLOR= (the same)].

Alec

IntegrationTools:-Parts(Int((z-x)*diff(F(x),x), x=L...z-1),z-x);
 
                                          z - 1
                                         /
                                        |
             F(z - 1) - F(L) (z - L) -  |       -F(x) dx
                                        |
                                       /
                                         L

IntegrationTools:-Parts(Int(x*diff(F(x),x),x=z..H),x);

                                         H
                                        /
                                       |
                    F(H) H - F(z) z -  |   F(x) dx
                                       |
                                      /
                                        z

Alec

The end of the animation itself, perhaps, gives a good idea about the PDF,

plot(eval(pdf, Dirac(10+x)=
    Statistics:-PDF(Normal(-10,0.01),x)),
    x=-20..40,y=-0.001..0.1,axes=boxed);

Alec

The last one, with square brackets, works for me, both in Classic and Standard Maple,

plot([exp(-t) + t, exp(-t) - t, t=-2..2]);  

Alec

Approximating the Dirac function with Normal PDFs, one can construct an animation. Something like

f:=piecewise(x<-10,0,Statistics:-CDF(Normal(0,10),x));

                 {           0                   x < -10
                 {
            f := {                  1/2
                 {               x 2
                 { 1/2 + 1/2 erf(------)        otherwise
                 {                 20

pdf:=diff(convert(f,Heaviside),x);

                                         2
                                        x     1/2
                                  exp(- ---) 2    Heaviside(10 + x)
                                        200
  pdf := 1/2 Dirac(10 + x) + 1/20 ---------------------------------
                                                  1/2
                                                Pi

                      1/2
                   x 2
         + 1/2 erf(------) Dirac(10 + x)
                     20

plots:-animate(plot,
[eval(pdf,Dirac(10+x)=Statistics:-PDF(Normal(-10,epsilon),x)),
x=-20..40,y=0..0.1],epsilon=0.1..0.01);

Alec

My guess is that it will be the usual PDF multiplied by the piecewise function which is 1 for x ≥ -10 and 0 for x < -10, plus the Dirac(x+10) multiplied by erfc(sqrt(2)/2)/2 = 0.158655254 .

Alec

A simple way is to apply max to solve,

max(3.241106003, -1.384090609);

                             3.241106003

Something like

plot(max(solve(x=y^2,y)),x=0..1);

To plot both solutions, one can include them in a list,

plot([solve(x=y^2,y)],x=0..1);

Specifying the positive range for y would plot only the positive one,

plot([solve(x=y^2,y)],x=0..1,y=0..1);

Alec

If you add output=list in the Eigenvectors command, the output will be a list containing the basis of the eigenspaces, similar to NullSpace. A set of vectors can be converted to a list of vectors and then to a Matrix. For example, for the matrix from the NullSpace help page,

A := <<6,3,0>|<4,2,0>|<2,1,0>>;

                               [6    4    2]
                               [           ]
                          A := [3    2    1]
                               [           ]
                               [0    0    0]

kern := NullSpace(A);

                                [-1/3]  [-2/3]
                                [    ]  [    ]
                       kern := {[ 0  ], [ 1  ]}
                                [    ]  [    ]
                                [ 1  ]  [ 0  ]

Matrix([kern[]]);

                            [-1/3    -2/3]
                            [            ]
                            [ 0       1  ]
                            [            ]
                            [ 1       0  ]

L:=Eigenvectors(A,output=list);

                          [-1/3]  [-2/3]            [2]
                          [    ]  [    ]            [ ]
            L := [[0, 2, {[ 0  ], [ 1  ]}], [8, 1, {[1]}]]
                          [    ]  [    ]            [ ]
                          [ 1  ]  [ 0  ]            [0]

map(x->Matrix([op(x[3])]),L);

                         [-1/3    -2/3]  [2]
                         [            ]  [ ]
                        [[ 0       1  ], [1]]
                         [            ]  [ ]
                         [ 1       0  ]  [0]

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

                            [-1/3    -2/3]
                            [            ]
                            [ 0       1  ]
                            [            ]
                            [ 1       0  ]

Finally, one can write a procedure for either of these methods. Using NullSpace, that can be done as

Eigenspace:=(A,lambda)->Matrix([LinearAlgebra:-NullSpace(
    A-lambda*LinearAlgebra:-IdentityMatrix(op(1,A)))[]]):

Eigenspace(A,0);

                            [-1/3    -2/3]
                            [            ]
                            [ 0       1  ]
                            [            ]
                            [ 1       0  ]

Alec

Here is a reference for the general case, which has some other references,

Tsoy-Wo Ma, Higher Chain Formula proved by Combinatorics, The Electronic Journal of Combinatorics 16 (2009), #N21

Alec

Or op([2,1], ...

op([2,1],ifactor(3^41+1))*op([2,1],ifactor(3^43+1));

                 22202243091419939970046193941356769

Also, numtheory:-factorset can be used instead of ifactor.

numtheory:-factorset(3^43+1);

                      {2, 82064241848634269407}

numtheory:-factorset(3^41+1);

                     {2, 33703, 270547105429567}

%[-1]*%%[-1];

                 22202243091419939970046193941356769

Alec

You could add one more variable, say z, one more equation, z=x-y, and one more range, z=0..1. For example,

eq1,eq2:=x^2+y^2=0.6,(1-x)^2+(1-y)^2=0.6:

fsolve({eq1,eq2,z=x-y},{x=0..1,y=0..1,z=0..1});

        {x = 0.7236067977, y = 0.2763932023, z = 0.4472135955}

Alec

For practical purposes, including export to Excel, it is better, probably, to drop the variable names and store solutions in a Matrix, with putting temp values in the first column, x1 values in the second column etc. In this case. one wouldn't have to worry about such things as scientific notation, using ExcelTools:-Export.

For example, if the solutions are named sol[1], ..., sol[7000], that can be done as

L:=[temp,x1,x2,y1,y2]:

M:=Matrix([L,seq(eval(L,sol[i]),i=1..7000)]);

Alec

To get rid of the scientific notation, and simultaneously to order the solutions, one can use sprintf,

sol:={x2 = 0.2397737599, x1 = 0.7633930119*10^(-14), 
    y1 = 1402.995860, y2 = -47636.68497, temp = 234.9000000}:

sprintf("temp=%f, x1=%f, x2=%f, y1=%f, y2=%f",
    eval([temp,x1,x2,y1,y2],sol)[]); 

  "temp=234.900000, x1=0.000000, x2=0.239774, 
    y1=1402.995860, y2=-47636.684970"

sprintf("temp=%f, x1=%.20f, x2=%f, y1=%f, y2=%f",
    eval([temp,x1,x2,y1,y2],sol)[]);

  "temp=234.900000, x1=0.00000000000000763393, 
    x2=0.239774, y1=1402.995860, y2=-47636.684970"

Alec

Sets were unordered prior to Maple 12 in contrast to lists. They are ordered in Maple 12 and Maple 13. So the simplest solution is to upgrade your Maple.

Or you can convert your sokutions to lists and sort them. For example,

sol:={x2 = 0.2397737599, x1 = 0.7633930119*10^(-14), 
    y1 = 1402.995860, y2 = -47636.68497, temp = 234.9000000}:

T:=table([temp=1,x1=2,x2=3,y1=4,y2=5]):

sort([sol[]],(x,y)->T[lhs(x)]<=T[lhs(y)]);

                                           -14
  [temp = 234.9000000, x1 = 0.7633930119 10   , x2 = 0.2397737599,

        y1 = 1402.995860, y2 = -47636.68497]

Alec

Random correlation matrix of size n×n can be considered as the Gram matrix of n random unit n-dimensional vectors. This way, it can be generated rather simple and fast,

RandomCorrelationMatrix:=proc(n)
    local C, M, i;
    uses Statistics, LinearAlgebra, ArrayTools;
    M:=Sample(Normal(0,1),n^2);
    for i from 0 to n-1 do 
        Normalize(Alias(M,n*i,[n]),2,inplace) 
    od;
    C:=Alias(M,[n,n]);
    Matrix(C^%T.C,shape=symmetric)
end;

For example,

R:=RandomCorrelationMatrix(1000);

                      [ 1000 x 1000 Matrix         ]
                 R := [ Data Type: float[8]        ]
                      [ Storage: triangular[upper] ]
                      [ Shape: symmetric           ]
                      [ Order: Fortran_order       ]

min(LinearAlgebra:-Eigenvalues(R));
                                             -7
                      0.381818645382836018 10

I used the minimal eigenvalue instead of IsDefinite because couldn't wait until IsDefinite stopped and interrupted it.

The values on the diagonal are slightly different than 1, even for 1×1 matrices,

RandomCorrelationMatrix(1);

                        [1.00000000069567996]

It looks as if Normalize was using lower precision for float[8] numbers than it should. Another bug?

Formally, one can argue that it is not a bug because the error is lower than acceptable by the Digits setting - and it is less after setting

Digits:=trunc(evalhf(Digits));

                             Digits := 15

but I still think that it is a bug, because otherwise what is the point in using float[8] datatype? And setting UseHardwareFloats to true doesn't fix that.

By the way, looking trough the corresponding LinearAlgebra:-LA_Main code, I've noticed an interesting detail - several procedures, including LinearAlgebra:-LA_Main:-VectorNorm, use floor(evalhf(Digits)) instead of trunc(evalhf(Digits)). What is the reason for that?

time(seq(trunc(1.0),j=1..10^6));

                                0.421

time(seq(floor(1.0),j=1..10^6));

                                30.061

Alec