@Carl Love 

2) Yes, if the Matrix does not have the symmetric indexing-function then the type-check against 'Matrix(symmetric)' will, alas, get to some interpreted Library code.

restart;

#showstat( `type/rtable/TwoDim`, 147..161 );

stopat(`type/rtable/TwoDim`, 148):

type( Matrix([[1,2],[2,3]], shape=symmetric, datatype=float[8]),
      'Matrix(symmetric, datatype=float[8])' );

                    true

type( Matrix([[1,2],[2,3]], datatype=float[8]),                 
      'Matrix(symmetric, datatype=float[8])' );
false
`type/rtable/TwoDim`:
 148*      if dims[1] <> dims[2] then
               ...
           end if;

DBG> cont

                    true

There should be lots of ways to speed that up that mechanism, whether float[8] or not.

3) I added the conversion back to Fortran_order just because it's the default for float[8] Matrices, and the LDLT code constructs a new triangular[lower,unit] Matrix already, for the return value. Back in 1999 the accelerated BLAS only supported Fortran_order, and LAPACK (and NAG, at that time) was designed for it, which is why it was made the default for float[8] Matrices in Maple 6. Nowadays the Intel MKL has an additional  C_order API (and reasonable support for it). LinearAlgebra could be improved to call out to that MKL C-interface instead of making Fortran_order copies, when passed C_order Matrices. But until then... it's mostly better to have one's float[8] Matrices be Fortran_order.

@Rouben Rostamian  

restart;

kernelopts(version);
    Maple 2018.1, X86 64 LINUX, Jun 8 2018, Build ID 1321769

s := x -> min(frac(x),1-frac(x)):
b := x -> sum(s(2^n*x)/2^n, n=0..4):

limit( ( b(x)-b(9.05) )/(x-9.05), x=9.05 );
                   0.

evalf(Limit( ( b(x)-b(9.05) )/(x-9.05), x=9.05 ));
               3.000000000

MultiSeries:-limit( ( b(x)-b(9.05) )/(x-9.05), x=9.05 );                  
                    3

foo := convert( ( b(x)-b(9.05) )/(x-9.05), rational, exact ):

lprint(foo);
  (min(1-frac(x),frac(x))+1/2*min(1-frac(2*x),frac(2*x))+1/4*min(1-
  frac(4*x),frac(4*x))+1/8*min(
  1-frac(8*x),frac(8*x))+1/16*min(1-frac(16*x),
  frac(16*x))-17/80)/(x-181/20)

MultiSeries:-limit( foo, x=905/100 );                        
                     3

evalf(Limit( foo, x=905/100 ));
                3.000000000

limit( foo, x=905/100 );
                     0

Using Maple 18.02 and 16.02 I get 3. from :-limit, but in Maple 2015.0 and 2015.2 and later I get 0.

@Carl Love Very nice. Vote up.

A few comments:

1) This relates only to the generation of the example Matrix A.
     For the call to RandomMatrix, it is quite a bit faster to use the syntax
     generator=0. .. 1. instead of generator=rand(0. .. 1).

2) The type-check in A::Matrix(symmetric, datatype= float[8]) on the first parameter of LDLT is, unfortunately, much slower than it ought to be in the case that the first (Matrix) argument happens to be symmetric in its values but does not actually have the symmetric indexing function. So it seems fair to pump in a shaped Matrix to LDLT. I see a timing improvement from about 2.12sec to about 1.56sec by that change of the input.

3) The compiled external function which does the Cholesky decomposition is actually DPOTRF from the Intel MKL, and not a generic LAPACK version such as NAG F07FDF. It is highly optimized to use SIMD or other fancy chipset extensions the MKL detects at runtime. But it is also tuned to minimize cache misses. One aspect of such optimization that can (sometimes) be partially achieved when using either Maple's Compiler or evalhf mode is improving the cache hits. The innermost loop of LDLT is the loop,
    for p to k1 do S:= S + A[i,p]*R[p] od
which accesses sucessive values along a row of the Matrix. Constructing the workspace Matrix AA of LDLT (which is passed as Matrix A of LDLT_internal) to be order=C_order means that its entries are stored by row in memory. This seems to improve the performance. For the 1024x1024 example, with shape=symmetric on the input Matrix, I see the timing improve from about 1.56sec to about 1.21sec.

