@PatrickT 

You won't see the same order of speedup by giving Cru option remember as you might for xt and yt, unless you plan on calling Cru itself in several circumstances. But the methodology is pretty similar.

It doesn't arise here, but in general you might have to be careful about scoping, sometimes, and ensure that the names used in the unapply and funciton calls (eg. t) aren't assigned.

One interesting subtlety is that, if your original Cru is a procedure which does not apply evalf to its final result and if it contains exact numerics like sqrt(3), then you might want a resulting procedure which does in fact apply evalf. This doesn't affect plotting, as `plot` itself applies evalf or evalhf to final results. But you can see the distinction below between ct(0.7) and anotherct(0.7).

restart:

SOL:=dsolve({diff(x(t),t)=y(t)+x(t),diff(y(t),t)=y(t)-x(t),x(0)=1,y(0)=1},
             numeric,output=listprocedure):
xt := subsop(3=remember,eval( x(t), SOL )):
yt := subsop(3=remember,eval( y(t), SOL )):

# Let's suppose your original Cru is a procedure that takes two arguments (a,b)

Cru := (a,b) -> sqrt(3)/3*sqrt(1/(a+1/10))+sqrt(3)/90*sqrt(1/(b+1/10)^3):

ct := unapply(subs(A=xt(t),B=yt(t), Cru(A,B)),
              t, 'proc_options'=['remember','operator','arrow']):

lprint(eval(ct)); # note absence of evalf

ct(0.7);

temp := subs(unassigneddummy=Cru(xt(t),yt(t)), t->evalf(unassigneddummy)):

anotherct := subsop(3=op({op(3,eval(temp)),'remember'}), eval(temp) ):

lprint(eval(anotherct)); # note presence of evalf

anotherct(0.7);

plot(ct, 0.0 .. 0.7);
plot(anotherct, 0.0 .. 0.7);

@PatrickT You can place `option remember` on the procedures returned from dsolve/numeric. (Or you could wrap them in procedures that have option remember.)

Compare both times taken to produce the second plot, P2, below.

restart:
SOL:=dsolve({diff(x(t),t)=y(t)+x(t),diff(y(t),t)=y(t)-x(t),x(0)=1,y(0)=1},
             numeric,output=listprocedure):
xt := eval( x(t), SOL ):
yt := eval( y(t), SOL ):

f := (a,b)->sin(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P1:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=red ):
time()-st;

                                   0.080

f := (a,b)->cos(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P2:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=blue ):
time()-st;

                                   0.070

plots:-display(P1,P2);

restart:
SOL:=dsolve({diff(x(t),t)=y(t)+x(t),diff(y(t),t)=y(t)-x(t),x(0)=1,y(0)=1},
             numeric,output=listprocedure):
xt := subsop(3=remember,eval( x(t), SOL )):
yt := subsop(3=remember,eval( y(t), SOL )):

f := (a,b)->sin(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P1:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=red ):
time()-st;

                                   0.070

f := (a,b)->cos(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P2:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=blue ):
time()-st;

                                   0.010

plots:-display(P1,P2);

You can also force your plots to use specific domain points (ie. pointplot instead of odeplot) and thus enforce re-use of the same `t` values (and thus benefit from option remember) in subsequent plots of different `f` compositions.

Array output directly obtained from dsolve/numeric is a convenience. But you can also produce it yourself, after calling dsolve. So it might just be an inconvenience that dsolve & parameters doesn't allow ouput=Array, which you can work around.

This page, now spammed (June 08, 2011), deserves a screengrab.

acer

@PatrickT It would be interesting to learn if you see performance improvements in your own examples. I did not investigate whether any relative performance gain was affected by how many distinct sets of ICs are supplied.

I'd like to turn (or see turned) the animation version into a procedure. I suspect that it might have to do things in this order: generate and save (to a table or whatever) all the odeplots, inspect the PLOT structure from each odeplot and take the max domain values from all of those, use those values for the xrange[i] and call fieldplot, then finally form allframes the call to `display` using the seq of odeplots and the fieldplot.

