Your second (1-line) way is a procedure call. The first way could also be made into a proc, and also be 1 line.

Your second way requires calling dsolve/numeric each time. That involves repeating all sorts of setup, initialization, deduction of problem class, algorithmic switching due to problem class, etc, etc. That is all unnecessary when one knows that the problem form will be the same, and that only some parameter will change. One purpose of the parameter option is to remove duplicated, unnecessary overhead.

The overhead costs both time and memory.

    |\^/|     Maple 13 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2009
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> ans:=(a,c)->dsolve({diff(y(x),x)+sin(x^2)*y(x)=a,y(0)=c},y(x),numeric):
> h:=rand(100):
> st,bu:=time(),kernelopts(bytesalloc):
> for ii from 1 to 1000 do
>   ans(h(),h())(5.7);
> end do: time()-st,kernelopts(bytesalloc)-bu;
                                5.424, 90161024
  
> quit
 
 
    |\^/|     Maple 13 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2009
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> ans:=dsolve({diff(y(x),x)+sin(x^2)*y(x\
> )=a,y(0)=c},y(x),numeric,parameters=[a,c]):
> h:=rand(100):
> st,bu:=time(),kernelopts(bytesalloc):
> for ii from 1 to 1000 do
>   ans(parameters=[h(),h()]);
>   ans(5.7);
> end do: time()-st,kernelopts(bytesalloc)-bu;
                                1.613, 38003920

   |\^/|     Maple 13 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2009
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> ans:=dsolve({diff(y(x),x)+sin(x^2)*
> y(x)=a,y(0)=c},y(x),numeric,parameters=[a,c]):
> h:=rand(100):
> nans:=proc(a,b,x) ans(parameters=[a,b]); ans(x); end proc:
> st,bu:=time(),kernelopts(bytesalloc):
> for ii from 1 to 1000 do
>   nans(h(),h(),5.7);
> end do: time()-st,kernelopts(bytesalloc)-bu;
                                1.632, 38003920

There may also be some technical reasons underneath. For example, using globals is not thread-safe.

acer

Hopefully, this is of some use. There may be easier ways.

> interface(warnlevel=0):
> with(DETools): with(plots): with(plottools):
> Diam := .4: m := 100: g := 9.8: c1 := 0.155e-3*Diam:
> c2 := .22*Diam^2: ynot := 50000*.3048: Wx := 60*.44704: Wy := 0:
> Vwx := diff(x(t), t)+Wx:
> Vwy := diff(y(t), t)+Wy:
> theta := arctan(Vwy/Vwx):
> Vwmag := sqrt(Vwx^2+Vwy^2):
> Fdragmag := -c1*Vwmag-c2*Vwmag^2:
> Fdragx := Fdragmag*cos(theta):
> Fdragy := Fdragmag*sin(theta):
> yeq := m*(diff(y(t), t, t)) = -m*g+Fdragy:
> xeq := m*(diff(x(t), t, t)) = Fdragx:
> ics := (D(x))(0) = ground, (D(y))(0) = 0, y(0) = 50000*.3048, x(0) = 0:
> ohthehumanity := dsolve({ics, xeq, yeq}, numeric, [y(t), x(t)],
>                         output = listprocedure):
> hite := rhs(ohthehumanity[2]):
> pos := rhs(ohthehumanity[4]):
>
> Z:=proc(GG,TT) global ground;
>   if not(type(GG,numeric)) then
>     return 'procname'(args);
>   else
>     ground:=GG;
>   end if;
>   pos(TT)^2 + hite(TT)^2;
> end proc:
>
> Zgrad := proc(V::Vector,W::Vector)
>   W[1] := fdiff( Z, 1, [V[1],V[2]] );
>   W[2] := fdiff( Z, 2, [V[1],V[2]] );
>   NULL;
> end proc:
>
> sol:=Optimization:-Minimize(Z, 'objectivegradient'=Zgrad,
>   'method'='nonlinearsimplex' ,initialpoint=[100,100]);
                                             [107.518346261289281]
           sol := [0.0000301543312111207006, [                   ]]
                                             [104.668066005965386]
 
>
> pos(sol[2][2]), hite(sol[2][2]);
                0.00454084653245279135, 0.00308788665268533435
 
> ground,sol[2][2];
                   107.518346261289281, 104.668066005965386
 
