@Aakanksha The code I posted runs as expected in the Standard GUI of my Maple 9.50. Does it not run for you? (The plots take about 20 seconds or so.)

Here is the second plot, in which the real component is red and the imaginary component is blue. This was done at default working precision, and I get pretty much the same at Digits=300.

@Aakanksha It does not seem that your `capacity` is purely real-valued for the range L=1..4. That kind of thing will prevent `plot` from being able to produce a single curve for your expression, as is.

If the code below I delayed assignment of the various parameters with floats, and kept the upper index of the summation as `NN` an unknown name so that it could do symbolic summation. (The results seem to concur with those obtained using numeric upper index 100, with `add` or `sum`.) This made `capacity` faster to plot, and more compact.

The code produces a red curve for the real part, and a blue curve for the imaginary part of `capacity`.

Do you have some reason -- physical grounds of the problem, say -- to think that `capacity` should be real-valued with L=1..4 ?

The float values given to `k` and `L` get into the various exponents in your expression, so I suppose then that you really do want principal roots and not real-valued surds.

restart:

N1:=sum(((yo^((b-a+2*p-1)/2)*z^((b-a+2*p-1)/2)*GAMMA((1-(b-a+2*p))/2))
        /(p!*GAMMA(p-a+1)*GAMMA(1+((1-(b-a+2*p))/2)))),p=0..NN):
N2:=sum(((yo^((b+a+2*p-1)/2)*z^((b+a+2*p-1)/2)*GAMMA((1-(b+a+2*p))/2))
        /(p!*GAMMA(p+a+1)*GAMMA(1+((1-(b+a+2*p))/2)))),p=0..NN):
N1:=eval(N1,NN=100):
N2:=eval(N2,NN=100):

x:=(Pi*csc(a*Pi))/(GAMMA(m*L)*GAMMA(k)):

c:=x*(N1-N2):

cc:=log[2](1+1/c):

me:=MeijerG([[],[1-((b+1)/2)]],[[-(b+1)/2,a/2,-a/2],[]],(z*yo)):

m1:=((z*yo)^((b+1)/2))/(GAMMA(m*L)*GAMMA(k)):
m:=2; k:=24.503; a:=k-m*L; b:=k+m*L-1;
z:=(k*m)/10^(0.1*15):

yo:=10^(0.1*10):

pp:=m1*me:

capacity:=cc*pp;

plot([Re,Im](capacity),L=-4..4, color=[red,blue],
     adaptive=false, numpoints=200);

plot([Re,Im](capacity),L=0..4, color=[red,blue],
     adaptive=false, numpoints=200);

Digits:=200:
plot([Re,Im](capacity),L=0..4, color=[red,blue],
      adaptive=false, numpoints=40);
Digits:=10:

By the way, do you have an old version of Maple? If so, which one is it?

@casperyc You may wish to recompute the rank using a (weakest, fastest) Normalizer that can correctly determine that all pivot choices were not hidden zeroes. Note that we could find other examples for which `simplify` is not strong enough (either as Normalizer or as check on the Determinant).

Also note that for some other examples you may get a division-by-zero error while simplifying the determinant computed (as you did it) via the RREF. Trouble starts in pivot selection, and in the most severe case your check might produce just an error message.