A good procedure would be able to accept (possibly separately) any valid options for the fieldplot as well as for the odeplots and/or the final display.

Sorry that I have only time to follow right now.

Let's suppose that you had a type-check in your code, like,

  type(f, procedure)

and that this worked ok for various inputs `f`, including `sqrt`. You could try changing that to, say,

  type(f, {procedure,And(`module`,'appliable')})

acer

@Andreas Madsen In your `residualsumofsquares` procedure, try,

  add(evalf((ln(Y[i])-ln(eval(f, T = X[i])))^2), i = 1 .. N);

@Andreas Madsen In your `residualsumofsquares` procedure, try,

  add(evalf((ln(Y[i])-ln(eval(f, T = X[i])))^2), i = 1 .. N);

@Christopher2222 It looked to me like you might have been trying to just type in all of it ( the d^2/dx^2 bit). I could be wrong, but I believe you can't just type that in letter by letter and have it be parsed a the differential operator. Ie, you might have to use command-completion or the palette entry (and raise the "d" on top to be d^2, to match the dx^2 on bottom). That way the parser knows how to accept it. It's similar for just d/dx, where typing that in letter by letter is not sufficient.

@Christopher2222 It looked to me like you might have been trying to just type in all of it ( the d^2/dx^2 bit). I could be wrong, but I believe you can't just type that in letter by letter and have it be parsed a the differential operator. Ie, you might have to use command-completion or the palette entry (and raise the "d" on top to be d^2, to match the dx^2 on bottom). That way the parser knows how to accept it. It's similar for just d/dx, where typing that in letter by letter is not sufficient.

What you did wrong was misinterpret a question about products as if it were a question about sums.

You wrote, "I have a maple exercise where it ask me to find the limit of a sum...".

But the picture is that of product, using capital Greek Pi instead of capital Greek Sigma. That's why those are the graphic symbols used for products and sums, because they are a kind of abbreviation ("P" as in Pi for product, and "S" as in Sigma for sum).

> f:= 2*k/(2*k+11);

                              2 k   
                            --------
                            2 k + 11

> g:= sum(f,k=1..n);

            6508   11    /    13\   11                 
        n + ---- - -- Psi|n + --| - -- gamma - 11 ln(2)
            315    2     \    2 /   2            
      
> limit(n^(11/2)*g,n=infinity);

                            infinity

> g:= product(f,k=1..n);

                                        (1/2)
                   10395 GAMMA(n + 1) Pi     
                   --------------------------
                                /    13\     
                        64 GAMMA|n + --|     
                                \    2 /     

> limit(n^(11/2)*g,n=infinity);

                         10395   (1/2)
                         ----- Pi     
                          64          

Or, using the inert forms of the two commands, to show the graphic symbols for product and summation,

> restart:

> f:= 2*k/(2*k+11);

                              2 k   
                            --------
                            2 k + 11

> g:= Sum(f,k=1..n);

                           n           
                         -----         
                          \            
                           )     2 k   
                          /    --------
                         ----- 2 k + 11
                         k = 1         

> limit(n^(11/2)*value(g),n=infinity);

                            infinity

> g:= Product(f,k=1..n);

                           n             
                       ,--------'        
                          |  |           
                          |  |     2 k   
                          |  |   --------
                          |  |   2 k + 11
                         k = 1           

> limit(n^(11/2)*value(g),n=infinity);

                         10395   (1/2)
                         ----- Pi     
                          64          

acer

What you did wrong was misinterpret a question about products as if it were a question about sums.

You wrote, "I have a maple exercise where it ask me to find the limit of a sum...".

But the picture is that of product, using capital Greek Pi instead of capital Greek Sigma. That's why those are the graphic symbols used for products and sums, because they are a kind of abbreviation ("P" as in Pi for product, and "S" as in Sigma for sum).

> f:= 2*k/(2*k+11);

                              2 k   
                            --------
                            2 k + 11

