They are not all not the same `i`. The important distinction is due to some of those `i`s occuring within operator definitions.

The operators in a[1], a[2], and a[3] are procedures with no variable `i` declared in them. Also, in those procedures the variable `i` is not assigned to. See the section "Implicit Local Variables" in the help-page ?procedure for an explanation of how the value of such names is inferred when the procedure gets run.

As a consequence of this, when you run it ( eg. a[1](t) ) you get t*4 as it find a value of 4 for i in the surrounding level. That's because, after your loop finishes, the value of the `i` at your current level is 4. If you had assigned 17 to i, afterwards, then subsequent invocation of your operators would produce 17*i.

One way to work with this behaviour of Maple is to create an operator with a dummy name in it, and then to substitute the loop variable i's current value for that dummy each time through the loop.

> for i from 1 to 3 do
>     a[i] := subs(ii=i,t -> ii*t);
> end do;
                                a[1] := t -> t
 
                               a[2] := t -> 2 t
 
                               a[3] := t -> 3 t

For your example, you could also define that template operator just once, outside of the loop. The call to eval() below is there just because operators have so-called last-name evaluation rules, see ?last_name_eval for details of this side issue.

> restart:
> myop:=t -> ii*t:
> for i from 1 to 3 do
> a[i] := subs(ii=i,eval(myop)); # or eval(myop,1)
> end do;
                                a[1] := t -> t
 
                               a[2] := t -> 2 t
 
                               a[3] := t -> 3 t

You could also do this sort of thing more directly, using an Array instead of an array.

> restart:
> myop:=t -> ii*t:
> A := Array(1..3, (i)->subs(ii=i,eval(myop,1)) );

                       A := [t -> t, t -> 2 t, t -> 3 t]

Lastly, some easier ways to do it.

> restart:
> A := Array(1..3, (i) -> t -> i*t );

                      A := [t -> 1 t, t -> 2 t, t -> 3 t]

> vector(3, (i) -> t -> i*t );

                        [t -> 1 t, t -> 2 t, t -> 3 t]
> lprint(%);
array(1 .. 3,[(1)=(t -> 1*t),(2)=(t -> 2*t),(3)=(t -> 3*t)])

What happens during the second example just above is essentially this.

> f := i -> t -> i*t:
> a := array(1..3):
> for i from 1 to 3 do
>     a[i] := f(i);
> end do;
                               a[1] := t -> 1 t
 
                               a[2] := t -> 2 t
 
                               a[3] := t -> 3 t

acer

Not only can you do this with a simple Maple commands, you can also make a procedure that can be re-used to do it for any function you want.

Here's a re-usable routine to do the job.

> T := proc(F,a,b,inc)
> local i;
> for i from a to b by inc do
> printf("%a %a \n",i,F(i));
> end do;
> end proc:

And here's an example of using it.

> f := x -> x^2:

> T( f, 1.0, 3.0, 0.5 );

1.0 1.00
1.5 2.25
2.0 4.00
2.5 6.25
3.0 9.00

Of course, if you just need to do it once you can get similar effects with calls like this.

seq( printf("%a %a \n",0.5+0.5*i,f(0.5+0.5*i)), i=1..5 );

acer

The task you are trying to do is not uncommon. Your attempt was on the right track. What should work here is the `if` operator.

Matrix(3, (i,j) -> `if`(P[i,j]=1, x[i,j], P[i,j]) );

acer

There are two ways to do that sort of thing. The first is to qualify the command by appending `assuming` on the end of it. The second is to place an assumption on the name such as `n` by calling assume(). That second method cuts down on the typing, but makes it harder to use the name `n` elsewhere without any assumptions on it.

> sin(n*Pi) assuming n::integer;
                                  0
> restart:
> assume(n::integer):
> cos(n*Pi);

                                        n~
                                    (-1)

acer

I suspect that what is happening is that the internal iteration limit of the least-squares solver is being hit before the optimality tolerance (goal) is met, or before some numerical effect takes over the computation.

I've extracted what Statistics:-NonlinearFit passes to Optimization:-LSSolve , thanks to Maple debugger.

Once the call to LSSolve has been produced, it can be modified to alter the iterationlimit and also the choice of numerical method. See ?LSSolve for more.

First I run with method=sqp and iterationlimit=10000. I also set infolevel[Optimization] to 10, to get more information displayed about what it's doing.