LDLT_a.mw

@David Sycamore You can utilize commands which are members of packages by using a long form of the name, or by loading the package and then using the short form of the name.

The NumberTheory package is newer than the numtheory package.

Both of those packages allow the long form syntax package:-member as well as the long form package[member] . I prefer the former as simple, safe and unambiguous. The latter can run into name-space issues depending on the content, unless you do more typing to arrive at package[':-member'] .

restart;

NumberTheory:-PrimeCounting(31);
                               11
NumberTheory[PrimeCounting](31);
                               11
numtheory:-pi(31);
                               11
numtheory[pi](31);
                               11

restart;

with(NumberTheory):

PrimeCounting(31);
                               11

PrimeCounting(47);
                               15

restart;

with(numtheory):

pi(31);
                               11

pi(47);
                               15

The LDL^T decomposition is suitable for symmetric indefinite mattices, and does pivoting to handle that.

Cholesky LL^T decomposition is for symmetric positive definite.

LDL^T implementations are usually considerably slower than Cholesky LL^T, which is not unexpected.

The OP wrote that his Matrix is "s.p.d". So why the need for LDL^T? Just because some instructor demands it? Or is there some numeric justification? (eg. close to indefinite, or...?) If there is in fact a numeric justification for preferring LDL^T over LL^T for the given Matrix then performing LL^T and then adjusting would defeat that goal.

I don't recall seeing LDL^T called "Cholesky" by an expert. It seems like an unnecessary and confusing equivocation.

@samen You have incorrectly changed the code I showed above. Your Reply shows,

plots:-animate(K->plots:-matrixplot(p(k)-p(K)),[k],k=0.0001..0.001);

while my code above showed,

plots:-animate(K->plots:-matrixplot(p(0.0001)-p(K)),[k],k=0.0001..0.001);

With your incorrect change it won't work properly for any N.

So first I reverted that change.

Next I replace the call to solve inside procedure p with a call to fsolve. That sped it up enough that N=19 was bearable, but N=49 was still too slow.

Then I replaced the use of fsolve (or solve) by use of LinearSolve, after generating the Matrix-form of the equations in sys, which are all linear in the variables in vars. That made the case N=49 run reasonably quickly. See the attachment.

crnk_ac.mw

The call to GenerateEquations is now the slowest part. At least it's outside the of any procedure like pmat or p, and so only gets called once. But you could avoid that time cost as well if you generated the Matrix-form directly instead of bothering to forms the eqs.

I don't know whether your discretization scheme is overall correct. I've focused only on the aspect of solving the equations you've targeted.

If your eventual goal is to numerically evaluate the eigenvectors at floating-point values for all the parameters then you would likely be better off with a procedure (function) that admitted those float values and did the eigenvector calculations in hardware double precision.

That is likely true whether you intend on operating purely within Maple or in C or Fortran. In the latter cases I mean a purely hardware float LAPACK eigen-solver like from the Intel MKL, say, as opposed to exporting the symbolic formulas for the eigenvectors.

Even just within Maple, an efficiently coded procedure that acted inplace on float[8] Matrices should be able to do millions of such eigenvector computations per second. (In a recent bit of work I was able to crank out several million 5x5 eigenvector computations per second in a highly optimized set of procedures in just Maple, and in C/Fortran there should be even less function-call overhead.)

But there is also numerical stability to consider. The eigen-solvers in LAPACK are quite good, but there's no guarantee that the floating-point instantiation of length symbolic eigenvector formulas would be so numerically robust.

And one more.

map( `^`, L1, 2 );

That should faster than using the user-defined operator u->u^2 for very long lists, since it only uses the kernel builtin `^`. 

[edited] And unlike the elementwise operator ~^ it will work in much older versions of Maple, back to MapleV R5 (released January 1998).

@jthress1 The only difference I see is in the order in which the 3 distinct eigenvalues might appear in the Vector, and in whether the (numerators and denominator) entries in the eigenvectors are factored. But scalar multiples of the eigenvectors would also be acceptable.

1D_Euler_Eigenvectors_ac.mw

I don't see anything to object to here. Just be careful to associate the eigenvectors with the right eigenvalues, and account for the order in which they appear.

