You say

There are still some few errors at the end of the program

Once the reference to [px, py] is removed (as in my previous post). There are no errors!

I have added some plots at the end of the attached worksheet - where are the errors you refer to??

CRSpline3.mw

Every time you post a question and I supply the code which fixes your problem, you immediately throw away this code and produce some more crap which dosn't work. Hence I will not be supplying you with any more code. It is a waste of my time. Your latest worksheet produces the error message

Warning, expecting only range variable t in expressions [px, py] to be plotted but found names [px, py]

Not really surprising since you have a statement

plot1 := plot({[px,py,t=0..6.99]}):

but px, py are defined nowhere in your worksheet, so you are asking for the plot of two unknown quahtities, pretty stupid really.

Ok - at great risk to my own credibility, I'd like to post a reality check here.

You claim that you want to deal with "square matrix under consideration is of size more than 2^30 times 2^30". In other words your matrix will have ~1.0*10¹⁸ entries. If I assume that you can use as few as 32-bits/entry, then the back of my cigarette pack says that you will need 2.3million, 2TB hard drives, just to store the matrix entries. I've seen some seriously big server farms in my time.....but this is definitely in the *ouch* category!!!

Possibilities

1. You may have a seriously sparse matrix
2. You may be able to generate entries "on-the-fly", and by some appeal to the structure of the resulting matrix, compute eigenvalues from partial matrix entries
3. Otherwise, you are f****d

On the specific subject of the algorithm which you have posted:

1. You should be aware that Maple (like Matlab and Mathematica) uses industry-standard (and free) librariess for solving problems in linear algebra. These libraries are known as lapack and blas (check wikipedia): none of these commercial software packages really want you to know that they are all using the same underlying routines - but they are!!
2. When you are using maple/matlab/mathematica, for linear algebra problems, what you essentially have is an "interface" to lapack/blas. For small problems, the "interface" time may be significant, but for large (linear algebra) problems, performance will be dominated by the underlying lapack routines, which is why each of these commercial packages will return essentially the same performance (whihc of course, no-one will admit)
3. When you are computing eigenvalues/eigenvectors in maple, the default commands - as in LinearAlgebra[Eigenvalues] or LinearAlgebra[Eigenvectors] will be calling the underlying lapack/blas library routines. No algorithm that you can write will be faster than the lapack/blas routines: so don't bother.
4. If the "Francis double step QR algorithm" is the fastest way to solve your problem, then lapack will figure this out and use it: or it will use something better/faster. In other words, the best solution for your example is

restart;
with(LinearAlgebra):
A := Matrix( [  [7, 3, 4, -11, -9, -2],
                      [-6, 4, -5, 7, 1, 12],
                      [-1, -9, 2, 2, 9, 1],
                      [-8, 0, -1, 5, 0, 8],
                      [-4, 3, -5, 7, 2, 10],
                      [6, 1, 4, -11, -7, -1]
                  ]
               ):
evalf(LinearAlgebra[Eigenvectors](A));

Since sum() attempts to find a closed form solution, and then perform the summation, this can occasionally be quicker than using a brute-force add(). For example, compare

CodeTools[Usage](add(j, j=1..100000000));

CodeTools[Usage](sum(j, j=1..100000000));

the latter is much quicker.

However for summations with relatively few terms, or where no closed-form solution for the sum() can be obtained, or where trying to figure out the closed form solution takes a while,  then add() will be quicker

If you wish to persist with the 2D math input, then , in the expression for the derivative, you cannot have a space between ∂y and its exponent - so ∂y ²∂t will be parsed incorrectly, but ∂y²∂t will work

pp:=n[0](t-i+1)*M[i-1]+n[1](t-i+1)*M[i]+n[2](t-i+1)*M[i+1]+n[3](t-i+1)*M[i+2]; 
return evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp,
i=1..Num-2
)
):