And all the above is only trying to deal with identifying pivots that simplify identically to zero. This topic can become far more complicated if you also wish to account for sets of values of the indeterminate names which make the rank even smaller (because they make all pivot choices zero for some given row's pivot determinantion step). Technically, your rank=8 answer is true only under a set of caveats (on the indeterminates) that none of the pivots used become zero.

@Stephan That is worth remembering. Thanks for sharing. To summarize with a smaller example,

restart:

P := piecewise( x<=t, 1/a, x>t, (1+b)/a ):

#simplify( P*a ); simplify( P*a, piecewise );
#expand( P*a ); combine( P*a );

convert( P*a, piecewise, x );
                             /   1        x <= t
                            {                   
                             \ 1 + b      t < x 

@oldstudent There are several ways to ensure that the origin is included in the view.

For example, given your DataTable3,

plots:-display(plot([[0,0]]), plots:-pointplot( DataTable3 ));

minx,maxx:=(min,max)(DataTable3[..,1]):
plots:-pointplot( DataTable3, view=[min(-0.1,miny)..maxx, default] ); 

@oldstudent Do you want,

plots:-pointplot( DataTable3 );

That routine accepts an m-by-2 Matrix, and interprets the two columns as x and y values.

@jonlg 

• The help-page for topic rtable_indexing may be of use, if you haven't seen it yet. It ought to be linked from more pages than it is, IMO.
  
  (If I used round brackets to index into R, it was probably by omission, as I was editing your own second attachment. I try not to use round-bracket "programmer indexing" unless I need to, so that I'll know there was a reason for it when I return to the worksheet after a long time.)
  
  
• If you haven't read it yet, the Programming Manual is worth a read. (Download from here.)

@fereydoon_shekofte bytesused is not right, since it can report the bytes processed by garbage collection that may have nothing to do with the storage of the expression just assigned to a name. It is not a representation of the memory used to store anything in particular. The "used" in "bytesused" does not refer to memory that is storing a thing.

Up until recently the difference in kernelopts(bytesalloc) was a somewhat OK measure. But it can be way off, in modern versions, as garbage collection can now release memory to the OS. It was never perfect for this, either.

Part of the difficulty is that the operation to construct a thing can involve memory for temporary objects that are no longer present upon completion of the construction. It seems to me that you don't want that counted, is that right? You want only the bytes used to store a thing internally?

If you want a measure that is somewhat near proportional to the memory needed only to store a value then perhaps you could look at its `length`, or the `length` of it sprintf'd to %m format (so-called dotm). Those are not always in bytes, but they might be reasonably close to proportional. It could be tricky to handle modules.

As best I can guess at present, you might be wanting to do something like this...

MapleSimulation_test_a0.mw

@jonlg You've used `A` in two different ways. Let's keep the line that assigns to `A` the uninitialized 4-dimensional Array.

Let's focus on the first loop. You don't want to assign that operator (taking 4 parameters n,b,theta,m) to the name `A` as that will just clobber the earlier Array assignment. Let's assign that operator to the name `Afunc` instead, for now.

You may not need Afunc to be defined inside the loop. It'll depend on how this all gets unmuddled.

You can use square-bracket indexing for R, on the left hand side of that assignment in that first loop. It's makes things more clear (IMO) if you only use round-brackets for indexing when it is necessary. It's less confusing when you read this years later. I'll use you round-bracketing style below, even though I'm pretty such it mostly?! means indexing.

Now, a key thing: In your assignment to R(n+1) you have the term A(n). What is does that mean? As you had it as an operator `A` takes 4 parameters. I suspect that your A(n) here is intended as a function call to your operator `A`. But if so then what are the 2nd, 3rd, and 4th parameters (ie. b, theta, and m) supposed to be for that term A(n) in that assignment to R(n+1)?

Are you trying to set it up so that for each and every any fixed set of values of b, theta, and m the sequence of R values and A values should be computed using those formulae?

Could you be more precise than just, "generalized symbolic answer"? What kind of terms would an acceptable answer have? What do you know about omega[1] and omega[2], relative to each other?

acer

@H-R The command that does this is LinearAlgebra:-GenerateMatrix.

But it takes a fast machine about 3-4 minutesfor GenerateMatrix to process your system of 2500x2500 equations and variables. This is the wrong way to go. As I wrote before, much better would be to never form the symbolic equations at all, if you need to deal with problems of this size. You seem insistent on ignoring this key point. You should instead write code that writes numeric coefficients directly to the lhs Matrix M and the rhs Vector b, to represent M.x=b.

Your problem has more serious issues. The numeric coefficients range from 10^1000 down to -10^1000 and this will not fit inside the range of hardware double precision. Casting your data to float[8] data type will result in many entries becoming infinity, and that will not solve using LU, QR, or SVD. Moving from hardware to software working precision will make it run too long - expectedly, I could get no factorisation after even hours and hours... So, either your code which generates the system is just wrong, or your problem is very badly scaled, and your code doesn't correct t that. This is a serious problem with your code/approach.

Another issue is that the rhs of the M.x=b representation, b, seems to be almost all zero. So the system may be inconsistent, or under specified, or perhaps wrongly generated. Hard to test his serious this is until the scaling problem of entries of M is resolved.

Written well a FEM calculation producing 2500x2500/systems should be able to compute from start to end in 2-20 seconds or so on a modern machine, even within Maple. It seems that you have at least three serious problems which will either make that speed impossible or even make any solution unattainable.

@H-R In that case you should separate the unknowns from the coefficients, to produce a purely numeric Matrix M and a purely numeric Vector rhs b such that M.x=b would hold. Here x stands for the vector of unknowns in which your system is linear. For decent petformamce M must be floats not exact rationals, and preferaby have datatype=float[8].

I suspect that it would have been much better to never have formed your symbolic system explicitly in the first place. Much better would have been for your code to have simply written purely float coefficient data to data type=float[8] Matrix M and rhs b in the first place. Then pass those to LinearAlgebra:-LinearSolve. That can do a 2500x2500 linear system with a single rhs Vector in a few seconds if the data type is right.

@H-R You will not be able to solve the symbolic system for your FEM problem to get a general solution (which you would then try to instantiate at numeric values to get a solution for a particular problem). It's the wrong way to try it. A workable way is to produce a procedure which accepts numeric values for a particular problem, generates a numeric Matrix, and then solves that.

What do you intend on doing with the result? It will likely be enormous.

acer