The structures you've named are indexable.

You should probably choose according to the kind of task at hand. There is a Portal page which describes some basic features of such structures -- but it is organized by structure rather than by quality.

I'll try and illustrate with some basic qualities. It's tricky, because the whole topic is complex and convoluted; the description always seems to double back to earlier references.

Mutability: an element can be changed, by indexed assignment to that element. The `table` and `rtable` have this quality. The `list` fakes it, and indexed assignment to an entry causes the Maple kernel to simply replace the whole list with the new list, which is quite inefficient.

Uniqueness: each element is unique, and duplicates are removed (or ignored) at creation time. The `set` has this.

Growability: the `table` structure has always had this. And for the past few releases all the `rtable` flavours (Array, Matrix, Vector) do as well, on account of so-called programmer indexing [ref. 1, 2].

Indexing by non-integer values: the `table` structure allows this, where elements can be indexed by name. Eg. a `table` T can have an entry indexed as T[p] for unassigned name `p`.

If you don't need to "grow" your structure (because you know in advance how many elements you need to provide values for), and if you don't need to reassign any of the elements after initial creation of the whole structure, then a  `list` is not a bad choice. It's construction syntax is easy and terse. Maple's memory management (both creation, and garbage collection) of small lists is fast.

For very large numeric data sets the hardware datatypes (eg. float[8], integer[4]) of rtable are memory efficient. For LinearAlgebra, computational Statistics, and mapping, zipping, etc, these can be relatively much more computationally efficient. They can often be used by precompiled external functions which can be passed the hardware precision data in memory (uncopied, by location).

There are other qualities, harder to explain quickly. For example, the nature of evaluation of the structures and their contents, especially as arguments or locals of procedures. Some structures, like `table`, are of type last_name_eval.

There are some other, more rare data structures used in Maple. The `Record` can have its elements indexed by non-numeric name, but is of fixed size. The unevaluated procedure call (like a `PLOT` structure) is not indexable, and people often use `op` to pick it apart. There is a new structure called `object`, which brings a few new qualities related to OOP. Most of these are not the kind of thing you are talking about, but we can keep their existence in mind because every now and then one of them might better suit a new task.

Of course it's always possible to choose a suitable structure and then abuse it with undesirable programming practices. One of the classic bad (but common!) techniques is to build up lists by repeated concatenation such as L:=[new,op(L)] done for each new element, over and over in a loop. That is O(n^2) in memory, while a [seq(...)] construction might be O(n) complexity. And so one scales poorly, as the size `n` of the problem increases. The `if` operator can quite often allow for constuction of lists via `seq` even when element admittance is restricted by some complicated conditional.

Maple 16 has a new mechanism for data type coercion. This is great news for routines which can take an argument as either Vector[row] or mx1 Matrix or 1D Array. Or for similar situations where the structures could be converted in-place without any need to copy the data. But I am worried that allowing this feature to be applied to all routines which accept, say, Vector and `list` interchangably, just for ease of use, will lead to lots more efficient usage patterns by Maple users on account of unnecessary copying.

You have made a good start with Maple, in asking these questions.

acer

See here, for some discussion of how memory addresses of the terms in a product (of symbolic terms) can affect their ordering within that product. In turn, that can lead to differences during numerical evaluation. The same is true for addition (which I keep meaning to describe in a post...).

The 18 decimal digits in your result may be the result of conversion from HFloat, in which case the differences you notice would be only a few ulps (in HFloat).

The differences you note look much smaller than what would be dictated by the default accuracy tolerance used in your pdsolve/numeric calls. If the result is only supposed to be accurate in a much more coarse way, then why is this smaller difference interesting? If the only bits that vary are ones that arenot even close to being guaranteed, then what's the matter?

So, in terms of pdsolve,numeric per se this question does not see very interesting. It may relate to more general questions about numerics and ordering issues in Maple, but that is a topic better covered with much clearer and smaller examples (arithmetic, say) so that we can see precisely what devitates and whem.

acer

