If you are calling dsolve,numeric inside a loop or process with many iterations, then memory allocation can be a problem. If that is the difficulty you're seeing, and if your ODEs are IVPs which differ only by some parameters in the formulae then you may be able to pull the dsolve call outside the loop, using its parameters option.

There is a simple example of the difference in performance and resources here. But I've seen a GlobalOptimization example where using the technique brought memory allocation down from many GB to a mere 40MB.

acer

Your expressions appear to have inert Int calls in them. Try value or evalf, not eval.

 (**) expr:=Int(x*M,x=0..1/(3*M));
                                         1
                                        ---
                                        3 M
                                       /
                                      |
                             expr :=  |     x M dx
                                      |
                                     /
                                       0

 (**) eval(expr,M=1);
                                     1/3
                                    /
                                   |
                                   |     x dx
                                   |
                                  /
                                    0

 (**) value(subs(M=1,expr));
                                     1/18

 (**) evalf(subs(M=1,expr));
                                 0.05555555556

acer

I'm sure that someone will provide a good response to your "area" query. I'd like to make a few side comments.

The code you posted for routine MyMandel looks similar to what is in one of the "code edit regions"of the FractalFun worksheet available here on the Application Centre.

That routine is (relatively) very slow when run using Maple's normal interpreter (which is what that worksheet does, using LinearAlgebra:-Map to populate a Matrix to be used to form an image). The whole construction would run considerably faster if edited to utilize Maple's evalhf interpreter to use MyMandel as initializer in a Matrix constructor call. And it would run faster still if it used Maple's Compiler:-Compile facility on the MyMandel routine. Much faster still would be a Compiled version of MyMandel which accepted a Matrix as a new argument (to be acted upon inplace).

It's surprising thing that the version of the demo up on the Application Centre is so slow, and doesn't make use of  Compile. This is precisely the sort of numeric routine for which the Compiler seems designed. (Maple on 32bit MS-Windows is bundled with the Watcom C compiler, which Compile uses. Every other platform can use a free compiler, such as gcc. So availability is not an issue.) And the demo doesn't even benefit from evalhf, which has been built into Maple versions for over 10 years and is a standard way to accelerate one's floating-point Maple routines.

Sure, one might alternatively use Maple's new Threads package to populate the Matrix (image data) for this embarrassingly parallelizable task. But that would get only a few-fold speedup (linear w.r.t. the number of cores, if one were lucky...). Compiling the whole thing, on the other hand, brings 100s of time speedup immediately, with little extra effort. I suspect that this is the sort of reasoning behind Alec's comment here.

Here is a version one might compare with that in the demo worksheet.

restart:

MyMandel := proc(N::integer[4], X::Vector(datatype=float[8]),
                 Y::Vector(datatype=float[8]),
                 M::Array(datatype=float[8]))
local ii::integer[4],jj::integer[4],i::integer[4],new::integer[4],
      zvalr::float[8],zvali::float[8],newzvalr::float[8],
      cvalr::float[8],cvali::float[8];
  for ii from 1 to N do
    for jj from 1 to N do
      zvalr:=X[jj]; zvali:=Y[ii];
      for i from 2 to 25 do
        newzvalr:=zvalr^2-zvali^2+X[jj];
        zvali:=2*zvalr*zvali+Y[ii];
        zvalr:=newzvalr;
        if zvalr^2+zvali^2 > 4 then
          M[ii,jj]:=i; break;
        end if:
      end do:
    end do;
  end do;
  NULL;
end proc:

t1, ba := time(), kernelopts(bytesalloc):

pts := 400:
(X1, X2, Y1, Y2) :=-2.0, 0.7, -1.35, 1.35:

X := Vector(pts, i->X1+(X2-X1)*(i-1)/(pts-1), datatype=float[4]):
X := Vector(X,datatype=float[8]):
Y := Vector(pts, i->(Y1+(Y2-Y1)*(i-1)/(pts-1)), datatype=float[8]):

pMyMandel := Compiler:-Compile(MyMandel):
Val := Array(1..pts, 1..pts, datatype=float[8]):
pMyMandel(pts,X,Y,Val):

time()-t1,kernelopts(bytesalloc)-ba;

img:=ImageTools:-Create(Val,fit=true):
ImageTools:-View(img);

For this particular fractal image's creation, another factor of two speedup could be attained by a cheap symmetry trick. But it's hardly necessary, I think. If anyone wants to post hard timing numbers, to compare with the original...

acer

