@Alec Mihailovs One of my example's points is that your use of evalhf not relevant to how the result is obtained. The `eval` call within the `evalhf` call that you made in Maple 15 does not inherit the Digits setting (of evalhf).

As I mentioned, a result from `eval` similar to that using hardware doubles may be due to the NLPSolve result being an HFloat (and cascading of that through the arithmetic of the `eval`). But it's not because it is called from within evalhf.

For some subtle reasons (to do with a small amount of corner case examples, and so subtle that I forget them over and over again, and always have to comb email archives to re-find them...) the indexed form of `evalhf` is slightly superior to the two-argument calling sequence. I believe that's why the help page for evalf lists the former and not the latter in recent versions. It's understood, you have your own preference, and I for one will not criticize it.

[edited: ..had to dig, for this] Here's one such corner case. Some users, after getting in the habit of using the non-indexed form like evalf(expression, n) , get confused why these two calls below don't work the same way. It's not that evalf(..., n) is wrong, but more that people who aren't in the habit of using it won't run into this confusing example. (So, it may not affect an expert like yourself, but there's still a reasonable argument for suggesting one over the other to less experienced users.)

> evalf(1/3, 7);

                           0.3333333

> g:=1/3, 7;

                              1   
                              -, 7
                              3   

> evalf(g);

                        0.3333333333, 7.

I'm not a big fan of how votes change the order of posts here, either. Together with the mix and mix-ups of answers vs comments, it makes a mess of lots of threads.

@Alec Mihailovs I do not wish to give offense, but it looks to me as if the stated problem was with the precision of evaluation of the objective at the returned solution (as opposed to say the question of whether that solution were globally or locally optimal, etc).

In Maple 10, which the submitter mentioned using, NLPSolve does not return a result with HFloats. So not even double precision will be used under `eval`. Also, `eval` cannot be called from within `evalhf` in Maple 10.

In Maple 11, it still wouldn't be adequate, if greater working precision that the current Digits setting (eg, default of 10) were wanted:

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                               [0., [x = 3.]]

> evalf(subs(sol[2],1/x));

                                0.3333333333

> evalhf(eval(1/x,sol[2]));

                            0.333333333299999979

So, for the stated issue, in earlier Maple it's probably better to either raise Digits first. Or invoke evalhf with its indexed form, like evalf[n](...).

In Maple 15, double precision may get used, due to the presence of HFloats in the NLPSolve result.

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                    [0., [x = HFloat(3.0)]]

> evalf(subs(sol[2],1/x));

                   HFloat(0.3333333333333333)

> evalhf(eval(1/x,sol[2]));
                      0.333333333333333315

> lprint(sol);
[0., [x = HFloat(3.)]]

> kernelopts(version);

   Maple 15.00, X86 64 WINDOWS, Mar 21 2011, Build ID 582305

But if even greater working precision is required, then raising Digits above 15, or calling evalf[n](..) with n>15, may be required.

@Alec Mihailovs I do not wish to give offense, but it looks to me as if the stated problem was with the precision of evaluation of the objective at the returned solution (as opposed to say the question of whether that solution were globally or locally optimal, etc).

In Maple 10, which the submitter mentioned using, NLPSolve does not return a result with HFloats. So not even double precision will be used under `eval`. Also, `eval` cannot be called from within `evalhf` in Maple 10.

In Maple 11, it still wouldn't be adequate, if greater working precision that the current Digits setting (eg, default of 10) were wanted:

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                               [0., [x = 3.]]

> evalf(subs(sol[2],1/x));

                                0.3333333333

> evalhf(eval(1/x,sol[2]));

                            0.333333333299999979

So, for the stated issue, in earlier Maple it's probably better to either raise Digits first. Or invoke evalhf with its indexed form, like evalf[n](...).

In Maple 15, double precision may get used, due to the presence of HFloats in the NLPSolve result.

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                    [0., [x = HFloat(3.0)]]

> evalf(subs(sol[2],1/x));

                   HFloat(0.3333333333333333)

> evalhf(eval(1/x,sol[2]));
                      0.333333333333333315

> lprint(sol);
[0., [x = HFloat(3.)]]

> kernelopts(version);

   Maple 15.00, X86 64 WINDOWS, Mar 21 2011, Build ID 582305

But if even greater working precision is required, then raising Digits above 15, or calling evalf[n](..) with n>15, may be required.

@Christopher2222 You do not have to have an existing (empty) worksheet saved, in order to open a new, empty one.

All you have to do is INTERFACE_WORKSHEET(':-display',':-text'=...) where the ... is a certain kind of long string of XML.

I don't think this avenue helps with Doug's question.

