In the earlier question which Marikyan has cited, Joe Riel showed a nice solution to this. But I'm curious, once converted to a Sum involving binomials, can it be more easily evaluated at particular posint n, without simplifying back to the original expression, I mean? And can we gracefully force Maple to curtail the summation at the point at which terms become zero?

restart:

z:=(x+1)^n;

                                   n
                            (x + 1) 

r:=convert(z,FormalPowerSeries,x);

               infinity                          
                -----                            
                 \         k                    k
                  )    (-1)  pochhammer(-n, k) x 
                 /     --------------------------
                -----         factorial(k)       
                k = 0                            

q := convert(r,binomial);

         infinity                                     
          -----                                       
           \                                          
            )        k                               k
           /     (-1)  binomial(-n + k - 1, -1 - n) x 
          -----                                       
          k = 0                                       

value(q);

                                   n
                            (x + 1) 

value(subs(n=4,r)); # ?

               infinity                          
                -----                            
                 \         k                    k
                  )    (-1)  pochhammer(-4, k) x 
                 /     --------------------------
                -----         factorial(k)       
                k = 0                            

value(subs(n=4,q)); # ?

             infinity                             
              -----                               
               \                                  
                )        k                       k
               /     (-1)  binomial(-5 + k, -5) x 
              -----                               
              k = 0                               

# How to get Maple to use this next fact, to curtail the summation?
simplify( binomial(-n+k-1, -1-n) ) assuming n::posint, k::nonnegint, k>n;

                               0

t:=subs(infinity=n,q); # since terms are zero for k>n

            n                                        
          -----                                      
           \                                         
            )       k                               k
           /    (-1)  binomial(-n + k - 1, -1 - n) x 
          -----                                      
          k = 0                                      

expand(subs(n=4,z));

                    4      3      2          
                   x  + 4 x  + 6 x  + 4 x + 1

value(subs(n=4,t));

                    4      3      2          
                   x  + 4 x  + 6 x  + 4 x + 1

One relatively easy way to save floating-point data inside the worksheet is to create a DataTable embedded component and to use the rtable (Matrix or Vector, say) which contains your data as the name of the location of the data for that component. This functionality was introduced in Maple 15.

Once the name of the rtable, X say, is associated with that component then you should be able to remove or comment out the lines which import that data and assign it to X. Upon saving and reopending the worksheet, the named rtable should still available for programatically accessing the stored entries.

You can even right-click and alter the properties of the component, to make it invisible (you can use DocumentTools:-SetProperty if you ever want to make it visible again). The data will also be available after restart.

datathing.mw

Do you mean that you want to return all indices k of a list L such that either the two entries L[k],L[k+1] are any pair of consecutive integers or they are 4 and 7?

restart: randomize():

Q:=proc(L)
  [seq(`if`(L[i+1]=L[i]+1 or (L[i]=4 and L[i+1]=7),i,NULL),i=1..nops(L)-1)];
end proc:

g:=rand(3..8):
LL:=[seq(g(),i=1..60)]:

ans:=Q(LL);

    [7, 12, 14, 16, 19, 21, 22, 24, 25, 30, 39, 40, 43, 55]

[seq([LL[x],LL[x+1]], x in ans)];

[[7, 8], [3, 4], [4, 5], [7, 8], [4, 7], [3, 4], [4, 5], [3, 4], 

  [4, 7], [4, 7], [4, 7], [7, 8], [3, 4], [7, 8]]

LL:=[4,7,3,13,10,14,4,7,3,13,10,14,4,7,3,13,10,14,4,7,3]:

ans:=Q(LL);

                         [1, 7, 13, 19]

[seq([LL[x],LL[x+1]], x in ans)];

                [[4, 7], [4, 7], [4, 7], [4, 7]]

You might need to adjust both or either of the working precision (15, used here) and tolerance (epsilon) of the call to evalf(Int()) inside procedure `f`.

with(DETools):

f:=proc(q) if not type(q,numeric) then return 'procname'(args); end if;
           evalf[15](Int(s->exp(cos(s)), 0..q, epsilon=1e-8)); end proc:

DEplot({diff(x(t),t)=f(t)},[x(t)],t=0..1,[[x(0)=1]],arrows=none);

