Maybe this bug could be alleviated by changing some instance of a call to `normal` inside the internal routines used by JordanForm to instead be a call to `Normalizer`. (I didn't debug in order to figure out for which lines, if any, this might be true.)

> restart:

> unprotect(normal):
> _orig_normal:=eval(normal):
> normal:=x->Algebraic:-Normal(x):

> A:=Matrix([[1,sqrt(2),-sqrt(6)],[sqrt(2),2,sqrt(3)],[-sqrt(6),sqrt(3),0]]):

> LinearAlgebra:-JordanForm(A);
                                [-3    0    0]
                                [            ]
                                [ 0    3    0]
                                [            ]
                                [ 0    0    3]
 
> Q1:=LinearAlgebra:-JordanForm(A, output='Q'):
> Q2,J2:=LinearAlgebra:-JordanForm(A, output=['Q','J']):

> normal:=eval(_orig_normal):
> protect(normal):

> simplify(Q1);
                 [             1/2  1/2             1/2  1/2]
                 [            2    3               2    3   ]
                 [ 1/3      - --------- + 2/3    - ---------]
                 [                2                    2    ]
                 [                                          ]
                 [   1/2           1/2                      ]
                 [  2             2                         ]
                 [- ----          ----                0     ]
                 [   6             6                        ]
                 [                                          ]
                 [  1/2              1/2                    ]
                 [ 6                6                       ]
                 [ ----         1 - ----              1     ]
                 [  6                6                      ]
 
> simplify(Q2);
                 [             1/2  1/2             1/2  1/2]
                 [            2    3               2    3   ]
                 [ 1/3      - --------- + 2/3    - ---------]
                 [                2                    2    ]
                 [                                          ]
                 [   1/2           1/2                      ]
                 [  2             2                         ]
                 [- ----          ----                0     ]
                 [   6             6                        ]
                 [                                          ]
                 [  1/2              1/2                    ]
                 [ 6                6                       ]
                 [ ----         1 - ----              1     ]
                 [  6                6                      ]
 
> simplify(Q2^(-1) . A . Q2 - J2);
                                 [0    0    0]
                                 [           ]
                                 [0    0    0]
                                 [           ]
                                 [0    0    0]

The idea is that, for a Matrix with radicals present, one might first assign `Normalizer` to be a routine capable of handling radicals (eg. `radnormal`, or Algebraic:-Normal, or `simplify`).

acer

Maybe this bug could be alleviated by changing some instance of a call to `normal` inside the internal routines used by JordanForm to instead be a call to `Normalizer`. (I didn't debug in order to figure out for which lines, if any, this might be true.)

> restart:

> unprotect(normal):
> _orig_normal:=eval(normal):
> normal:=x->Algebraic:-Normal(x):

> A:=Matrix([[1,sqrt(2),-sqrt(6)],[sqrt(2),2,sqrt(3)],[-sqrt(6),sqrt(3),0]]):

> LinearAlgebra:-JordanForm(A);
                                [-3    0    0]
                                [            ]
                                [ 0    3    0]
                                [            ]
                                [ 0    0    3]
 
> Q1:=LinearAlgebra:-JordanForm(A, output='Q'):
> Q2,J2:=LinearAlgebra:-JordanForm(A, output=['Q','J']):

> normal:=eval(_orig_normal):
> protect(normal):

> simplify(Q1);
                 [             1/2  1/2             1/2  1/2]
                 [            2    3               2    3   ]
                 [ 1/3      - --------- + 2/3    - ---------]
                 [                2                    2    ]
                 [                                          ]
                 [   1/2           1/2                      ]
                 [  2             2                         ]
                 [- ----          ----                0     ]
                 [   6             6                        ]
                 [                                          ]
                 [  1/2              1/2                    ]
                 [ 6                6                       ]
                 [ ----         1 - ----              1     ]
                 [  6                6                      ]
 
> simplify(Q2);
                 [             1/2  1/2             1/2  1/2]
                 [            2    3               2    3   ]
                 [ 1/3      - --------- + 2/3    - ---------]
                 [                2                    2    ]
                 [                                          ]
                 [   1/2           1/2                      ]
                 [  2             2                         ]
                 [- ----          ----                0     ]
                 [   6             6                        ]
                 [                                          ]
                 [  1/2              1/2                    ]
                 [ 6                6                       ]
                 [ ----         1 - ----              1     ]
                 [  6                6                      ]
 
> simplify(Q2^(-1) . A . Q2 - J2);
                                 [0    0    0]
                                 [           ]
                                 [0    0    0]
                                 [           ]
                                 [0    0    0]

The idea is that, for a Matrix with radicals present, one might first assign `Normalizer` to be a routine capable of handling radicals (eg. `radnormal`, or Algebraic:-Normal, or `simplify`).

acer

It is mostly explained above. If the Matrix has the symmetric or hermitian indexing function, and if the data is all of type complex(numeric) with some float present (or complex[8] or float[8] datatype, etc), then the underlying method used by the Eigenvectors command will sort the eigenvalues in the returned Vector.

And the eigenvectors are columns in the returned Matrix which correspond to entries in that returned Vector of sorted eigenvalues, so the columns are thus sorted for these cases.

This is all due to which external clapack routine is used.

The key thing is that the Eigenvectors and Eigenvalues commands do not try to automatically detect symmetry -- the symmetric of hermitian indexing function must be specificed in the data Matrix.

I think that that is ok behaviour, because with the lack of a nonsymmetric indexing function there would otherwise be no way to force use of the nonsymmetric method in the case that float data only appeared close to hermitian. An automatic test would likely involve a tolerance in the float case.

acer

It was bothering me, that I got the wrong result when I first tried to use Quantile for this. (This is what Quantile is for, in general. So it's a problem if it is weak at this.) Here's a bit of analysis.

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

> # This next goes wrong because `evalf/RootOf` doesn't tell fsolve to use a real range!
> evalf[50](Quantile(X,1-10^(-18)));               -3.9377187255341233799589551305177585550470808666755

            - 4.5754487927888755526460423197296824251197610676463 I

> # Observe that for exact arguments Quantile returns a RootOf that lacks "realness" info.
> # Really, that is a bug of Quantile (certainly not fsolve's fault).
> forget(evalf):
> T:=Quantile(X,1-10^(-18));
                                                                (1/2)
       RootOf(500000000000000000 erf(_Z) - 499999999999999999) 2     

> # Sometimes this next call (which gets to `fsolve`) will "get lucky"
> # and find the real root. But often not!
> evalf(%);
                                 8.757290342

> # Let's fix that, by injecting a real-range into the RootOf
> forget(evalf):
> subsindets(T,specfunc(anything,RootOf),t->RootOf(op(t),1..20));evalf[50](%);
                                                                     (1/2)
   RootOf(500000000000000000 erf(_Z) - 499999999999999999, 1 .. 20) 2     
             8.7572903487823150638811286221420828183378493880373

Now for some easier general purpose workarounds (read: better use of Statistics).

# Gosh, silly me. Statistics can do its own real root-finding,
# but we have to let it know that it's a float problem up front.
> forget(evalf):
> Quantile(X,1.0-10^(-18)); # notice 1.0 instead of exact 1
                                 7.391904314

> # And at default precision that answer is poor.
> # But it seems reasonable to expect to have to raise Digits above
> # 18 (at least), if we throw terms like 1.0-10^(-18) around.
> Digits:=20:
> Quantile(X,1.0-10^(-18));
                            8.7572903487832079711

> # Similarly, the 'numeric' option also tells Quantile to use floats.
> forget(evalf):
> Digits:=20:
> Quantile(X,1-10^(-18),numeric);
                            8.7572903487832079711

I ought to have come up with that last pair before.

acer

It was bothering me, that I got the wrong result when I first tried to use Quantile for this. (This is what Quantile is for, in general. So it's a problem if it is weak at this.) Here's a bit of analysis.

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

> # This next goes wrong because `evalf/RootOf` doesn't tell fsolve to use a real range!
> evalf[50](Quantile(X,1-10^(-18)));               -3.9377187255341233799589551305177585550470808666755

            - 4.5754487927888755526460423197296824251197610676463 I

> # Observe that for exact arguments Quantile returns a RootOf that lacks "realness" info.
> # Really, that is a bug of Quantile (certainly not fsolve's fault).
> forget(evalf):
> T:=Quantile(X,1-10^(-18));
                                                                (1/2)
       RootOf(500000000000000000 erf(_Z) - 499999999999999999) 2     

> # Sometimes this next call (which gets to `fsolve`) will "get lucky"
> # and find the real root. But often not!
> evalf(%);
                                 8.757290342

> # Let's fix that, by injecting a real-range into the RootOf
> forget(evalf):
> subsindets(T,specfunc(anything,RootOf),t->RootOf(op(t),1..20));evalf[50](%);
                                                                     (1/2)
   RootOf(500000000000000000 erf(_Z) - 499999999999999999, 1 .. 20) 2     
             8.7572903487823150638811286221420828183378493880373

Now for some easier general purpose workarounds (read: better use of Statistics).

# Gosh, silly me. Statistics can do its own real root-finding,
# but we have to let it know that it's a float problem up front.
> forget(evalf):
> Quantile(X,1.0-10^(-18)); # notice 1.0 instead of exact 1
                                 7.391904314

> # And at default precision that answer is poor.
> # But it seems reasonable to expect to have to raise Digits above
> # 18 (at least), if we throw terms like 1.0-10^(-18) around.
> Digits:=20:
> Quantile(X,1.0-10^(-18));
                            8.7572903487832079711

> # Similarly, the 'numeric' option also tells Quantile to use floats.
> forget(evalf):
> Digits:=20:
> Quantile(X,1-10^(-18),numeric);
                            8.7572903487832079711

I ought to have come up with that last pair before.

acer

Since there has been exporting to image files in this thread, Robert may have had in mind something close to the trusty, dusty 640x480 (4:3) vga standard graphics mode. See also vesa.

acer

Since there has been exporting to image files in this thread, Robert may have had in mind something close to the trusty, dusty 640x480 (4:3) vga standard graphics mode. See also vesa.

acer

Did you mean you instead pass [result] as the first argument to pointplot, ie. the expression sequence result wrapped in a list.

And did you want to instead create result so that it took the rhs() of f(i)? You want result to be a list of lists, with each inner list having two purely numeric entries (not equations).

result := seq([i, rhs(f(i))], i = 0 .. 25);

plots:-pointplot([result]);

acer

Did you mean you instead pass [result] as the first argument to pointplot, ie. the expression sequence result wrapped in a list.

And did you want to instead create result so that it took the rhs() of f(i)? You want result to be a list of lists, with each inner list having two purely numeric entries (not equations).

result := seq([i, rhs(f(i))], i = 0 .. 25);

plots:-pointplot([result]);

acer

It is a shame that the newer NumericalAnalysis:-IsMatrixShape doesn't behave like LinearAlgebra:-IsDefinite does with respect to nonsymmetric Matrices and the positive-definite check. It might be a recent coding oversight (read: bug) in the NumericalAnalysis package.

As mentioned above, the definition used by IsDefinite is that the scalar quantity computed by x^%T.M.x , for any Vector x not all zero, is strictly greater than 0 for a positive-definite Matrix. That's a fine definition, regardless of whether the Matrix is symmetric.

In the case that the Matrix M happens to be symmetric the eigenvalues will also be [not only real but also] strictly positive. That doesn't hold in the nonsymmetric case, and so isn't the general definition. OK.

Now, any real Matrix M can be split into two parts, one of which is symmetric and the other of which contributes nothing to x^%T.M.x . Since the two parts add up to be equal to M the whole, one need only test the symmetric part of M in order to ascertain whether M is positive definite. That`s what LinearAlgbera:-IsDefinite does.

> M:=Matrix([[1,0],[1/2,1]]):

> (M+M^%T)/2 + (M-M^%T)/2 - M; # two parts, add to equal the whole
                                 [ 0 0 ]
                                 [     ]
                                 [ 0 0 ]

> N:=(M+M^%T)/2: # so-called "symmetric part" of M

> Student:-NumericalAnalysis:-IsMatrixShape(N,positivedefinite); # good
                                    true

> <a,b>^%T . (M-M^%T)/2 . <a,b>; # contributes nothing
                                      0

> LinearAlgebra:-IsDefinite(M); # good
                                    true

> Student:-NumericalAnalysis:-IsMatrixShape(M,positivedefinite); # oh dear
                                    false

acer

In this particular case, I think that `simplify` could be expected to do better.

If there is no `ds` then it can get shortened to a form using both sin(x)^2=1-cos(x)^2 and cos(x)^2=1-sin(x)^2. But that trig simplification gets missed when the terms appear in a coefficient of a polynomial in ds.

> simplify((sin(x)^2-1)*sin(x)*(1-cos(x)^2));
                                     3       2
                              -sin(x)  cos(x) 

> T:=(sin(x)^2-1)*sin(x)*(1-cos(x)^2)*ds+dt;
                 /      2    \        /          2\        
                 \sin(x)  - 1/ sin(x) \1 - cos(x) / ds + dt

> simplify(T); simplify(T,trig); # worth an SCR
                                  2                4        
                 -sin(x) ds cos(x)  + sin(x) cos(x)  ds + dt
                                  2                4        
                 -sin(x) ds cos(x)  + sin(x) cos(x)  ds + dt

> Tsimp:=map(simplify,(collect(T,ds)));
                                 3          2     
                          -sin(x)  ds cos(x)  + dt

> coeff(Tsimp,ds);
                                     3       2
                              -sin(x)  cos(x) 

> coeff(T,ds);simplify(%);
                     /      2    \        /          2\
                     \sin(x)  - 1/ sin(x) \1 - cos(x) /
                                     3       2
                              -sin(x)  cos(x) 

acer

In this particular case, I think that `simplify` could be expected to do better.

If there is no `ds` then it can get shortened to a form using both sin(x)^2=1-cos(x)^2 and cos(x)^2=1-sin(x)^2. But that trig simplification gets missed when the terms appear in a coefficient of a polynomial in ds.

> simplify((sin(x)^2-1)*sin(x)*(1-cos(x)^2));
                                     3       2
                              -sin(x)  cos(x) 

> T:=(sin(x)^2-1)*sin(x)*(1-cos(x)^2)*ds+dt;
                 /      2    \        /          2\        
                 \sin(x)  - 1/ sin(x) \1 - cos(x) / ds + dt

> simplify(T); simplify(T,trig); # worth an SCR
                                  2                4        
                 -sin(x) ds cos(x)  + sin(x) cos(x)  ds + dt
                                  2                4        
                 -sin(x) ds cos(x)  + sin(x) cos(x)  ds + dt

> Tsimp:=map(simplify,(collect(T,ds)));
                                 3          2     
                          -sin(x)  ds cos(x)  + dt

> coeff(Tsimp,ds);
                                     3       2
                              -sin(x)  cos(x) 

> coeff(T,ds);simplify(%);
                     /      2    \        /          2\
                     \sin(x)  - 1/ sin(x) \1 - cos(x) /
                                     3       2
                              -sin(x)  cos(x) 

acer

This is the first parser difference I mention here.

acer

This is the first parser difference I mention here.

acer

It's a bug in the documentation only, I think, that it doesn't describe the nonsymmetric case.

The term "positive definite" refers to what sign the scalar value x^%T.A.x can take on for any nonzero column Vector x, rather than (always) what sign the eigenvalues have. Hence Robert's example below is not a counter-example to match your bracketed/ie description, I think.

Now, some texts state the definition of the term "positive definite" by stating that "A symmetric/hermitian matrix is positive-definite if..." and so of course with that wording the symmetric aspect is built right into the concept. And for a symmetric/hermitian Matrix the positive-definiteness matches whether the eigenvalues are all positive. But the definition needn't be taken as always including the symmetric quality.

Modulo bugs, you ought to be able to rely on IsDefinite for nonsymmetric Matrices.

Having said that, I sort of remember that there was a short-lived (1 release, then fixed, maybe in Maple 12?) glitch introduced in the recent past for the nonsymmetric case. It's Sunday afternoon,... and I hope that I haven't remembered that situation all wrong.

acer