Normally, I guess these lines enable the definition of the curve with the Catmull-Rom splines. But, I didn't understand the operation of it namely the 2 following points :
- why it used the index (t-i+1) in the defintion of pp ?
- how the return function works? I have seen the help, but i didn't seen any examples which makes me difficult to understand for me.

Well, it doesn't "use the index (t-i+1) in the defintion of pp".because (t-i+1) is not an "index" -it is an argument

Each of the n[] is a function which may be evaluated for any argument. So think of (say) n[0] as a function - as a "toy" example, suppose n[0] is x->x^2+1, then n[0](t-i+1) would return (t-i+1)^2+1.

At the point in the code where pp is evaluated, neither 't' nor 'i' is defined, so (for example) n[0](t-i+1) will return an expression in terms of 't' and 'i', rathr than any numeric value. In your code, each of the M[] values is a simple two element list of floats so n[](t-i+1)M[] returns an expression in t and i, which multiplies a list of two floats. As a toy example, suppose that n[0]:=x->x^2+1 and M[i]:=[0.1, 0.2], then n[0](t-i+1) would be ((t-i+1)^2+1)*[0.1,0.2], which you could also interpret as [0.1*((t-i+1)^2+1),0.2*((t-i+1)^2+1)]

At the point in your code where pp is evaluated, nether t nor i is known, so pp will be a (moderately) complicated expression involving both 't' and 'i'.

The expression evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp, i=1..Num-2 will essentially substitute integer values for 'i' in the expression (1+signum(-t+i))*(1+signum(t-i+1))/4*pp, including values for 'i' in pp, and then sum them: this will return a (moderately complicated) expression in terms of t. The evalm() part just ensures that if you have expression1*[a,b]+expression2*[c,d], then this will be combined as [expression1*a+expression2*c, expression1*b+expression2*d], so that the result of

evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp, i=1..Num-2

will be a two-element list whihc could be summarised as

[ complicated expression of t, another complicated expression of t].

All that the return statement does is to "return" this 2-element list of more-or-less complicated expressions. So ultimately, the procedure returns a 2-element list [x(t), y(t)], and in subsequent plot commands, for given values of 't', you are plotting y(t) values against x(t) values

#############################################################

I think that the use of the function "piecewise" should be adapted to the definition of this polynomial curve with Catmull-rom splines.

May you help me to implement this polynomial curve with a piecewise function in Maple ?

This statement probably demonstrates your complete and utter misunderstanding of the use of splines (of any knid) for fitting data. The objective of using splines is to produce a continuous function which is valid for all the supplied data. Having generated a continuous function which is valid for the domain of your domain - please exaplin why you wish to split this continuous function into a piecewise one

I have fixed another pile of syntax errors

It would be helpful if you tell the difference between the lower case letter 'l', and the integer '1', Similarly the upper case letter 'O' is not the same as the integer '0'

I have made suggestions about fixing other logical errors in your procedure, but since I don't really know what you are trying to achieve - I'm largely guessing - see the attached. At least all section of the code seem to execute, although I am pretty sure the answers are not what you want.

I have copied/executed the original textbook code as a new section in the attached - this runs, although I have no way of telling whether the answers it produces are correct

CRSpline2.mw

The following "works" (ie returns a numeric value) for nvals=1,2,3,4

restart;
nvals:=3:
int(  exp( -add(x[i],i=1..nvals)^3),
      seq(x[i]=0..1, i=1..nvals),
      numeric=true
   );

However for nvals=5 (or greater), the result is returned unevaluated - interesting - and I'm still thinking

Fixed numerous, minor, syntax errors, indicated in the attached.

Could not complete fixes because OP's worksheet calls the function/procedure crom_2d(), which is not defined in the supplied workheet

CRSpline.mw

Running OP's script in Maple2015.2 on Win7, 64-bit, with no restart() between cases, the reported memory allocation according to the status bar, goes in the following steps

Case              Reported Memory Allocation (MB)

10^4             20.13

10^5             20.13, 38.18