I think that this below shows that it is the kernel (not the GUI) which is preventing quit/done/stop from succeeding. It might be possible to bamboozle the kernel into thinking it was running in the commandline interface, but obvious attempts at that (like forcing a lie in the remember table of `IsWorksheetInterface`) don't seem to work here.

FromInert(_Inert_PROC(_Inert_PARAMSEQ(), _Inert_LOCALSEQ(),
                      _Inert_OPTIONSEQ(), _Inert_EXPSEQ(),
                      _Inert_STATSEQ(_Inert_STOP()),
                      _Inert_DESCRIPTIONSEQ(), _Inert_GLOBALSEQ(),
                      _Inert_LEXICALSEQ(), _Inert_EOP(_Inert_EXPSEQ())))();

Warning, done/quit/stop disabled.  Please use File->Close

What else could be tried? Hmm. Has anyone ever experimented with external-calling to Java on the libs bundled with the product? (Is there a public/private issue, with that? I'm also asking because I'd like to find the routine the GUI uses to encode embedded images -- it's not just base64.) I tried an OpenMaple help example once, and it acted as it it were CLI. Just wild thinking.

@Christopher2222 You do not have to have an existing (empty) worksheet saved, in order to open a new, empty one.

All you have to do is INTERFACE_WORKSHEET(':-display',':-text'=...) where the ... is a certain kind of long string of XML.

I don't think this avenue helps with Doug's question.

I think that this below shows that it is the kernel (not the GUI) which is preventing quit/done/stop from succeeding. It might be possible to bamboozle the kernel into thinking it was running in the commandline interface, but obvious attempts at that (like forcing a lie in the remember table of `IsWorksheetInterface`) don't seem to work here.

FromInert(_Inert_PROC(_Inert_PARAMSEQ(), _Inert_LOCALSEQ(),
                      _Inert_OPTIONSEQ(), _Inert_EXPSEQ(),
                      _Inert_STATSEQ(_Inert_STOP()),
                      _Inert_DESCRIPTIONSEQ(), _Inert_GLOBALSEQ(),
                      _Inert_LEXICALSEQ(), _Inert_EOP(_Inert_EXPSEQ())))();

Warning, done/quit/stop disabled.  Please use File->Close

What else could be tried? Hmm. Has anyone ever experimented with external-calling to Java on the libs bundled with the product? (Is there a public/private issue, with that? I'm also asking because I'd like to find the routine the GUI uses to encode embedded images -- it's not just base64.) I tried an OpenMaple help example once, and it acted as it it were CLI. Just wild thinking.

Thanks for adding the favorites. I am guessing that reading each others' favorites might become a popular pasttime.

I'm having trouble adding favorites to posts/answers which I've previously up-voted. (When I try, it just tells me that I've already voted, but doesn't add it as a favorite. It looks like interference with the up-vote icon right above the add-favorites icon.)

On my favorites page, any author's names seems to just be linked to www.maplesoft.com instead of to that user's profile page.

All the "Added" dates are not correct for anything submitted before "Primes v.2", ie. April 27, 2010.

Could the favourites please be sorted by the date of original submission of the respective posts?

The new Moderator badge is listed twice.

I'm sure that any issues can be fixed, so please don't take the above as criticism.

acer

@PatrickT Sure, it will work. But just like for other dsolve-numeric results with parameters, you'll have to set a parameter value before numerically computing.

Note the use of `forget`, since all of xt, yt, ft, and fft have been given option remember. This should be done each time you change any parameter value. If you don't clear their remember tables (by hand, or using `forget`) then their results will get remembered to correspond to the previous parameter value.

restart:

SOL := dsolve({ diff(x(t),t) = y(t)+x(t),
                diff(y(t),t) = y(t)-x(t),
                x(0) = p, y(0) = 1},
                'type' = numeric, 'output' = listprocedure,
                'parameters' = [p]):

xt := subsop(3=remember,eval( x(t), SOL )):
yt := subsop(3=remember,eval( y(t), SOL )):

f := (x,y) -> x^2+y^2:

ft := unapply( subs( X=xt(t), Y=yt(t), f(X,Y) ),
              t, 'proc_options' = ['remember','operator','arrow']):

temp := subs( unassigneddummy = f(xt(t),yt(t)), t -> evalf(unassigneddummy) ):
fft := subsop(3=op({op(3,eval(temp)),'remember'}), eval(temp) ):

ft(0.7);
Error, (in xt) parameters must be initialized before solution can be computed

fft(0.7);
Error, (in xt) parameters must be initialized before solution can be computed

SOL(parameters=[p=2.7]): #seems to work, but printed output looks odd

ft(0.7);

                              33.6176069516257

fft(0.7);

                              33.6176069516257

SOL(parameters=[p=3.4]):

ft(0.7); # rememebered!

                              33.6176069516257
fft(0.7); # remembered!

                              33.6176069516257


forget(ft):
forget(fft):
forget(xt):
forget(yt):

ft(0.7);

                              50.9333113513925
fft(0.7);

                              50.9333113513925

I wonder why doesn't searching for the words minimal and curvature find this interesting comment by DJKeenan.

acer

The posted code makes it easy to adjust the slopes at the knots. Thanks, Alec.

Let's recall that cubic spline interpolation, in a commonly used sense, minimizes curvature while producing a curve which is twice continuously differentiable. I guess people like this, if they hope to model smooth phenomena. There are lots of variants which can use piecewise cubics to produce something which may only be once continuously differentiable.

Christopher asked for the extra constraint that the slope be zero at certain points. Alec cleverly interpreted this to mean zero only at knots which are locally extreme (w.r.t the discrete adjacent knot values). This is a nice refinement on an instance of a Cardinal spline with parameter c=1, and also more sophisticated that what's called a Catmull–Rom (c=0) variant of the cubic Hermite spline.

Alec, I believe that your code does in fact use the average of the secants on either side, ie. the average of the slopes of the two line segments from knot i-1 to i and from knot i to i+1. (Your code gets this because the x-differences are all equal, I think... I tried to alter this to be true also for irregularly spaced knots, by using a scaled version of the so-called "finite difference" variant, see attached.)

But there's still a lot of leeway as to what tweaks one can make to the slope values at the knots. I'm attaching a worksheet in which the finite difference variant is applied to all internal knots which are not locally extreme. I keep Alec's twist to satisfy Chris, that slopes be zero at knots whose immediate neighbours are either both greater or lower in value. I scale each of the nonzero internal knots by parameter `c`. And then (unless I computed it wrongly) I find the `c` value which minimizes the total absolute curvature.

Here's a worksheet. The two animations should each produce the curve Alec showed above when c=1 (if I did it all right). Except that I changed the value at x=9 to be 7, instead of 6, so as to accentuate the hiccup that Christopher asked about in his comment above.

alec_modified_hermit.mw

Lots more variations are possible. One could use a set of c[i] and optimize all of them at once w.r.t the total absolute curvatue (each c[i] affects a pair of segments, just as the slopes themselves do). That might reduce hiccups and inflexions even more. One can think of even more interesting variants -- more than just those which have fancy sounding names.

I forgot to scale the slopes at the two end-knots by `c`, but changing the code is easy.

I didn't put any thought into how to compute the curvature more quickly (given that I know it's a piecewise, and I know the formula, etc). Corrections and comments are always welcome. The Optimization took about 10sec, and the plot a bit longer. There must be a fast way, but I didn't really try.

Alec's code layout is simple and easy to modify. I'm not sure whether the easier path to efficiency is to implement it using linear algebra, or to edit the procs to be Compilable.

acer

The spammed posts which have been cluttering up the Recent page make clear some other unfortunate aspects of the "new" Mapleprimes: a signigicant number of old posts have their formatting ruined, and are misclassified as Posts instead of Questions (and the replies as Comments instead of Answers).

The mis-format typically consists of each entire post or reply, including any pre-tag code, being crammed into one illegible paragraph.

I think that the main value of Mapleprimes is that it is a repository of great answers. But the broken search mechanism (indexing the "summary pages" as one problem has been called. Yikes!), the lack of good math typesetting, the difficulty in inserting 2D Math input except as attached full worksheets, the mis-format of old posts, the misclassification of old posts, and so on, take away from this value.

Mapleprimes could be more than just a free, one-stop (throw-away) Maple Support site, which is the epitome of a short-term benefit.

Fixing this site could transform it into an ongoing and easily accessed repository of Maple knowledge and solutions, and it could once again foster high quality technical discussion. Those are things which can bring additional long-term benefits such as increased web traffic, technical reputation, and good will.

Trying to think ahead: solving the problems don't (and really, really should not) need yet another re-implementation. Nobody wants another version 0.1 completely new site, twitter/facebook/cloud access, etc. Like for so many other things, the best course of action is to just amend what's already here with careful fixes to specific problems. Solid reliability using simple browser based access is better than version 0.1 implementation of new tech or even of old tech merged with new.

ps. I find that setting up Google Reader, with RSS feeds for the Question and Posts and Maplesoft Blog forums, is a convenient way to see all new postings without the spam getting too much in the way.

acer

@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