It's an interesting question: how to project the 3D spacecurve onto the x-y plane along with the field arrows corresponding to 3D points (only) along that curve. I'm not sure how to do that easily, without just building and inserting a bunch of arrows by hand.

The 3D field arrows corresponding to 3D points along a given vertical line (x and y both constant) do not all point in a constant x-y direction. So there is no obvious correct choice to make for the field arrows elsewhere in the x-y plane (except maybe for the x-y points of the projection of the 3D spacecure).

Here's the 3D plot. Perhaps it will convey the message you hoped for from the 2D portrait.

restart:

DE:={D(x)(t)=y(t)-z(t),D(y)(t)=z(t)-x(t),D(z)(t)=x(t)-y(t)}:

dsn:=dsolve(DE union {x(0)=1,y(0)=0,z(0)=2}, numeric, output=listprocedure):

plots:-display(
  plots:-odeplot(dsn, [x(t),y(t),z(t)], t=-2..2),
  plots:-fieldplot3d(subs(x(t)=x,y(t)=y,z(t)=z,
                        eval([diff(x(t),t),diff(y(t),t),diff(z(t),t)],
                             convert(DE,diff))),
                   x=-0.5..2.5,y=-0.5..2.5,z=-1.0..3.0) ,axes=box);

acer

As given, the worksheet attempts to do symbolic solving of the DE. This results, after using quite a lot of time and memory, in a very long solution which could be of less practical use than you might want. (Trying to do symbolic solving of equations with floating-point coefficients is often a warning that something is amiss, IMO.)

You could try something like this, at the end,

sol3 := dsolve([eqB, eqC, eqD, eqE, eqF,
                speed(0) = 0, V(0) = 0.7047824398e-1],
               [speed(t), Tm(t), Tw(t), V(t), omega(t)],
               numeric, output = listprocedure):

spd := subs(sol3, speed(t)):

pos := proc(T) evalf(Int(spd,0..T,epsilon=1e-3)); end proc:

st:=time():
plot([spd,pos], 0..20);
time()-st;

That's computing both the speed and position numerically. Now, didn't Markiyan recently answer a related earlier question by you, citing an old post in which Joe Riel suggested augmenting the DE system so as to allow dsolve/numeric itself to do the integration (as opposed to using evalf/Int to do it with numeric quadrature)? That was probably good advice, because dsolve/numeric is supposed to be clever about generating a sequence of computed values, whereas the above scheme is going to repeat effort when integrating numerically from T=0 out to T=blah, then once more from T=0 to further point T=blech, etc, etc for all the points needed for the plot.

If you choose to stick with a numeric scheme for calculating position like I show here, then Patrick T's related advice (also to you) seems to save about 1/4 the computational effort while plotting (but this may involve luck, that the repeated numeric quadratures re-use points in time -- instead of using all ugly unre-usable gaussian points, say...). For example, to compare, this seems to take 3/4 the time of the above plot,

spd := subsop(3=remember, subs(sol3, speed(t))):
st:=time():
plot([spd,pos], 0..20);
time()-st;