> g:= sum(f,k=1..n);

            6508   11    /    13\   11                 
        n + ---- - -- Psi|n + --| - -- gamma - 11 ln(2)
            315    2     \    2 /   2            
      
> limit(n^(11/2)*g,n=infinity);

                            infinity

> g:= product(f,k=1..n);

                                        (1/2)
                   10395 GAMMA(n + 1) Pi     
                   --------------------------
                                /    13\     
                        64 GAMMA|n + --|     
                                \    2 /     

> limit(n^(11/2)*g,n=infinity);

                         10395   (1/2)
                         ----- Pi     
                          64          

Or, using the inert forms of the two commands, to show the graphic symbols for product and summation,

> restart:

> f:= 2*k/(2*k+11);

                              2 k   
                            --------
                            2 k + 11

> g:= Sum(f,k=1..n);

                           n           
                         -----         
                          \            
                           )     2 k   
                          /    --------
                         ----- 2 k + 11
                         k = 1         

> limit(n^(11/2)*value(g),n=infinity);

                            infinity

> g:= Product(f,k=1..n);

                           n             
                       ,--------'        
                          |  |           
                          |  |     2 k   
                          |  |   --------
                          |  |   2 k + 11
                         k = 1           

> limit(n^(11/2)*value(g),n=infinity);

                         10395   (1/2)
                         ----- Pi     
                          64          

acer

@Alejandro Jakubi I think that using 2D Math input for summation can be OK, as long as one knows what one is getting. This can be achieved by deliberately typing it in. There's already enough subtlety in the distinctions amongst the `sum`, `Sum`, and `add` commands that any more difficulty from the palette makes it all too much for the inexperienced user.

The main palette problem is that the Expression palette only has an entry for the `sum` command. We've see more than enough posts on this forum from members using `sum` to do summation with a literal (non-symbolic) finite upper value for the index, and then running amok of the lack of special evaluation rules for `sum`.

So I'd suggest that anyone wanting 2D Math input for summation do it by typing, as a reminder of what one is asking for. Type in the first three characters only (and not yet the left-bracket), and then hit the command-completion keystrokes (Ctl-space on MS-Windows). That way one can get the 2D input form of inert `Sum`, or active `sum` and be aware (at time of authoring, at least) of what is being requested.

What else does 2D Math need, to not hammer the new user so hard on this topic? Allowing for 2D Math representation of calls to `add` would be nice. But in order for this to not create more confusion it would have to be accompanied by a hover-over/balloon/tool-tip for 2D summation symbols already appearing in a worksheet, so that one could tell at a glance what underlies the graphic. Having to convert 2D Math input to 1D and back again using context-menu actions, just to discover what is being represented, is very user-unfriendly. This would also help with the inert/active forms, since the distinction between black and gray summation graphics is visually weak and not generally known.

There will be people who want to do finite non-symbolic summation, and use 2D Math input for it, and not know/want to handle premature evaluation due to the lack of special evaluation rules on `sum`.

@Alejandro Jakubi I think that using 2D Math input for summation can be OK, as long as one knows what one is getting. This can be achieved by deliberately typing it in. There's already enough subtlety in the distinctions amongst the `sum`, `Sum`, and `add` commands that any more difficulty from the palette makes it all too much for the inexperienced user.

The main palette problem is that the Expression palette only has an entry for the `sum` command. We've see more than enough posts on this forum from members using `sum` to do summation with a literal (non-symbolic) finite upper value for the index, and then running amok of the lack of special evaluation rules for `sum`.

So I'd suggest that anyone wanting 2D Math input for summation do it by typing, as a reminder of what one is asking for. Type in the first three characters only (and not yet the left-bracket), and then hit the command-completion keystrokes (Ctl-space on MS-Windows). That way one can get the 2D input form of inert `Sum`, or active `sum` and be aware (at time of authoring, at least) of what is being requested.