And if any eigenvalue were associated with more than one independent eigenvector (not the case in this followup) then various linear combinations would also be acceptable.

@GoBills3 The key thing is your use of OS X.

I have seen this reported before, and alternatives suggested if I'm not mistaken. Searching this site might find the older Questions I recall seeing.

@jthress1 Perhaps you are having trouble expressing your idea mathematically, but this paragraph you wrote above contains logical contradictions and some of the sentences do not make sense.

"The eigenvalues are exactly what I expected. In fact, I know that the eigenvalues are correct. What I expected was for each Eigenvalue to have only one eigen vector associated with it. Since the last four eigenvalues are identical, they should have the same eigenvector(s) as mentioned earlier. What I did not expect is for each eigenvalue to have four eigenvectors associated with it. What i expected was to have six total eigenvectors (one eigenvector for each eigenvalue). Does this clear things up? Sorry for the confusion"

It does not make sense to say that the same eigenvalue (or multiplicity four) can be considered as four eigenvalues each with only one eigenvector, as in this sentence, "What I expected was for each Eigenvalue to have only one eigen vector associated with it".

Generically, the Matrix A as given can be shown to have six linearly independent eigenvectors. But there are only three distinct eigenvalues. One of the eigenvalues is associated with an eigen-space having four linearly independent eigenvectors as its basis.

It makes no sense to say that you expect one of the eigenvalues to have multiplicity greater than one, while both objecting to that eigenvalue's associated eigen-space having rank greater than one as well as expecting a full linearly independent set of eigenvectors for the Matrix altogether.

@jthress1 

Leaving eigenvectors aside for the moment, your symbolic Matrix A appears to have an exact eigenvalue with multiplicity of four. That is, there are only three distinct eigenvalues. Are you saying that is not what you expect? If that contradicts your expectations then perhaps you should resolve that issue before proceeding to eigenvector computation.

For many numeric values of the parameters there can still be six linearly independent eigenvectors, however. It's not clear whether that alone satisfies all your requirements, irrespective of how many distinct eigenvalues your Matrix might have.

@dharr Please see the edit in my followup Comment to my Answer.

It seems to me that the manual solution of the Nullspace (of the CharacteristicMatrix instantiated at each eigenvalue) provides six independent eigenvectors, as computed in my Answer.

It may be only the computation using Eigenvectors -- with forced Normalizer -- that returns the two zero-vectors. Given the presence of radicals in two of the eigenvalues then forcing Normalizer=normal may have been a mistake of mine in the followup..

It seems as if the slow behavior of the LinearAlgebra:-Eigenvectors command (for this example) can be improved by assigning a (straightforward) particular simplifier to Normalizer -- and thus bypassing any "smart" deduction of the normalizer.

Some commands use Normalizer and Testzero to keep down expression swell during row reduction, as well as guard against inadvertant division by hidden zeros while pivoting.

For this example plain old normal is adequate to obtain reasonably short expressions for the eigenvectors.

[edited]
Member dharr points out that two columns of evecsA as computed in the attachment below -- by the Eigenvectors command with Normalizer=u->normal(u) -- are all-zero, which seems to indicate a dearth of independent eigenvectors and a defective Matrix. Yet in the earlier Answer I gave above there appear to be 6 independent non-zero Vectors computed, which each satisfy the definition appropriately. This suggests that the zero-vectors provided by Eigenvectors in the attachment below are incorrect, and those in the Answer above are correct.
Download eigs_LA.mw
[end of edit]

@jthress1 No, the eigenvectors associated with each of the evalsA[ii] are in the sets assigned to evecsA[ii].

There are 4 eigenvectors associated with one of the eigenvalues, which is the same eigenspace for the evalsA[ii], ii=3..6 which are all equal. On the other hand, there is just one eigenvector associated with evalsA[1], and also just one associated with evalsA[2].

That last part of the worksheet is just a verification that the (associated pairs of) eigenvalues and eigenvectors satisfy the definition:
   A . evec = lambda*evec
which all produce the expected zero vectors.

I computed the eigenvectors evecsA[ii] associated with each eigenvalue evalsA[ii] by computing the nullspace (kernel) of the characteristic Matrix (ie. A-evalsA[ii]*IdentityMatrix).