10^6             38.18, 54.18, 74.19, 90.19, 106.19,

                     122.19, 250.19, 266.19, 282.19, 298.19

So the "effect" would seem to be real.

My desktop has 16GB of RAM, and the total allocated memory never got above 6GB, so this particular memory usage is not really an issue for me.

Maybe Maple's garbage collection only runs when allocated memory starts to approach "available" memory. This would seem to be the time-efficient approach. The memory-efficient approach which would be to run gc, as often as possible, ignoring the time overhead?.

On the other hand if I create the loop as a procedure, so that I can use CodeTools, and then rerun the same thing with restarts, between each case, then CodeTools[Usage] reports

Using CodeTools[Usage] for the three cases

          n      memory used   alloc change  cpu time  real time    gc time

       10^4      0.70GiB          0 bytes         7.16s       7.17s      249.60ms

       10^5      7.03GiB          100MiB       70.04s      70.02s        2.11s

       10^6    70.27GiB          0.51GiB      11.88m     11.67m     29.86s

1. In order to demonstrate root values, you are going to have to increase the accuracy of the calculation. Digits:=50 is barely adequate, but Digits:=100 produces vanishingly small residuals
2. The "root" values which you supply, ie  0.165237712988657e-1, .103583272213766  and  .290071279318035 are not roots, whatever you might think.

See the attached

Roots.mw

Like Rouben has said, your request is very unclear: consider for example

1. You say that you want to "be able to do fast fourier transforms". Well, the fast fourier transform is an algorithm for computing the discrete fourier transform - but I'm guessing that you really, really do not want to do this. In any case it would be rather pointless for the SHM
2. You say that "I can then graph" - well no you can't. In order to graph anything, you would have to be able to evaluate everything to numeric quantities. For the SHM example, thsi would ean having two initial conditions which evaluate to numbers (and a numeric value for the natrual frequency). Otherwise no graph is possible

Having said the above the attached worksheet shows very simple (non-numeric) solutions for the SHM - what more do you want/need?

SHM.mw

Certain vertical lines ...

Since you have opened a new thread on this topic, I assume that this question is no longer relevant

So far as I am aware there is no way to automatically indent code in the standard worksheet, unless you use codeEditRegion(s) where a limited(?) capability is available

I rarely uses codeEditRegions because the result is bound to the worksheet in which it is created - ie not transportable across worksheets.

My "default" coding approach is to write a lot of procedures in an external editor (which has all the usual features - bracket matching, indentation etc). Because I have been doing it this way for a while, I find it much easier to read/understand code, when it is indented in the way I expect.

When responding to posts on this site, I often have to read OP's code very slowly/carefully to understand what is going on. Analysing other people's code for other people's problems is much harder than reading/understanding code I have written myself for my own work.

I have adopted the habit of manually indenting OP's code as part of this reading process - it actually helps my understanding (and doesn't slow me down much!)

I think you need to provide some clarification of exactly what you are trying to achieve. In particular your statement

how to separate the ODE system into disjoint cases

I may be wrong, but this *reads* as if you want to change the essential form of the ODE system, based on events. As the relevant help page states very clearly

Note that events are restricted to operations that act on variable values, and cannot be used to change the essential form of the ODE system or exchange equations.

Consider the bouncing ball problem on the same help page: the same ode system is valid for all time, but events[] are used to immediately reset diff(y(t),t) to -diff(y(t),t) when y(t)=0 - in other words to reverse the sign of the velocity when the ball hits the ground. An obvious way to interpret this is

1. From t=0, solve the diffential equation subject to the initial conditions diff(y(t),t,t)=-9.8, y(0)=1, D(y)(0)=0.
2. When y(t)=0, at time= te (say), then solve the same diffential equations subject to the "initial" (ie at t=te) conditions, diff(y(t),t,t)=-9.8, y(te)=0, D(y)(te)=-D(y)(te)

The solution looks disjoint, but this has been achieved by operating on the variable values, not changing the essential form of the of the ODE system