> fnew := proc (w, z, needfi, ruser) local i; for i to 4 do z[i] := evalf(proc ( p, r, i) p[1]*(1-p[2]*cos(Pi/(r[i,1]+1)))^(1/2)-r[i,2] end proc(w,ruser,i)) end do end proc:
>
> fjac := proc (w, z, ruser, iuser) local i; for i to 4 do z[i,1] := evalf(proc (p, r, i) (1-p[2]*cos(Pi/(r[i,1]+1)))^(1/2) end proc(w,ruser,i)); z[i,2] := eval f(proc (p, r, i) -1/2*p[1]/(1-p[2]*cos(Pi/(r[i,1]+1)))^(1/2)*cos(Pi/(r[i,1]+1)) end proc(w,ruser,i)) end do end proc:
>
> infolevel[Optimization]:=10:
>
> sol := Optimization:-LSSolve(
  [2, 4],
  fnew,
  datapoints = (Matrix(4, 2, [[8.,-1.73120199199999991],[12.,-1.9665823319999999 9],[16.,-2.12740288800000021],[20.,-2.18699629199999990]], datatype = float[8])) ,
  output = tables,
  initialpoint = (Vector(2, [-9.,.699999999999999956], datatype = float[8])),
  objectivejacobian = fjac,
  method = sqp,
  iterationlimit = 10000 );

LSSolve:   calling nonlinear LS solver
SolveGeneral:   using method=sqp
SolveGeneral:   number of problem variables   2
SolveGeneral:   number of residuals   4
SolveGeneral:   number of nonlinear inequality constraints   0
SolveGeneral:   number of nonlinear equality constraints   0
SolveGeneral:   number of general linear constraints   0
PrintSettings:   feasibility tolerance set to   .1053671213e-7
PrintSettings:   optimality tolerance set to   .3256082241e-11
PrintSettings:   iteration limit set to   10000
PrintSettings:   infinite bound set to   .10e21
SolveGeneral:   trying evalhf mode
Warning, limiting number of major iterations has been reached
E04USA:   number of major iterations taken   10000
                       sol := [settingstab, resultstab]
 

> eval(sol[2]);

table(["visible" = ["objectivevalue", "solutionpoint", "residuals",
    "evaluations", "iterations"], "iterations" = 10000,
    "objectivevalue" = 0.0492868830043316011,
                      [-0.0488744033423635441]
    "solutionpoint" = [                      ],
                      [ -1732.28244861786970 ]
    "residuals" = [-0.241298347220388010, -0.0384274262517967369,
    0.110016320471836959, 0.163611172536819893]
    ])

> eval(a*sqrt(1-b*cos(Pi/(v+1))),[a=sol[2]["solutionpoint"][1],
>                                 b=sol[2]["solutionpoint"][2]]);
                               /                             Pi   \1/2
        -0.0488744033423635441 |1 + 1732.28244861786970 cos(-----)|
                               \                            v + 1 /

I notice that the residuals and the objective function above are decreasing quite slowly as the iterationlimit is increased. You may wish to play around with this a bit more. You could set the method to modifiednewton, or adjust the toleration of iterationlimit some more. If you decide that roundoff error is preventing it from attaining the best solution then you could increase Digits. What is interesting is that, depending on how I alter those options, very different solution points can be produced which each give quite close to the same residuals.

However, with method=modifiednewton and iterationlimit=100000 I got result without seeing the warning.

With those settings, I saw userinfo messages like this,

E04FCF:   number of major iterations taken   1474
E04FCF:   number of times the objective was evaluated   15018

and a final result that included this,

    "objectivevalue" = 0.0492807204246035,


                      [-0.0192100345224023965]
    "solutionpoint" = [                      ],
                      [ -11218.7180037680664 ]

> eval(a*sqrt(1-b*cos(Pi/(v+1))),[a=sol[2]["solutionpoint"][1],
>                                 b=sol[2]["solutionpoint"][2]]);
                               /                             Pi   \1/2
        -0.0192100345224023965 |1 + 11218.7180037680664 cos(-----)|
                               \                            v + 1 /