You have declared `d` and `a[1]` as locals in your creation of the procedure `Sorozatok_1`. And the return value of calling that procedure contains those local names.

Such are known as "escaped locals", and are not the same as the names `d` and `a[1]` at the higher scope. In particular, when you call your procedure from the top-level the escaped names are not the same as the global names `d` and `a[1]`. Hence the sets of equations Solseq1 and Solseq2 are not identical. Their addresses differ, since the addresses of their indeterminate dependent names differ. And so the evalb test returns false.

Why did you declare `d` and `a[1]` as local? Perhaps you have not realized that it may not be necessary (or desirable). Your procedure does not assign to either of those names, incidentally.

You could remove `d` and `a[1]` from the local declaration statement inside Sorozatok_q1, and the evalb test would succeed.

Did you declare them as locals because you wanted to make sure that the names used in eq1 and eq2, as used inside the procedure body, were not already assigned? In that case, you could alter Sorozatok_q1 to accept two arguments, and then pass a pair of unassigned names whenever you call it.

restart:

Sorozatok_q1:= proc(A, B)
    local solved, n, l, k, s_n, eq1, eq2, MO, ret_length:
    solved := false:
    #randomize(): # commented out, b/c I hard code my Solseq2 below
    while ( (not solved) or hasfun([MO], RootOf) ) do
      n   := RandomTools:-Generate(integer(range=5..20)):
      l   := RandomTools:-Generate(integer(range=2..n)):
      k   := RandomTools:-Generate(integer(range= 1..(l-1) )):
      s_n := RandomTools:-Generate(integer(range=5..30)):
      eq1 := abs( A +(k-1)*B ) = abs( A+(l-1)*B );
      eq2 := n * ( A + (A + (n-1)*B) )/2 = s_n;
      MO  := solve( [eq1,eq2], [A,B], Explicit=true ):
      if not ( MO = [] ) then
         solved := true:
       end if:
    end do:
    ret_length := 4 + nops([MO]):
    return [ret_length, n, k, l, s_n, MO];
  end proc:

Out := Sorozatok_q1(`a[1]`,d);

 [                 [[       -130      20]  [       30       ]]]
 [5, 17, 7, 8, 30, [[a[1] = ----, d = --], [a[1] = --, d = 0]]]
 [                 [[        17       17]  [       17       ]]]

Solseq1 := op(map(convert,Out[-1],set));

          /       -130      20\    /       30       \ 
         { a[1] = ----, d = -- }, { a[1] = --, d = 0 }
          \        17       17/    \       17       / 

Solseq2 := {`a[1]` = -130/17, d = 20/17}, {`a[1]` = 30/17, d = 0};

\          /       -130      20\    /       30       \ 
         { a[1] = ----, d = -- }, { a[1] = --, d = 0 }
          \        17       17/    \       17       / 

evalb( Solseq1[1] = Solseq2[1] );
                              true

So if you'd already assigned to `d` and `a[1]` then you could call the procedure and pass it two other (unassigned) names of your choice.

I thought that I had a way for Maple 11, by temporarily converting to name, that would get a list of the sets in the same order as Maple 14. But it was flawed, so I deleted it.

Maybe an easier way is to build them as you want them, rather than to use `powerset`.

Do you mean this?

The first thing to realize is that `subs` may not be what you want for this kind of thing. Consider this example, which is only slightly more complicated:

restart:

subs( a*z=t, a*l*z );

                                    a l z

There just happens to be no "a*z" in this instance of the internal representation of a*l*z.

Now compare with the behaviour of the `algsubs` command.

restart:

algsubs( a*z=t, a*l*z );

                                     l t

Another point is that the reciprocal is stored differently. One commonly used trick is thus to substitute for both. The `algsubs` command doesn't allow multiple substitutions, but you can nest calls instead.

restart:

Q := 1/(a*l)*exp(a*l);

                                     exp(a l)
                                Q := --------
                                       a l   

algsubs( a*l=t, algsubs( (a*l)^(-1)=t^(-1), Q ) );

                                   exp(t)
                                   ------
                                     t   

Of course, I'm sure you realized that you could have simply substituted a=t/l in your examples, with the `subs` command. It happens to work for them, but might do too much for other examples. You might not want stand-alone instances of `a` (meaning: not products a*l) in your expression to all get replaced with t/l.