Your attached worksheet seems to contain some mysterious hidden call which starts the ODE assistant (maplet), if run using the triple-exclam from the menubar. Weird. Here's an edited version using what I've mentioned. (I didn't hunt for the mysterious thing.)

Do you still need to resolve the data table interpolation question, afer trying this? On-the-fly interpolation shouldn't be faster than your degree 6 polynomial fit. The only practical benefit I can think of for still trying it is if you feel it would interpolate the data better. Less wiggle, or what have you.

Gear_Ratio_Study-ne.mw

acer

Yes.

You can use procedures for this. Put all the computational code inside a procedure whose parameters are the variables you wish to vary. Then, instead of assigning some values to those variables at the top-level (ie. instead of outside the proc), call your procedure with those values. Or call it again with different values. And so on.

But it's not a trick! It's the natural progression from the basic assign-values-then-loop-while-doing-stuff-then-plot-and-rely-on-menubar-triple-exclam-to-execute-all approach to using Maple.

Sure, you might have to cut some of the material and paste it into just one execution block to contain the procedure.

If you run into problems in accomplishing this, you could upload the worksheet in a new comment to this thread and tell us what you want to vary.

acer

Try command-completion, while typing. That also works for me in Maple 12 (which I believe you still use.)

In other words, do it like so:

• type in diff
• hit the command-completion key sequence (Ctl-space on Windows)
• Choose either the third or second entry in the list the pops up, ie. diff (inline), or diff (inline partial)
• Don't put the expression to be differentiated in the numerator of the "derivative symbol", put it afterward.

You want to end up with some like

and not something like

If you want to control how much of an expression gets the applied to it, then use round brackets to delimit the scope. Note the parsed meaning,

rather than being equal to,

You can also use the appropriate entry from the Expression palette.

acer

DEtools[DEplot] is really just out of date and inefficient in its technique.

A leaner and faster approach is to use plots:-fieldplot alongside a sequence of (with varying ICs) of plots:-odeplot calls.

I'm going to branch off this Answer, and show both plots:-animate and Embedded Component implementations of both the old DEplot and that new "combined" approach.

As for saving the worksheet and huge memory use by the GUI, well you can solve that by removing all the output from the worksheet before ever saving it. You might object and say, "but my animation takes 45 sec to create on the very fastest machine! I don't want to have to wait every time I open it!" To which I'd respond, just look how much faster the combined approach is: maybe 40-60 times faster to create the whole animation, for the ODE system in this thread.

And the Embedded Component implementation saves, with output included, with no problem. And the combined approach allows it to display smoothly and quickly, while allowing for nice manipulation of two parameters.

acer

> sum(1/(2*z)^2, z = 1 .. infinity);


                             1    2
                             -- Pi 
                             24    

It's pretty standard, to use a form like 2*k and 2*k-1 (or 2*k+1 for the odd posints if the summation begins from 0) for representing either even or odd integers. And such representations may be taken advantage of by `sum` or `product`. That's not really true of some other techniques such as multiplying the general term by a function call which evaluates to 1 or 0 according to oddness, since `sum` would mostly see such as a black box (or get it wrong).

For the example above, the factor 1/4 can be taken outside the summation, reducing it easily to the Pi^2/6 result for the sum over all posints. So this example doesn't really illustrate `sum` "understanding" the even-posint term per se. Another example might illustrate the advantage more.

> sum(x^(2*k-1)/factorial(2*k-1), k = 1 .. infinity);

                            sinh(x)

Compare with, say,

> sum('`mod`(k,2)'*x^k/factorial(k), k = 1 .. infinity); # wrong

                            exp(x) x

> iver:=proc(e)
>    if e::{constant,relation(constant)} then
>       if e then 1 else 0 end if;
>    else 'procname'(e);
>    end if;
> end proc: 

> sum(iver(k::odd)*x^k/factorial(k), k = 1 .. infinity); # no result

                    infinity               
                     -----                 
                      \                   k
                       )    iver(k::odd) x 
                      /     ---------------
                     -----   factorial(k)  
                     k = 1                 

So only one of the above three attempts gets the right answer. But if you replace `infinity` by literal numbers (in which case it may not do symbolic summation) then they concur.

> sum(x^(2*k-1)/factorial(2*k-1), k = 1 .. (7+1)/2);

                      1  3    1   5    1    7
                  x + - x  + --- x  + ---- x 
                      6      120      5040   

> sum('`mod`(k,2)'*x^k/factorial(k), k = 1 .. 7);

                      1  3    1   5    1    7
                  x + - x  + --- x  + ---- x 
                      6      120      5040   

> sum(iver(k::odd)*x^k/factorial(k), k = 1 .. 7);

                      1  3    1   5    1    7
                  x + - x  + --- x  + ---- x 
                      6      120      5040   

(I wrote this up, but didn't post it, after reading this posting.)

acer

You didn't say what dimension is your data set.

For the case of 3-d data (z-values for each 2-d x-y point in a grid) you might find this of some use.

If you have the simpler situation of just a 1-d set of input points (and thus 2-d data, with y-values for each 1-d point) then it should be pretty simple to modify that approach.

acer

You didn't say what dimension is your data set.

For the case of 3-d data (z-values for each 2-d x-y point in a grid) you might find this of some use.

If you have the simpler situation of just a 1-d set of input points (and thus 2-d data, with y-values for each 1-d point) then it should be pretty simple to modify that approach.

acer

As far as I can see, the points of significance claimed in this post are all factually wrong. Please explain, if I have misundertood the central points.

1.) The output doesn't look like that, in either Document or Worksheet, in the Standard GUI of Maple 15, 14, 13, or 12. In a Worksheet, there is output like the following (the closest I can see to what was claimed in this Post). Note that those below are echoed assignments (with :=), and not equations (with =) as in the posted claim.