>
> Digits:=30:
> sol:=Optimization:-Minimize(Z, 'objectivegradient'=Zgrad,
>   'method'='nonlinearsimplex' ,initialpoint=[100,100]);
sol := [
 
                                       -13  [107.518073058989098792559364875]
    0.898257246759777661510402862893 10   , [                               ]]
                                            [104.668080856614961089921927934]
 
>                                                                                 
> pos(sol[2][2]), hite(sol[2][2]);
                                    -6                          -6
            -0.405056976848783279 10  , -0.339596677889630882 10
 

EDIT: the objective gradient routine `Zgrad` is not necessary when using the nonlinearsimplex method (which doesn't use derivatives, if I recall). I created it only because it helps when using some other 'method's of Minimize. I had to try various methods before I found one that got the result nicely. Minimize's ProcedureForm or MatrixForm calling sequences sometimes run into trouble unless an objective gradient (and/or constraintjacobian) routine is provided explicitly for "black-box" objective (and/or constraint) procedures which cannot be automatically handled by codegen[GRADIENT] (and/or codegen[JACOBIAN]. The same result returns from simply,

Optimization:-Minimize(Z,
  'method'='nonlinearsimplex' ,initialpoint=[100,100]);

The nature of routine Z is to find a matching pair 'ground' and time values for which both `hite` and `pos` (Y and X) are simultaneously zero.

acer

Hopefully, this is of some use. There may be easier ways.

> interface(warnlevel=0):
> with(DETools): with(plots): with(plottools):
> Diam := .4: m := 100: g := 9.8: c1 := 0.155e-3*Diam:
> c2 := .22*Diam^2: ynot := 50000*.3048: Wx := 60*.44704: Wy := 0:
> Vwx := diff(x(t), t)+Wx:
> Vwy := diff(y(t), t)+Wy:
> theta := arctan(Vwy/Vwx):
> Vwmag := sqrt(Vwx^2+Vwy^2):
> Fdragmag := -c1*Vwmag-c2*Vwmag^2:
> Fdragx := Fdragmag*cos(theta):
> Fdragy := Fdragmag*sin(theta):
> yeq := m*(diff(y(t), t, t)) = -m*g+Fdragy:
> xeq := m*(diff(x(t), t, t)) = Fdragx:
> ics := (D(x))(0) = ground, (D(y))(0) = 0, y(0) = 50000*.3048, x(0) = 0:
> ohthehumanity := dsolve({ics, xeq, yeq}, numeric, [y(t), x(t)],
>                         output = listprocedure):
> hite := rhs(ohthehumanity[2]):
> pos := rhs(ohthehumanity[4]):
>
> Z:=proc(GG,TT) global ground;
>   if not(type(GG,numeric)) then
>     return 'procname'(args);
>   else
>     ground:=GG;
>   end if;
>   pos(TT)^2 + hite(TT)^2;
> end proc:
>
> Zgrad := proc(V::Vector,W::Vector)
>   W[1] := fdiff( Z, 1, [V[1],V[2]] );
>   W[2] := fdiff( Z, 2, [V[1],V[2]] );
>   NULL;
> end proc:
>
> sol:=Optimization:-Minimize(Z, 'objectivegradient'=Zgrad,
>   'method'='nonlinearsimplex' ,initialpoint=[100,100]);
                                             [107.518346261289281]
           sol := [0.0000301543312111207006, [                   ]]
                                             [104.668066005965386]
 
>
> pos(sol[2][2]), hite(sol[2][2]);
                0.00454084653245279135, 0.00308788665268533435
 
> ground,sol[2][2];
                   107.518346261289281, 104.668066005965386
 
>
> Digits:=30:
> sol:=Optimization:-Minimize(Z, 'objectivegradient'=Zgrad,
>   'method'='nonlinearsimplex' ,initialpoint=[100,100]);
sol := [
 
                                       -13  [107.518073058989098792559364875]
    0.898257246759777661510402862893 10   , [                               ]]
                                            [104.668080856614961089921927934]
 
>                                                                                 
> pos(sol[2][2]), hite(sol[2][2]);
                                    -6                          -6
            -0.405056976848783279 10  , -0.339596677889630882 10
 

EDIT: the objective gradient routine `Zgrad` is not necessary when using the nonlinearsimplex method (which doesn't use derivatives, if I recall). I created it only because it helps when using some other 'method's of Minimize. I had to try various methods before I found one that got the result nicely. Minimize's ProcedureForm or MatrixForm calling sequences sometimes run into trouble unless an objective gradient (and/or constraintjacobian) routine is provided explicitly for "black-box" objective (and/or constraint) procedures which cannot be automatically handled by codegen[GRADIENT] (and/or codegen[JACOBIAN]. The same result returns from simply,

Optimization:-Minimize(Z,
  'method'='nonlinearsimplex' ,initialpoint=[100,100]);

The nature of routine Z is to find a matching pair 'ground' and time values for which both `hite` and `pos` (Y and X) are simultaneously zero.

acer

That's an interesting find.

That routine works by computing (approximations, while Digits=8, to) the zeros, extrema, and inflections of the expression. Then it reverts Digits and determines convexity/concavity by taking signum(eval(diff(expr,x,x),x=mid)) for 'mid' the midpoint between each pair of those points of interest.

I see using signum as mildly superior to using `is` directly, since signum can sometimes use evalr and a few other techniques as well as `is`.

I note that these Calculus1 routines (including `Roots`) demand a finite range. I'm not so sure about the wisdom of setting of Digits to 8 while finding critical points.

While determining subintervals on which the expression is (piecewise) convex is interesting, that routine doesn't quite bring all the same inferences. For one thing, it uses only a finite range. It also doesn't actually make an statement about adjoining subintervals, I suspect. Sure x^2 is convex from a..0 and also from 0..b, which it shows. But that's not as strong a statement as saying that, yes, "x^2 is convex from a..b" or yes, "x^2 is convex".

acer

In order to use the test of the 2nd derivative being nonnegative, the expression should be twice differentiable. That is one way of seeing why it falls down for 1/(x*x), which is not even once differentiable.

The code might be altered to check differentiability, running limit of the derivative(s) against values returned by discont or fdiscont. (There is an existing "property" named 'differentiable' known to is&assume&assuming, but it's for "functionals".)

If the method is edited to force use of the definition, avoiding the 2nd derivative test, then one gets behaviour like,

> isconvex(1/(x*x));
isconvex:   falling back to the definition
                                     false
 
>
> isconvex(-sqrt(x^2));
isconvex:   falling back to the definition
                                     false
 
> isconvex(1/(x*x),0,infinity);
isconvex:   falling back to the definition
                                     true
 
> isconvex(-sqrt(x^2),-infinity,0);
isconvex:   falling back to the definition
                                     true
 
> isconvex(abs(x));
isconvex:   falling back to the definition
                                     FAIL

10-15 years ago I used to hear experts say that it was only a misdemeanor if Maple's answer was wrong only on a set of measure zero. Nowadays more Maple routines return and accept piecewise results, and in turn people have come to expect more of the same. Given that not all of Maple is prepared to handle piecewise(), what havoc would be caused by having diff(1/x^2,x) return piecewise(x=0,undefined,-2/x^3) ?

acer

I value your comments, Alec.

acer

As is, it fails for lots of things. Obviously I was implying that when I wrote that there is lots of room for improvement.

Note also that there was a comment in the code about whether/how to normalize. Testing for differentiability (another tricky area) is also possible.

Testing for convexity is similar to lots of computer algebra tasks in that there are lots of theoretical and practical holes. Unless someone starts it off, even if only for discussion, there will always be nothing.

"Pragmatism. Is that all you have have to offer?" - Guildenstern

acer

I was just minimizing a squared expression, so that if NLPSolve found values where the objective was (approximately zero) then that could be a zero of the non-squared expression. (I used NLPSolve as an alternative approach, because fsolve requires as many expressions/functions as variables. You might have just one function but several variables.)

I have to rush off home. Since I wasn't sure of your exact query, I didn't make use of x(t) in my earlier code. I used just y(t). I'll look again, if nobody posts a methodology for your clarified task before then.

acer

I was just minimizing a squared expression, so that if NLPSolve found values where the objective was (approximately zero) then that could be a zero of the non-squared expression. (I used NLPSolve as an alternative approach, because fsolve requires as many expressions/functions as variables. You might have just one function but several variables.)

I have to rush off home. Since I wasn't sure of your exact query, I didn't make use of x(t) in my earlier code. I used just y(t). I'll look again, if nobody posts a methodology for your clarified task before then.

acer

There's a lot of scope for improving this code below. (A probabilistic fallback might also be worthwhile. I suspect that a fallback that using testeq & signum might not work; testeq is rather limited.)

Note that `is` can depend upon Digits.

> isconvex := proc(expr,A::realcons:=-infinity,B::realcons:=infinity)
> local a,b,t,x,use_expr,res;
>   x := indets(expr,name) minus {constants};
>   if nops(x) = 1 then
>     x := op(x);
>   else error;
>   end if;
>   # For more than one variable, one could test
>   # positive-definiteness of the Hessian.
>   res := is( diff(expr,x,x) >=0 )
>            assuming x>=A, x<=B;
>   if res<>FAIL then
>     return res;
>   else
>     userinfo(1,'isconvex',`falling back to the definition`);
>     # Is it better to expand, or simplify, or put
>     # Normalizer calls around the two-arg evals or
>     # their difference, oder...?
>     use_expr:=expand(expr);
>     is( eval(use_expr,x=(1-t)*a+t*b) <=
>         (1-t)*eval(use_expr,x=a)+t*eval(use_expr,x=b) )
>       assuming a<=b, t>=0, t<=1, b<=B, a>=A;
>   #else
>   ##Is global optimization another fallback?
>   end if;
> end proc:

> # note: this shows that none of the examples
> # below had to fall back to the defn after failing
> # a 2nd deriv test. But many such will exist.
> infolevel[isconvex]:=1:

> isconvex( y^2 );
                                     true

>
> isconvex( y^3 );
                                     false

> isconvex( y^3, 0 );
                                     true

>
> isconvex( (y-3)^2 );
                                     true

>
> isconvex( (y-3)^4 );
                                     true

>
> isconvex( exp(y) );
                                     true

> isconvex( -exp(y) );
                                     false

>
> isconvex( -sin(y), 0, Pi );
                                     true

acer

Perhaps someone will add it as the block comment syntax for Maple on wikipedia's language syntax comparison page (both here, and here).

acer

Perhaps someone will add it as the block comment syntax for Maple on wikipedia's language syntax comparison page (both here, and here).

acer

I'm not sure what you mean by "solve for ground". What are the variables?

Are you trying to find the value of 'ground' that satisfies a given y(t) at a given t?

> Z:=proc(TT,GG,YY) global ground; ground:=GG; hite(TT)-YY; end proc:

> fsolve( 'Z'(33.,gr,11853.) );
                                  106.7378019

Or are you trying to find both a 'ground' and a 't' that satisfy a given y(t)?

> Optimization:-Minimize( 'Z'(t,gr,11700)^2, gr=0..200, t=0..100 );
                        -16
[0.240100000000000009 10   ,

    [gr = 75.8416010127991882, t = 33.6050437589785673]]

> fsolve( 'Z'(33.605044,gr,11700.) );
                                  75.84162069

You mentioned `minimize`. What do you wish to minimize, in terms of what else?

acer

I'm not sure what you mean by "solve for ground". What are the variables?

Are you trying to find the value of 'ground' that satisfies a given y(t) at a given t?

> Z:=proc(TT,GG,YY) global ground; ground:=GG; hite(TT)-YY; end proc:

> fsolve( 'Z'(33.,gr,11853.) );
                                  106.7378019

Or are you trying to find both a 'ground' and a 't' that satisfy a given y(t)?

> Optimization:-Minimize( 'Z'(t,gr,11700)^2, gr=0..200, t=0..100 );
                        -16
[0.240100000000000009 10   ,

    [gr = 75.8416010127991882, t = 33.6050437589785673]]

> fsolve( 'Z'(33.605044,gr,11700.) );
                                  75.84162069

You mentioned `minimize`. What do you wish to minimize, in terms of what else?

acer

Despite the new error messages you mention, I believe that one of those approaches will still be needed.

The "procedure form" of Maximize is somewhat prone to doing that (Test-1) especially if a procedure to compute the objective gradient (and constraint jacobian, if constraint procs are supplied) is not explicitly supplied. A "known" bug.

As for your Test-2 and Test-3 attempts, it might be nice to have your .xls data, to figure out what went wrong. (Could it be that that `ln` call is producing a non-real?)

acer