Is this of some use? (You didn't mention whether variables K, s, i1, i2 and i3 were to be part of the solution or not. It may be that those cannot all be eliminated.)

> eqns:=[Vo = K*(Vb-Vo), i3 = C3*Vb*s, Va-Vb = R3*i3,
>        i2 = C2*s*(Vo-Va), Vi-Va = R1*i1, i3 = i1+i2]:

> Vo=eval(Vo,eliminate(eqns,{Vo,Vb,Va,i1,i2,i3})[1]);

Vo = K Vi/(1 + s R3 C3 + R1 C3 s + K + s R3 C3 K + K R1 C3 s + C2 R1 s
 
           2                2
     + C2 s  R3 C3 R1 + C2 s  R3 C3 K R1)

acer

Aside from other answers, I would suggest that it's generally not good practice to use a protected name such as `sum` as a local variable.

For fun, compare with,

f := proc(a)
local s;
  s := a mod 9;
  if s=0 then 9 else s; end if;
end proc:

acer

It's Matrix(), not MATRIX(), to construct them.

The "end if" and "end do" statements are all mixed up and interleaved. Those have to be sorted into a logically valid layout.

Loops like this simply repeat the exact same command each time through, since the left hand side of the assignment, X_0, always refers to the same thing. It seems very likely that's not what you intend.

for i from 1 to n do
  X_0:=0
end do:

acer

plot(sin(x),labels=["my x-axis","my y-axis"],
     labeldirections=[horizontal,vertical]);

acer

If you are using 2D Math input mode (the default, so you likely are) then the "OverBarUnder" item in the Layout palette can do this.

Look for the Layout palette, as a panel at the left of the Standard Maple interface (GUI). Open it. One of the items appears as an A with a bar under it. In 2D Math mode, simply double-click on that. Doing so should insert the appropriate doubly-underlined template, for which you can overwrite the content A.

acer

The string "i" doesn't change just because you changed the value of i. It'll always be simply "i". That's no bug.

You can fix your code in several ways. In this first way below, your procedure test is unchanged while the way to call it using i is altered.

> restart:
> test:=proc(s::string)
> print( cat("Matrix_",s) );
> end proc:

> for i from 3 to 4 do
> test(convert(i,string));
> od;
                                  "Matrix_3"
 
                                  "Matrix_4"

In this next way, both the procedure `test` and the way it gets called are altered.

> restart:
> test:=proc(s)
> print( cat("Matrix_",s) );
> end proc:
>
> for i from 3 to 4 do
> test(i);
> od;
                                  "Matrix_3"
 
                                  "Matrix_4"

Notice that,in that last example, `cat` is automatically converting the integer s to a string so that it can concatenate with it. But that wouldn't work if s were 9.3 a float. So you could bolster it like so,

> restart:
> test:=proc(s)
> print( cat("Matrix_",convert(s,string)) );
> end proc:

> for i from 3 to 4 do
> test(i);
> od;
                                  "Matrix_3"
 
                                  "Matrix_4"
 
> test(9.3);
                                 "Matrix_9.3"

acer

Replace pi with Pi.

In Maple, lowercase pi doesn't mean anything special.

acer

Your integration range is being supplied as f=100..1000 where the range endpoints 100 and 1000 are exact integer (not floating-point) values. In that case, calling value() will bring about exact integration rather than numerical quadrature.

To get (here, much faster) numerical quadrature call evalf() instead of value(), or replace 100..1000 with 100.0..1000.0 as a floating-point range.

See the ?evalf,Int help-page.

acer

This has the advantage that it doesn't involve having to repeat actions with the mouse pointer (ie. is programmatic). And it is sometimes easier than using frontend (which can sometimes be trickier, when it freezes too much).

As mentioned elsewhere,

gdiff := (f,x)->thaw(value(subs(x=freeze(x),Diff(f,x)))): #


gdiff(5*x^2+y,x^2);

(The pre tag is acting up! It looks good in the editor, but is lost on submission.)

acer

First, you should lose the habit of loading a package using `with`, and then using its exports inside a proc.

Instead, either use the full name inside the proc, or use `uses` or `use.

For example you could put the line,

uses LinearAlgebra:-Modular;

at the beginning of your proc, and remove that `with` call altogether.

Next, avoid a mismatch in datatypes. Try changing the last argument of `Create` from ('integer')[] to simply 'integer'. Alternatively, insert datatype=integer[kernelopts(wordsize)/8] as a last argument to the Matrix call where TestMatrix is created.

acer

It looks as if you may be getting bitten by "premature evaluation".

You pass X(a,b) to Optimization:-Maximize. But, as is Maple's usual behaviour, that argument X(a,b) gets evaluated up front before a or b are assigned numeric values. Your routine X does not appear to be prepared for that symbolic input a,b. The unintended evaluation of X(a,b) prior to parameters a and b getting values (which happens only subsequently, inside the routine X), is premature.

There are a few ways to handle that. One is to have X return unevaluated if a and b are not both numeric.

X := proc(a,b)
local ...
  if not( type(a,numeric) and type(b,numeric) ) then
    return 'procname'(args);
  end if;
  ...
end proc:

Another way is to pass arguments to Maximize in so-called operator form. (See the ?Optimization,General,OperatorForm help-page, and in particular the last Example.) For that form, you'd be passing just X rather than X(a,b), and [.6,.2] rather than [a=.6,b=.2], etc. You'd also need to read the section on passing the constraints as a list of procedure(s).

Another way that might work better is to pass 'X'(a,b) so as to delay evaluation.

This also comes up with plot, fsolve, evalf/Int, etc.

acer

> fsolve(unapply(hite(t)-11853,t),0..1000);
                               33.00303938

> fsolve(t->hite(t)-11853);
                               33.00303938

In Maple 13, you could call dsolve/numeric using 'ground' as a 'parameter'. See ?dsolve,numeric,IVP for explanation of the parameters option.

You might also have a look at the ?dsolve,numeric,Events help-page. That might be what you'd call the "elegant" way to do this sort of thing.

acer