What do you plan on doing with such an explicit set of 13 formulas for the roots? What if it were 100s of pages long?

Is it the case that, eventually, you will try and evaluate any of those formulas (or some mathematical function of them) at numeric values of theta?

Let's suppose that, eventually, you would try and use numeric root (as a function of real-values theta). Instead of explicitly inverting that degree 13 polynomial in q[b] you could construct a procedure that takes a numeric value of theta and returns the roots in q[b].

You could then pick off the smallest root, or do something else with them, for the given value of theta.

eg.

expr := 472392*q[b]^13+(-10392624*theta+12597120*theta^2)*q[b]^12+(90069408*theta^2+110697192*theta^4-207222624*theta^3)*q[b]^11+(318718800*theta^6-394931376*theta^3-1213796664*theta^5
+1269661392*theta^4)*q[b]^10+(-3766288752*theta^5+4657375260*theta^6
+957739788*theta^4-1782380376*theta^7+14700393*theta^8)*q[b]^9+(207696*theta^10
+6288766560*theta^6-36914598*theta^9-1374185736*theta^5
+4234180392*theta^8-9185729832*theta^7)*q[b]^8+(19428192*theta^10-5820095184*theta^9
+10977610320*theta^8-6491850624*theta^7-116480*theta^11+1225901952*theta^6)*q[b]^7
+(-8848975200*theta^9+5424037888*theta^10+4454542848*theta^8-709841664*theta^7
+12076160*theta^11+25600*theta^12)*q[b]^6+(284107776*theta^8-2169162752*theta^9
+5185955328*theta^10-14362624*theta^12-3764956416*theta^11)*q[b]^5+(819093504*theta^10+2011156480*theta^12-95922176*theta^9-2311908352*theta^11
+5013504*theta^13)*q[b]^4+(-818036736*theta^13-276840448*theta^11
+810954752*theta^12-655360*theta^14+36864000*theta^10)*q[b]^3
+(-221405184*theta^13+240451584*theta^14+78905344*theta^12
-10616832*theta^11)*q[b]^2+(-45088768*theta^15+39190528*theta^14
-9437184*theta^13)*q[b]-3670016*theta^15+4194304*theta^16:

getqb := proc(th)
local ans;
    ans := fsolve(eval(expr, theta = th), q[b], args[2 .. -1]);
    if type([ans], list(complex(numeric))) then ans else NULL end if
end proc:

expr:=convert(expr,'horner',q[b]):

getqb(7.2);

getqb(2.3);

getqb(2.3, q[b]=0..1);

getqb(2.3, complex);

plot(t->[getqb(t)][1],1..5); # 30-40 secs or so

Look at march and/or LibraryTools.

The modern format is a single file with .mla as filename extension. The much older format is a pair of files that go together, with .lib and .ind filename extensions.

If you are content to act on all plots at once, after you have saved the worksheet/document to a file, then you might find Preben's script helpful.

G:=proc(n::posint)
  local t;
  if n>2 then
    t:=[h1,h$(n-3),h2]:
    LinearAlgebra:-BandMatrix([t,
                              [e$(n-1),e-b],
                               t]):
  end if;
end proc:

G(6);

                          [e   h1  0  0  0     0  ]
                          [                       ]
                          [h1  e   h  0  0     0  ]
                          [                       ]
                          [0   h   e  h  0     0  ]
                          [                       ]
                          [0   0   h  e  h     0  ]
                          [                       ]
                          [0   0   0  h  e    h2  ]
                          [                       ]
                          [0   0   0  0  h2  e - b]

eval(%, [b=0.5,e=2,h=1,h1=1.5,h2=1.5]);

                         [  2  1.5  0  0    0    0]
                         [                        ]
                         [1.5    2  1  0    0    0]
                         [                        ]
                         [  0    1  2  1    0    0]
                         [                        ]
                         [  0    0  1  2    1    0]
                         [                        ]
                         [  0    0  0  1    2  1.5]
                         [                        ]
                         [  0    0  0  0  1.5  1.5]


G:=proc(n::posint,b,e,h,h1,h2)
  local t;
  if n>2 then
    t:=[h1,h$(n-3),h2]:
    LinearAlgebra:-BandMatrix([t,
                              [e$(n-1),e-b],
                               t]):
  end if;
end proc:

G(6,0.5,2,1,1.5,1.5);

                         [  2  1.5  0  0    0    0]
                         [                        ]
                         [1.5    2  1  0    0    0]
                         [                        ]
                         [  0    1  2  1    0    0]
                         [                        ]
                         [  0    0  1  2    1    0]
                         [                        ]
                         [  0    0  0  1    2  1.5]
                         [                        ]
                         [  0    0  0  0  1.5  1.5]

You could pass only the subset of data which you want fitted. If unsorted then you could sort the data first, and extract the subset for which the xdata entries were within range.

It's not optimally efficient to convert to listlist, just to sort the data. But it's convenient. You could optimize that aspect, which I have not done. If your xdata and ydata are already sorted then you just have to extract the blocks of subdata.

An alternative is to weight the data, using weight values of 1 or 0 to mask out the influence of xdata entries which are out of the desired subrange.

restart:
randomize():
N:=100:

origp,origq:=0.0,3.0: # original range

xdata := LinearAlgebra:-RandomVector[row](N,generator=origp..origq):
ydata := Vector(N,i->evalf(exp(xdata[i]))):

c:=Statistics:-Fit( a*x^2+b, xdata, ydata, x );

p,q:=1,2: # desired subrange

L:=[seq([xdata[i],ydata[i]],i=1..N)]:
newL:=select(a->a[1]>=p and q>=a[1], sort(L,(a,b)->b[1]>a[1])):

newc:=Statistics:-Fit( a*x^2+b, newL, x );

newerc:=Statistics:-Fit( a*x^2+b, xdata, ydata, x,
                         weights=Vector(100,i->`if`(xdata[i]>=p
                                                      and q>=xdata[i],
                                                    1,0)) );


plots:-display(
  plot(L, style=point),
  plot(c,x=origp..origq,color=green),
  plot(newc,x=p..q,color=blue,thickness=3),
  plot(newerc,x=p..q,color=gold,thickness=3)
               );

linsplit:=E->`if`(type(E,'specfunc(anything,{Int,int})')
                    and type(op(1,E),`+`),
                  map(op(0,E),op(1,E),op([2..-1],E)),
                  E):

e := Int(x*y+y, [x,y]);
                                                      
                                 /   /
                                |   |
                          e :=  |   |  x y + y dx dy
                                |   |
                               /   /

> linsplit(e);
                                                                 
                        /   /               /   /
                       |   |               |   |
                       |   |  x y dx dy +  |   |  y dx dy
                       |   |               |   |
                      /   /               /   /

> value(%);

                                  2  2          2
                             1/4 x  y  + 1/2 x y

> value(e);

                                  2  2          2
                             1/4 x  y  + 1/2 x y

> f := Int(Int(x*y+y,x),y);                             

                                 /   /
                                |   |
                          f :=  |   |  x y + y dx dy
                                |   |
                               /   /

> # Could also use type anything in next line, as linsplit checks the type
> subsindets(f,'specfunc(anything,{Int,int})',linsplit);                  

                        /   /               /   /
                       |   |               |   |
                       |   |  x y dx dy +  |   |  y dx dy
                       |   |               |   |
                      /   /               /   /

> value(%);                                                               

                                  2  2          2
                             1/4 x  y  + 1/2 x y

Suppose that you precomputed N values of the time dependent "function". By interpolating, you should be able to get approximate values for t taken between those data points. That way, you could vary the number of frames, without having to construct `p` with the number of animation frames hard coded.

For example

restart:
N:=200:
L:=3.0:

Vt:=Vector(N,i->evalf((i-1)*2.5/(N-1))): # time values

T:=2*Pi:
Vv:=Vector(N,i->evalf(cos(Pi/2+T*(i-1)*2.5/(N-1)))): # time dep. values

p:=proc(t)
   if not type(t,numeric) then 'procname'(args);
   else CurveFitting:-ArrayInterpolation(Vt, Vv, t, method=spline);
   end if;
end proc:

W:=sin(x)/10:

plots:-animate(p(t)*W, x=0..L, t=0..2.5, frames=50, view=-1..1);

So, you can now change the frames=n value, and leave the rest as is. The number of frames should still evenly cover the time range of 0..2.5.

Note, however, that the time range is incorporated into Vv the Vector of time dependent values. (I just made up the function.)