What else does 2D Math need, to not hammer the new user so hard on this topic? Allowing for 2D Math representation of calls to `add` would be nice. But in order for this to not create more confusion it would have to be accompanied by a hover-over/balloon/tool-tip for 2D summation symbols already appearing in a worksheet, so that one could tell at a glance what underlies the graphic. Having to convert 2D Math input to 1D and back again using context-menu actions, just to discover what is being represented, is very user-unfriendly. This would also help with the inert/active forms, since the distinction between black and gray summation graphics is visually weak and not generally known.

There will be people who want to do finite non-symbolic summation, and use 2D Math input for it, and not know/want to handle premature evaluation due to the lack of special evaluation rules on `sum`.

@Christopher2222 Your new version which uses StringTools:-Parsetime makes the very common mistake of building up a list (or set) by repeated concatenation. This results in a collection of increasing length lists being created, and having to be subsequently memory-managed and garbage collected.

Think of the collection of longer and longer instances of list `a` as if they were all stacked and arranged in a triangle. That has area, which is to say that this implementation incurs a quadratic cost (of resources). Fortunately, there are some common remedies for this common malady. Often, `seq` or `map` can get used instead, to form just a single list, and they can often reinstate (near) linear cost. Maybe see here, and for fun here too.

restart:

newdaytoday := proc (d1, m1, y1, d2, m2, y2)
local i, j, k, a, startinterval, endinterval;  
uses StringTools:
  a := []:
  for k from y1 to y2 do
    for j from 1 to 12 do
      for i from 1 to 31 do
        if ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-weekDay = -1 then break: end if:
        a := [op(a), cat("", ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-monthDay,
                         "/", ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-month,
                         "/", ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-year)]:
      end do:
    end do:
  end do:
  startinterval := ParseTime("%Y-%m-%d", cat(y1, "-", m1, "-", d1)):-yearDay:
  endinterval := ParseTime("%Y-%m-%d", cat(y2, "-12-31")):-yearDay
                 - ParseTime("%Y-%m-%d", cat(y2, "-", m2, "-", d2)):-yearDay+1:
  a[startinterval .. -endinterval];
end proc:

newerdaytoday := proc (d1, m1, y1, d2, m2, y2)
local i, j, k, a, startinterval, endinterval;  
uses StringTools:
  a := [seq(seq(seq(`if`(ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-weekDay = -1, NULL,
                    cat("", ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-monthDay,
                        "/", ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-month,
                        "/", ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i)):-year)),
                    i=1..31), j=1..12), k=y1..y2)];
  startinterval := ParseTime("%Y-%m-%d", cat(y1, "-", m1, "-", d1)):-yearDay:
  endinterval := ParseTime("%Y-%m-%d", cat(y2, "-12-31")):-yearDay
                 - ParseTime("%Y-%m-%d", cat(y2, "-", m2, "-", d2)):-yearDay+1:
  a[startinterval .. -endinterval];
end proc:

S1:=newdaytoday(25,12,1900,15,3,1950):
S2:=newerdaytoday(25,12,1900,15,3,1950):
evalb(S1=S2);

plots:-display(
   plot([seq([i,time(newdaytoday(25,12,1900,15,3,1900+i^2))],i=1..15)],color=red),
   plot([seq([i,time(newerdaytoday(25,12,1900,15,3,1900+i^2))],i=1..15)],color=green));

Having a better order of complexity is more important (maybe only eventually...) than a constant factor speedup.In your revision you could also optimize multiple instances of the very same ParseTime invocation by doing it just once, assigning its value, and then reusing that assigned varaiable. Ie, inside you inner loop,

        T:=ParseTime("%Y-%m-%d", cat(k, "-", j, "-", i));
        if T:-weekDay = -1 then break: end if:
        a := [op(a), cat("", T:-monthDay,
                         "/", T:-month,
                         "/", T:-year)]:

But for a duration of 200 days, that code change only brings a constant factor 15% speedup. The overall performance still grows quadratically, and that is what is characterizing the performance here.

acer

@DuncanA