In the fourth DE there is derivative of f, which looks like it should be f(y) instead of just f.

In your Question, you wrote dsolve({sys,bcs}type=numeric) which is missing a comma just before the word `type`.

Pasting plaintext code here means using 1D Maple Notation for input instead of 2D Math input. (That's actually separate from whether you use a Worksheet or a Document).

acer

There are three common ways to get around premature evaluation of the arguments.

One way is to use a procedure which returns unevaluated if the arguments do not satisfy some type requirements, eg. all numeric, say.

Sometimes creating such a procedure can be made easier, by using the `unapply` command with its 'numeric' option. You may or may not find that helpful here.

Another way is to use unevaluation quotes, as Georgios suggested. These first two ways involve the expression-form calling sequence, where the variable ranges are each given like name=range(numeric).

A third way is to use the operator-form calling sequence, for which the various ranges of integration (or fsolving, or plotting, or Optimizing) are just range(numeric) instead of being name=range(numeric). Unfortunately this does not seem to work for the Compiled multivariate proc, unless it is wrapped in yet another proc.

f:= proc(kF1, kF2, m, k, omega, q)
      kF1+kF2+m+k+omega+q;
    end proc:

# First, using the original procedure `f`.

evalf(Int('f'(q,r,s,t,u,v),[q=0..1,r=0..2,s=0..3,t=0..4,u=0..5,v=0..6])); # Georgios
                          7560.000000

evalf(Int(f,[0..1,0..2,0..3,0..4,0..5,0..6]));
                          7560.000000

evalf(Int((q,r,s,t,u,v)->f(q,r,s,t,u,v),[0..1,0..2,0..3,0..4,0..5,0..6]));
                          7560.000000

# Now, using Compiled procedure `cf`.

cf:=Compiler:-Compile(f):

evalf(Int('cf'(q,r,s,t,u,v),[q=0..1,r=0..2,s=0..3,t=0..4,u=0..5,v=0..6])); # Georgios
                          7560.000000

CF:=unapply('cf'(q,r,s,t,u,v),[q,r,s,t,u,v],numeric):
evalf(Int(CF(q,r,s,t,u,v),[q=0..1,r=0..2,s=0..3,t=0..4,u=0..5,v=0..6]));
                          7560.000000

evalf(Int(cf,[0..1,0..2,0..3,0..4,0..5,0..6])); # a pity
Error, (in evalf/int) incompatible dimensions

# The following workaround (to the "pity" above) incurs the cost of an extra
# function call (of which there may be many, for multidimensional integration!).
evalf(Int((q,r,s,t,u,v)->cf(q,r,s,t,u,v),[0..1,0..2,0..3,0..4,0..5,0..6]));
                          7560.000000

And what's life without a few kludges to make it interesting. Somehow evalf/Int is poking at the operator-form integrand and deciding that it doesn't demand the same number of parameters as there have been supplied numeric ranges of integration. We could try and force the operator-form to work, but without having to wrap it in another (potentially costly) procedure call.

newcf := FromInert(subsindets(ToInert(eval(cf)),
                     specfunc(anything,_Inert_PARAMSEQ),
                     t->op(indets(ToInert(eval(f)),
                                  specfunc(anything,_Inert_PARAMSEQ))))):

evalf(Int(newcf,[0..1,0..2,0..3,0..4,0..5,0..6]));
                          7560.000000

acer

Of course you could try and compute the higher derivatives numerically, after calling dsolve/numeric. You could do such numerical approximation in several ways: `fdiff`, or finite differences, or applying to interpolant polynomials, etc. But I doubt that such approaches are best.

It might be better to let dsolve/numeric itself handle the whole computation instead. Apart from other considerations, that way it might strive to obtain for you approximations for all your desired outputs that satisfy the same tolerances. It can be done by augmenting your system with new targets, which happen to equal the higher derivatives which you want.

A simple example,

restart:

sys,IC := {diff(x(t),t)=x(t)*t^3}, {x(0)=1}:

sol:=dsolve( sys union IC
             union {d2x(t)=diff(x(t),t,t), D(x)(0)=0,
                    d3x(t)=diff(x(t),t,t,t), (D@@2)(x)(0)=0},
             numeric, output=listprocedure ):

X:=eval(x(t),sol):
DX:=eval(diff(x(t),t),sol):
D2X:=eval(d2x(t),sol):
D3X:=eval(d3x(t),sol):

#plot([X,DX,D2X,D3X], 0..2, view=0..10, legend=[X,DX,D2X,D3X]);

# How can we deduce those initial conditions for D(x)(0) and (D@@2)(x)(0) ?

newIC := IC union convert(eval(sys, [op(IC), t=0]),D):
eq2:=map(diff,diff(x(t),t)=x(t)*t^3,t):
newIC := newIC union convert(eval({eq2},[x(0)=1,t=0]),D);
            newIC := {x(0) = 1, D(x)(0) = 0, @@(D, 2)(x)(0) = 0}

acer

It's not clear whether you want to augment with all zeroes, or with specific data.

If you just want to create a larger Matrix whose entries in rows/columns 1..4 are the original M, then it is easy to do it in several ways. Here are two:

restart:
with(LinearAlgebra):
M:=RandomMatrix(4,4);

N:=copy(M):
N(7,7):=0:
N;

restart:
with(LinearAlgebra):
M:=RandomMatrix(4,4);

N:=Matrix(7,7,M);

Of course, you could easily modify either of those so that the result was M, and not N.

acer

The ArrayTools:-Copy command can be used to assign to the diagonal of a rectangular storage Matrix efficiently.

M:=LinearAlgebra:-RandomMatrix(3,2);

                                    [ 13  92]
                                    [       ]
                               M := [-51  59]
                                    [       ]
                                    [-65  20]

V:=LinearAlgebra:-RandomVector[row](2); # or column Vector

                                V := [-48, 0]

# Acts on newM. Use M if you want M overwritten.
newM:=copy(M):

ArrayTools:-Copy(min(3,2),V,newM,0,3+1);

newM, M;

                            [-48  92]  [ 13  92]
                            [       ]  [       ]
                            [-51   0], [-51  59]
                            [       ]  [       ]
                            [-65  20]  [-65  20]

You can also create a procedure to do it. Here's one version, with boilerplate to try and make it robust for mixed datatypes and storage orders.

replacediag := proc( A::Matrix(storage=rectangular),
                     B::Vector(storage=rectangular),
                     {inplace::truefalse:=false} )
local Anew, Bnew, m, n;
  if inplace then Anew := A; else Anew := copy(A); end if;
  if rtable_options(B,':-datatype')<>rtable_options(A,':-datatype') then
    Bnew := Vector(B,':-datatype'=rtable_options(A,':-datatype'));
  else Bnew := B;
  end if;
  (m,n) := LinearAlgebra:-Dimensions(Anew);
  ArrayTools:-Copy(min(m,n), Bnew, Anew, 0,
    `if`(rtable_options(Anew,':-order')=':-Fortran_order', m+1, n+1));
  Anew;
end proc:


replacediag(M,V), M; # M is not overwritten. Result is a new Matrix.

                            [-48  92]  [ 13  92]
                            [       ]  [       ]
                            [-51   0], [-51  59]
                            [       ]  [       ]
                            [-65  20]  [-65  20]

replacediag(M,V,inplace), M; # M is overwritten. This saves time/memory.

                            [-48  92]  [-48  92]
                            [       ]  [       ]
                            [-51   0], [-51   0]
                            [       ]  [       ]
                            [-65  20]  [-65  20]

acer

I did something like that here, embedding "plot" axes onto an image.

The basic idea was to produce an empty plot containing just the axes (or in your case, a textplot without axes perhaps), export that as an image, read in that image using ImageTools, and finally overlay it onto your other (background) image. The overlaying is done just by overwriting the Array entries -- copying the overlay image's Array's nonzero entries onto the background image's Array.

There was some trickiness with adjusting for whitespace border around the exported plot. The referenced link might help with that.

acer

I changed val_M0 to vec_M0 (which you probably had in your original, a transcription error, perhaps).

That allowed the posted code to produce the nonreal result,

fsolve({mat_M_V[2,1]=0,mat_M_V[1,2]=0,Im(mat_M_V[1,1])=0,Im(mat_M_V[2,2])=0},
       {alpha11=0,alpha22=0,alpha12=0,theta=0},'fulldigits');

           {theta = -1.666216272 + 0.4818505195 I, alpha11 = -1.570796327, 

            alpha12 = -9.424777961, alpha22 = 3.141592654}

So then I specified a real range for theta, and got (with or without option 'fulldigits'),

fsolve({mat_M_V[2,1]=0,mat_M_V[1,2]=0,Im(mat_M_V[1,1])=0,Im(mat_M_V[2,2])=0},
       {alpha11=0,alpha12=0,alpha22=0,theta=0},theta=-infinity..infinity);

             {theta = 3.997794044, alpha11 = -1.570796327, 

              alpha12 = 1.570796327, alpha22 = -1.570796327}

Interestingly, 1.570796327 = Pi/2, which makes me wonder whether some exact approach is also possible (before assigning float values to parameters).

acer

This is an important question. It's actually one of the most basic instances of a more general set of computational tasks, for which dsolve/numeric has been augmented with a whole slew of facilities. The functionality is known as Events.

The basic notion is like this: a rootfinder (like fsolve, etc) operating on the result from a call to dsolve/numeric is acting a bit in the dark. It needs to be told special information about the dependent functions (like x(t) or y(t) say) in order to get the best chance of finding the desired root. The seperate rootfinder doesn't even know that the solution exists (which can be helpful knowlegde). It doesn't know the range, or the stiffness, etc, that characterizes the curve.  But dsolve/numeric itself is computing the solved approximations to the dependent functions, so it is in the best position to be able to invert that computation.

In practice, using dsolve/numeric with "events" can get such inversion results using less resources (memory and time), and quite often more accurately.

The accuracy issue is subtle. The accuracy of the dsolve/numeric solution is (generally) only as good as what is specified by its abserr and relerr tolerances. Suppose that you obtain a procedure from dsolve/numeric, which computes y(t) for any given value of `t`. Bumping up Digits, so as to try the separate rootfinding (eg. fsolve) to obtain a high accuracy value of the independent variable `t` that satisfies y(t)=K for a given numeric value K, is not going to get you a more accurate answer. The rootfinding results are only as good as the interpolating spline that dsolve/numeric uses to compute y(t) in the first place.

So (I feel that) there is some risk in using a separate rootfinder like fsolve, of being beguiled into believing that one has control over the accuracy of the computed roots.

It's a little more work to set up, for involved examples, but it can pay off for larger and more difficult or computationally intensive projects.

Here is an illustration, using Markiyan's example. I've also used output=listprocedure, because that makes subsequent computation of x(t) and y(t) easier (without needing `op`).

restart:

sys := (D(x))(t) = y(t)-x(t), (D(y))(t) = x(t)*(D(x))(t)-2*y(t):
ic := x(0) = 1, y(0) = -1:

# Default for nonstiff rfk45 is abserr=1e-7, relerr=1e-6
sol := dsolve({ic, sys}, numeric
              , parameters=[Py]
              , events=[[y(t)-Py,halt]]
              , output=listprocedure):

solt:=eval(t,sol):
solx:=eval(x(t),sol):
soly:=eval(y(t),sol):

interface(warnlevel=0):

solt(parameters=[-0.5]):
ans := CodeTools:-Usage( solt(1000) );

memory used=213.85KiB, alloc change=0 bytes, cpu time=15.00ms, real time=12.00ms
                          ans := 0.623905007259310

soly(ans);

                             -0.500000000000000

restart:

sys := (D(x))(t) = y(t)-x(t), (D(y))(t) = x(t)*(D(x))(t)-2*y(t):
ic := x(0) = 1, y(0) = -1:

# Default for nonstiff rfk45 is abserr=1e-7, relerr=1e-6
sol := dsolve({ic, sys}, numeric
              , parameters=[Px,Py]
              , events=[[x(t)-Px,halt],[y(t)-Py,halt]]
              , output=listprocedure):

solt:=eval(t,sol):
solx:=eval(x(t),sol):
soly:=eval(y(t),sol):

interface(warnlevel=0):

solt(parameters=[1e10, -0.5]):

ans := CodeTools:-Usage( solt(1000) );

memory used=216.88KiB, alloc change=127.98KiB, cpu time=0ns, real time=12.00ms
                          ans := 0.623905007259310

soly(ans);

                             -0.500000000000000

solt(parameters=[0.95, 1e10]):
ans := CodeTools:-Usage( solt(1000) );

memory used=31.67KiB, alloc change=0 bytes, cpu time=0ns, real time=2.00ms
                          ans := 0.0253259108029954

solx(ans);

                              0.950000000000000

Here's a worksheet, comparing using Events and fsolve for this example. Note both memory allocation as well as speed differences. This is just expository and, being based on a single example, is not supposed to prove anything. For one's own serious projects, it'd be best to try out the alternatives.

dsolveinv0.mw

acer

Have you tried CodeGeneration[Matlab] ? That will work on some Maple procedures or expressions (but not on entire Documents or Worksteets).

acer

I have seen it done on the same host, using the commandline interface. That is, starting the kernel and the client (interface) separately.

I don't know whether it could be done across a network. I'll give it a try (later, sorry).

By "X-window interface" I presume that you mean the Classic GUI. If anything, I would expect that the Standard (Java) GUI would be more difficult than Classic here, as far as accepting parameters when launched and not automatically firing up its own kernel.

acer

I found an ubuntu 12.04 machine and reproduced the problem in the commandline interface. I will look harder, and submit a bug report. You don't need to fiddle with libname and other stuff that I suggested earlier (now deleted).

%maple15.01 -s
    |\^/|     Maple 15 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2011
 \  MAPLE  /  All rights reserved. Maple is a trademark of
   Waterloo Maple Inc.
      |       Type ? for help.

> kernelopts(version);           
            Maple 15.01, X86 64 LINUX, Jun 1 2011, Build ID 635833

> ssystem(`more /etc/*release*`);
[0, "DISTRIB_ID=Ubuntu
    DISTRIB_RELEASE=12.04
    DISTRIB_CODENAME=precise
    DISTRIB_DESCRIPTION="Ubuntu 12.04 LTS"

    "]

> int(x^n, x=2..3);
                                     3
                                    /
                                   |    n
                                   |   x  dx
                                   |
                                  /
                                    2

I strongly suspect that the underlyng problem is another issue which has been previously reported. On that same system,

> log[2](8);
Error, (in iroot) powering may produce overflow

It's quite plausible that `int` is generating that error internally, but catching and suppressing it, with the end result being that it returns unevaluated.

Here is a comment to a previous report, suggesting that you contact Maplesoft's Tech Support ( support@maplesoft.com )

Also, the example succeeds in 64bit Maple 16.00 on that same machine.

acer

M:=LinearAlgebra:-RandomMatrix(5,generator=0..2);

                                 [0  0  1  1  1]
                                 [             ]
                                 [2  0  1  0  0]
                                 [             ]
                            M := [2  1  2  2  0]
                                 [             ]
                                 [1  0  2  0  2]
                                 [             ]
                                 [0  1  0  0  1]

rtable_num_elems(M,'NonZero');

                                     14

add(`if`(x<>0,1,0),x=M);

                                     14

You can subtract the above numbers from the total number of entries (the product of dimensions, 5*5 or 18*18, etc) to get the number of zeros.

acer

This is not short and pretty,

> expr := -2*Pi*sin(Pi*a)/(-1+cos(2*Pi*a)):

> combine(simplify(convert(convert(expand(expr),tan),sincos)));

                                      Pi
                                   ---------
                                   sin(Pi a)

acer