> a:=[2,3,4]:
> printlevel:=2:
> for i in a do
>   i;
> od;

                                    i := 2

                                       2

                                    i := 3

                                       3

                                    i := 4

                                       4

2.) The output is not "doubled" -- at least not in a Worksheet. And it's not even duplicated. Successive lines only happen to look similar (or be exactly the same, in a Document) because the penultimate line is `i;`.

If instead that line were i+k, then it wouldn't look like a duplicate anywhere. It would look like,

> a:=[2,3,4]:
> printlevel:=2:
> for i in a do
>   i+k;
> end do;


                                   i := 2

                                    2 + k

                                   i := 3

                                    3 + k

                                   i := 4

                                    4 + k

3.) Where's the bug? All I see is that the incremented assignment of the loop index variable is being echoed. What's wrong with that? (It does actually get assigned, as following exit from the loop the name `i` has value. That is normal. (In the commandline interface, taking into account the syntax change of "od;" versus "end do;", the results are the same as here, in Maple V R4 as they are in Maple 14 and Maple 15.)

If you don't want assignment of the loop indexed echoed as output, then don't raise the printlevel from its default.

acer

Try not making that previous assignment to R. That may necessitate other changes to the whole code (unseen here).

We cannot guess all previous and subsequent situations and tasks of the OP, but 2-argument eval, etc, could help. Of course, this approach may require revision and altering of previous and subsequent subtasks (so cooking up such is a spurious counter-point, unless they have no similarly themed solution!).

Assignment is overrated.

> m1 := Matrix([[X/R, Y/R], [Z/R, t/R]]);

                                      [X  Y]
                                      [-  -]
                                      [R  R]
                                m1 := [    ]
                                      [Z  t]
                                      [-  -]
                                      [R  R]

> m2 := eval(m1, [X = 1, Y = 1, Z = 1, t = 1]);

                                      [1  1]
                                      [-  -]
                                      [R  R]
                                m2 := [    ]
                                      [1  1]
                                      [-  -]
                                      [R  R]

> eval(m2, R = sqrt(X^2+Y^2+Z^2));

                 [         1                    1         ]
                 [-------------------  -------------------]
                 [              (1/2)                (1/2)]
                 [/ 2    2    2\       / 2    2    2\     ]
                 [\X  + Y  + Z /       \X  + Y  + Z /     ]
                 [                                        ]
                 [         1                    1         ]
                 [-------------------  -------------------]
                 [              (1/2)                (1/2)]
                 [/ 2    2    2\       / 2    2    2\     ]
                 [\X  + Y  + Z /       \X  + Y  + Z /     ]

acer

"Doctor, it hurts when I do this."

You've accidentally used plain = in two places, when trying to update amax and tmax. But that just creates equations and doesn't assign anything to either amax or tmax.

You probably want := instead, so that those are assignments.

You may also want to edit the printf's format, as amax will not be an integer, as things stand.

acer

Fortran(subs(unknown=cg,[codegen[optimize](f)]));

acer

I am not a big fan of this much code abstraction. But I suspect that this eventually goes out to external, compiled code. (It is not only kernel builtins for which one cannot see the source.)

Sometimes a module is implemented to be appliable [sic], by putting code inside a procedure which is its local ModuleApply member. This seems true of Finance:-BrownianMotion. You can call it with arguments, and it is a module, and its ModuleApply procedure is where the work first starts getting done. (Sorry,if you knew that already.)

> restart:
> kernelopts(opaquemodules=false):

> print(Finance:-BrownianMotion:-ModuleApply);

proc(x0::{Vector, realcons})
  if type(x0, ':-realcons') then
    return BrownianMotion1D(args)
  else
    return BrownianMotion(args)
  end if;
end proc;

> print(Finance:-BrownianMotion:-BrownianMotion1D);

 module () local GetStochasticVariables, GetStochasticFactors, ModuleApply, 

    ProcessOptions, StochasticFactorMap; end module

> print(Finance:-BrownianMotion:-BrownianMotion1D:-ModuleApply);

proc()
local parms, rest;
  parms, rest := ProcessOptions(args);
  if 0 < nops(rest) then error "unexpected parameters: %0", op(rest) end if;
  try
    return Utilities:-MakeStochasticProcess(':-_X', parms[1], parms[2])
  catch:
    error 
  end try;
end proc;

> print(Finance:-Utilities:-MakeStochasticProcess);

     proc()
     local x, p, y, t;
       if type([args[1 .. 2]], [':-name', {
         ':-_FinanceObject'(':-MapleMarkovProcess'), 
         ':-_FinanceObject'(':-MapleMarkovProcess1D')}]) then
         y := `tools/genglobal`(args[1]);
         t := `tools/genglobal`(':-_StochasticProcess');
         assign(t, args[2]);
         setattribute(y, t);
         protect(y);
         return y;
       elif type([args[1 .. 3]], [':-name', ':-string', ':-table']) then
         return procname(args[1], MakeFinanceObject(args[2 .. -1]))
       end if;
     end proc;

> showstat(Finance:-Utilities:-MakeFinanceObject);

Finance:-Utilities:-MakeFinanceObject := proc(extfunc::string, parameters::table, printfunc::procedure := proc () module () end module end proc)
local m, i, p, q, f, data, external;
   1   data := Finance:-Utilities:-GetImplementation(parameters);
   2   p := [op(sort(map2(op,1,[indices(parameters)]))), ModulePrint];
   3   q := map(convert,p,':-name');
   4   f := Finance:-ExternalSupport:-InitializeExternalFunction(extfunc);
   5   external := f(data);
   6   m := eval(subs(('__Y__') = ('external'),('__Z__') = op(q),parse("module() option implementation = __Y__; export __Z__; end module")));
   7   for i to nops(p)-1 do
   8     m[q[i]] := parameters[p[i]]
       end do;
   9   m[':-ModulePrint'] := eval(printfunc);
  10   return eval(m,1)
end proc

acer

You'll want to check this reasoning.

It seems to me that you are trying to compute something like a sum of terms of OEIS A001511. As Axel mentioned, finding the greatest common factor of the form 2^k in N! involves a sum of the exponents of the factors of 2 in each even term of the factorial. Hence the summation. And at first glance this is nasty because we want N=2011!, a number too large to take the factorial a second time.

But there is a formula (Daniele Parisse) given for terms in that A001511 sequence, as a function of terms of the OEIS A000129 sequence, namely:

2-A000120(n)+A000120(n-1)

And so it looks to me as if a sum of such A000120 terms will be a telescoping sum. Hoo-rah.

If I've done it right, then procedure `PP` below computes that telescoped sum. You can enjoy checking whether it is OK. It might be off by one, or...

And there in `PP` you can see n! since we wish to find the whole 2^k factor for 2011! as the input to `PP`.

PP:=proc(n) n! - Q(n!/2); end proc:

# Q and QQ both produce OEIS A000120, http://oeis.org/A000120
Q:=proc(n) local w, m, i;
   w := 0; m := n;
   while m > 0 do
     i := m mod 2;
     w := w+i;
     m := (m-i)/2;
   end do;
   w;
 end:

# a (slower?) alternative to QQ (Benoit Cloitre)
QQ:=proc(n)
     log[2](((2*n)!/(n!)^2)/numer(binomial(2*n,n)/4^n));
   end proc:

seq(PP(i),i=2..10);

        1, 4, 22, 116, 716, 5034, 40314, 362874, 3628789

for j from 2 to 10 do
  st:=time():
  if irem((j!)!,2^PP(j)) = 0 and
     irem((j!)!,2^(PP(j)+1)) > 0 then
    print(j,"pass",time()-st);
  else
    print(j, "fail",time()-st);
  end if;
end do:

                         2, "pass", 0.
                         3, "pass", 0.
                         4, "pass", 0.
                         5, "pass", 0.
                         6, "pass", 0.
                         7, "pass", 0.
                        8, "pass", 0.125
                        9, "pass", 2.547
                       10, "pass", 56.125

# PP(2011); # if you want a result to look at (don't worry, it is fast, less than a page)

acer

Using Maple 12.02 or 14.01 (Windows 7, 32bit Maple),

s:=Sockets:-Open("Tomorrowsgaspricetoday.com",80):
Sockets:-Write(s,"GET /rssfeed.php?city=16 HTTP/1.1\nContent-Type: text/plain\nUser-Agent: Maple\nHost: tomorrowsgaspricetoday.com:80\n\n"):
R:="": count:=0:
while count<100 do
  r:=Sockets:-ReadLine(s,30);
  if r=false then break;
  else R:=cat(R,"\n",r);
       count:=count+1;
  end if;
end do:
R;

In Maple 15, (Windows 7, 32bit Maple),

HTTP:-Get("http://tomorrowsgaspricetoday.com/rssfeed.php?city=16");