One possible conclusion from this (not discounting Georgios' comments in a post below about the merit of the nonlinear form that you have chosen to fit to) is that Statistics:-NonlinearFit could benefit from taking a new option like iterationlimit=<value>

acer

What dharr is trying to tell you, I believe, is that the inert form uppercase-e Eval() will display like a usual notation for "y evaluated at x=0". That uppercase Eval is the so-called inert form of eval. You should be able to assign calls like Eval(y,x=0) to variables, or simply use them in your Document, and so on.

On the other hand, the lowercase-e eval() can actually do the evaluation and is not inert. The two often look good in combination, on either side of an equation for example. Try executing this in Maple

y := cos(x);
Eval(y,x=0) = eval(y,x=0);

acer

There are probably better ways. It might also be possible to get the relational symbols from one of the palettes, and subs into that beast (lprinted to see it, say).

> f := (t,s) -> convert(StringTools:-SubstituteAll(
        cat("-T <= ",s," <= T"), "T", convert(t,string)),
                    name):

> f(2,"x");

                                -2 <= x <= 2

acer

Would you accept using T::ppositive , where ppositive were a new type?

Something like this, perhaps.

TypeTools:-AddType( ppositive,
  x->evalb(type(x,positive) or is(x,positive)=true) );

acer

If one could copy & paste your ode into Maple, as displayed in your post here on mapleprimes, then it would be easier to help.

But the ode appears as an image, in my firefox, from some maplenet server.  Does anyone here know how to get that into a new Maple session, without having to enter all that by hand?

acer

Here's an example.

sleep := proc(delay::nonnegative,expr::uneval,$)
local x, endtime;
   print(expr);
   if kernelopts(':-platform') = "unix" then
      ssystem(sprintf("sleep %d",delay));
   else
      endtime := time() + delay;
      while time() < endtime do end do;
   end if;
   seq(eval(x), x in args[2..-1]);
end proc:
                                                                                
sleep(5,int(x,x));

It could be even nicer if the seq Modifier described on the help-page ?parameter_modifiers allowed a sequence of uneval parameters. (Ambiguity could be resolved by having each argument taken separately and not as a sequence.)

Sidebar: if system/ssystem calls are disabled only in the Tools->Options section Standard GUI of Maple 11 then attempting to use them (like above, for Unix) results in a Warning but not necessarily an Error. But attempting to use system calls after -z is passed an an option to the Maple launcher (TTY or Standard) results in an Error.

acer

How about creating new procedures.

Sin := x->sin(Pi*x/180);
plot(Sin,0..360);

Or you could create and load a module which defines its own routines, so that the function names could appear the same.

my_trig := module()
option package;
export sin, cos;
  sin := x -> :-sin(Pi*x/180);
  cos := x -> :-cos(Pi*x/180);
end module:

with(my_trig);

plot(evalf@sin,0..360);
plot(sin(t),t=-180..180);
plot(sin,0..360); # bug, due to what??

You could add more routines to that module, even though above I only put in sin and cos. Simply define and export whichever you want.

If you feel really brave, you could also try this next trick below. I don't particularly recommend it, because apart from the danger of unprotecting sin, it doesn't get much more. It even seems to suffer from same problem as the module, for the last example.

restart:

Sin:=copy(sin):
unprotect(sin):
sin:=x->Sin(Pi*x/180):

plot(evalf@sin,0..360);
plot(sin(t),t=-180..180);
plot(sin,0..360); # again, why?

acer

The z in op(2,a) is an escaped local variable name.

You should be able to use convert(op(2,a),`global`) instead.

a:=FunctionAdvisor( sum_form, tan);
convert(tan(z),Sum) assuming convert(op(2,a),`global`);

acer

Your instincts are right. It can be done using seq(), and a loop which looks like,

S := S, ...

is to be avoided.

But I ought to ask, do you really need to create all 10 operators, which act by evaluating a piecewise?

Or do you just need that f[i](j) returns the value in M[i,j] , or possibly that f[i](y) returns M[i,trunc[y]] ?

There are two parts of the above question. The first is this.  Do you really need 10 operators, f[1]..f[10] (where f was presumably just a Maple table) or would it suffice to have one procedure f which was able be called in that indexed manner?

The second is this. Since the operators (10 individual ones, or 1 indexed proc, no matter) are apparently only going to be used to evaluate a piecewise (in x) at x=y then do you really need to utilize piecewise constructs for that? Maybe these operators don't need to use piecewise constructs at all (since those are, relatively speaking, expensive to evaluate at a specific point).

Some ideas follow, which you might find useful.

M := LinearAlgebra:-RandomMatrix(10,18):

f := proc(y)
#global M; # or by lexical scoping
local i;
  if type(procname,indexed) and type(op(1,procname),posint) then
    i:=op(1,procname);
  end if;
  M[i,trunc(y)];
end proc:

f[2](1);
f[2](1.3);
M[2,1];

Perhaps you really did need what you described, for some other reasons other than just the returned value M[i,j]. You noticed that this invocation below doesn't work, since seq() sees all the arguments flattened out (and hence invalid).

seq( j < x, MM[j,j]*`if`(j=1,0,M[i,j]), j=1..3 );

But maybe you could do something like this,

B := j < x, MM[j,j]*`if`(j=1,0,MM[i,j]);
seq( B, j=1..3 );

acer

Very nice, Axel.

But it seems to me that you have helped Maple out more than you should have had to, by "giving" it the FTOC result.

How about this: instead of,

D(F)(x): int(%,x=0..t) = F(t) - F(0);

you could have,

D(F)(x): Int(%,x=0..t) = int(%,x=0..t,continuous);

so that Maple does more of the work.

 acer

You forgot to define beta[1] as  value*deg .

By the way, I really prefer Units[Standard] to Units[Natural]. The Units[Natural] environment robs the global namespace of too many common symbols which one might want otherwise to use as variables.

> restart:

> with(Units[Standard]):

> alpha[1]:=4.25 * Unit(deg):
> theta[1]:=11.31 * Unit(deg):
> theta[2]:=78.69 * Unit(deg):
> alpha[2]:=-9.25 * Unit(deg):
> a:=1400.00:
> b:=509.90:

> beta[1]:=33.33 * Unit(deg): # I made that up.

> -b*sin(alpha[2]+alpha[1]-beta[1]-theta[2])+b*sin(beta[1]-theta[2]);

                              91.4313519

acer