Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 358 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Since the exercises were already due, I thought that there'd be no harm in posting solutions now. I urge you to experiment with various changes to this code. My level of explanation below assumes that you have some experience programming in some other language, although I don't know if that's true. Maybe someone else will add more explanation. Or, feel free to ask specific questions. I also assume that you'll read the help pages for all the commands that I use.

Write a Maple proc L that takes as input two lists and returns as output a set containing the terms they have in common (not necessarily in the same position). For example, L([1, 3, x, 5],[2, x, 5, 2, 1, 8]) should return [1, x, 5].

There's a contradiction in the problem statement. He says that the output should be a set, but he shows it being a list. I'll just assume that it should be a set.

L:= (A::list, B::list)-> {A[]} intersect {B[]};

or

L:= proc(A::list, B::list)  {A[]} intersect {B[]}  end proc;

Note: If A is a set or list, A[] returns the sequence of elements; so {A[]} converts a list to a set.

This procedure can be easily generalized to take an arbitrary number of sets/lists as arguments:

L:= (L::seq({set,list}))-> `intersect`(seq({_l[]}, _l= L));

 

Write a Maple program that implements the Euclidean algorithm to compute the gcd(x, y) for any natural numbers x, y (except x = y = 0).

I saw no reason to treat (0,0) as a separate case. I just have it return 0, which is what Maple's gcd and igcd both do.

GCD:= proc(X::integer, Y::integer)
local r, x, y;
     (x,y):= (max,min)(abs~([X,Y])[]);
     while y <> 0 do
          r:= irem(x,y);
          (x,y):= (y,r)
     end do;
     x
end
proc;

If you really insist on not allowing (0,0) as input, change it to

GCD:= proc(X::integer, Y::integer)
local r, x, y;
     (x,y):= (max,min)(abs~([X,Y])[]);
     if x = 0 then  error "Invalid input: Both arguments can't be 0"  end if;  
     while y <> 0 do
          r:= irem(x,y);
          (x,y):= (y,r)
     end do;
     x
end proc;

or, better yet, replace the error statement with infinity or undefined.

Notes:

  1. abs~([X,Y]) is short for [abs(X), abs(Y)].
  2. (max,min)(x,y) is short for max(x,y), min(x,y).
  3. (x,y):= (y,r) is short for x:= y; y:= r.
  4. In lieu of a return statement, the last expression evaluated in a procedure will be its return value. In this case, that's the x before the end proc.

For a two-line version (as promised), make an inner procedure which acts recursively after the outer procedure has sorted the absolute values of the arguments:

GCD:= (X::integer, Y::integer)->
    ((x,y)-> `if`(y=0, x, thisproc(y, irem(x,y))))((max,min)(abs~([X,Y])[]));

 

Write a Maple program to plot the Pythagorean spiral. The program should take a single positive integer n as an argument and plot the spiral until the side of length n is reached.

Since the sides have lengths that are square roots of integers, I'm surprised that he didn't say "until the side of length sqrt(n) is reached." Anyway, that doesn't add any complexity other than that we need to loop up to n^2 rather than n.

Hopefully you've read the Wikipedia page so you know that the nth radial side is sqrt(n) long and the nth central angle is arctan(1/sqrt(n)).

My "long" version is two procedures. The first loads a two-column matrix with the polar coordinates of the points; the second plots them.

PShf:= proc(n::posint, Pts::Matrix)
local k, r, theta:= 0;
     Pts[1,1]:= 1;
     for k from 2 to n do
          r:= evalf(sqrt(k-1));
          theta:= theta + evalf(arctan(1/r));
          Pts[k,1]:= r;  Pts[k,2]:= theta
     end do
end proc:

PythagoreanSpiral:= proc(n::posint)
local k, Pts:= Matrix((n^2,2), datatype= float[8]), O:= <0|0>;
     evalhf(PShf(n^2, Pts));
     plot([Pts, seq(<Pts[k,..], O>, k= 1..n^2)], coords= polar, scaling= constrained)
end proc:

Notes:

  1. The angle brackets < > are Matrix/Vector constructors; so <0|0> is a row vector of two zeroes (in this case it represents the origin on the plot).
  2. A dimension of a Matrix/Vector/Array specifed by simply .. indicates that the whole of the dimension is to be extracted; so Pts[k,..] extracts the kth point as a row vector.
  3. <Pts[k,..], O> stacks the two row vectors into a 2x2 Matrix.
  4. plot interprets a two-column Matrix as a sequence of points to be plotted and connected.
  5. datatype= float[8] declares that the Matrix will only contain eight-byte double precision values. These are usually the most efficient numbers to work with.
  6. scaling= constrained means that the unit measure on the x- and y-axes will be visually the same (see ?plot,options).
  7. The plot itself is the return value of the procedure.

My three-line version:

PythagoreanSpiral:= (n::posint)->
     (Pts-> plot([Pts[], map2(op, 1, Pts)], coords= polar, scaling= constrained))
          ([seq([[sqrt(_k), evalhf(add(arctan(1/sqrt(_j)), _j= 1.._k-1))], [0,0]], _k= 1..n^2)]);

The evalhf just makes them faster. You can remove the evalhf in either case and they'll still run. If you're going to be plotting any significantly detailed fractals, then you'll need to learn evalhf, which is a tricky command to use correctly.

 

The command fnormal doesn't automatically map over Matrices. Change your command

fnormal(C12f, Digits, 1e-6);

to

fnormal~(C12f, Digits, 1e-6);

(There seems to be no rhyme nor reason to which commands automatically map themselves over which container objects. The command fnormal automatically maps overs lists and sets, but not tables or rtables.)

There is a command essentially for this purpose: plottools:-transform. Here's an example of its use. I use  Kitonum's last example---his plane mapping applied to a unit square---and I show the first three iterations.

#Unit square parametrized on 0..4:
Sq:= t-> piecewise(t < 1, [t,0], t < 2, [0,t-1], t < 3, [t-2,1], [1,t-3]):
#Split it into coordinate functions (required by `plot`):
Sq||(x,y):= seq(subs(_k= k, t-> Sq(t)[_k]), k= 1..2):
#plottools:-transform requires that F return a list (in square brackets):
F:= (x,y)-> [x^2/2 - 3/2*y^3, -x^3/2 + y^2/3]:
T:= plottools:-transform(F):
P:= plot([Sqx, Sqy, 0..4], thickness= 3):
plots:-display(
     Matrix(
         (2,2),
         [seq(plots:-display((T@@k)(P), title= cat("k = ", k)), k= 0..3)]
     ), scaling= constrained
);

Sorry, MaplePrimes can't display array plots from the Standard GUI. You'll need to copy-and-paste the code to see it.

How's it possible that you've been on MaplePrimes for nearly three years and have posted 366 Questions and yet you don't know the answer to this very simple and fundamental question?

p:= x^2*y + x*y + x:
[op(p)];

Please read the help ?op---it's one of the most commonly used low-level commands.

The following could be made faster by pregenerating arrays of random numbers (like you did in your code), but this is a good exposition of the Maple animation technique.

2D animation:

U:= RandomTools:-Generate(float(method= uniform, range= 0..1), makeproc):
frames:= 2^5:
points:= 2^9:
plots:-display(
     ['plot(x-> abs(x)*U(), -5..5, adaptive= false, numpoints= points)' $ frames],
     insequence, thickness= 0
);

3D animation:

frames:= 2^4:
points:= 2^7:
plots:-display(
     ['plot3d((x,y)-> abs(x)*U()+abs(y)^2*U(), -5..5, -5..5, grid= [points$2])' $ frames],
     insequence, shading= zhue, thickness= 0, style= wireframe, transparency= 0.5,
     axes= frame
);

I think that the 3D one may be too bulky to animate here in the post, but you will see it animated in a worksheet. You should use the toolbar controls to slow the Frames Per Second (FPS) to 1 to give the rendering a chance to keep up.

Here is the complete and simplified code. This procedure will print the complete contents of a module---all locals and all exports---and apply itself recursively to all submodules. I was able to simplify from the earlier code by eliminating the need for unwith and parse.

ModulePrint:= proc(M::`module`)
uses ST= StringTools;
local
     N,
     SIV:= interface(verboseproc= 2),
     SKO:= kernelopts(opaquemodules= false),
     EVs:= eval~(op~([1,3], eval(M)))
;
     for N in ST:-Split(sprintf(" %q", op~([1,3], eval(M))[]), ",") =~ EVs do
          printf("%s", lhs(N));
          print(rhs(N));
          if rhs(N)::'`module`' then  thisproc(rhs(N))  end if;
     end do;
     kernelopts(opaquemodules= SKO);
     interface(verboseproc= SIV);
     return
end proc:

combinat:-permute([0$3,1$3], 3);

     [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]]

I'm guessing that the above list is what you really want. If you actually want the output as you presented it, then add (as Kitonum did)

(cat@op)~(%)[];

     000, 001, 010, 011, 100, 101, 110, 111

Here's a procedure for it. It's a bit quick-and-dirty because I'm in a hurry right now. Ideally, it'd do exports as well as locals, but the specific module that you mentioned doesn't have any exports. Also, ideally, it'd apply itself recursively to submodules.

The printout of the Maplet code is unfortunately not indented. This would be difficult to achieve, but I'll think about it.

restart:

ModulePrint:= proc(M::`module`)
uses ST= StringTools;
local
     Ns:= ST:-Split(sprintf("%q", op(3, eval(M))), ","),
     N
;
     interface(verboseproc= 2);
     kernelopts(opaquemodules= false);
     for N in Ns do
          lprint(N);
          parse(sprintf("print(eval(%s))", N), 'statement')
     end do
end proc:

ModulePrint(Student:-Calculus1:-DiffTutor);

Edit: There is a bug in the above procedure. Please use the corrected code in the third Reply below.

If your evaluations of the parameters are rational numbers, then an exact solution can be obtained instantaneously, as in

Vars:= {A||(1..8)}:
Rand:= rand(1..9):
Params:= indets(sys1, name) minus Vars:
Vals:= zip(`/`, '['Rand()'$nops(Params)]'$2):
sys2:= eval(sys1, Params=~ Vals):
<solve(sys2, Vars)[]>;

Even exact rational solutions with thousands of digits can be obtained instantaneously. Maple/GMP is very fast with rational arithmetic (state of the art, perhaps?).

If you have any procedure, and you want to expand its domain to types that it wasn't designed for, or you want to change its behavior for certain types of input, then you can use overload. In this case, we want to extend the definition of numelems to include input of type name.

unprotect(numelems);
numelems:= overload([
     proc(C::name &under (C-> eval(C))) option overload; 0 end proc,
     numelems
]):

numelems(a);
     0

numelems([1,2,3]);

     3

numelems(a+b);

Error, invalid input: numelems expects its 1st argument, t, to be of type indexable, but received a+b

(A curiosity about the above code: I was forced to use C-> eval(C) rather than the simpler eval. The latter tries to use two-argument eval, as seen in the following:

type(a, name &under eval);

Error, (in type/&under) invalid input: eval received NULL, which is not valid for its 2nd argument, eqns

This is a very curious error message because I thought that a procedure parameter could never be passed NULL! But what I'm more curious about is why it chooses to use two-argument eval. Here's a simpler example of the same thing:

eval(a, NULL);

Error, invalid input: eval received NULL, which is not valid for its 2nd argument, eqns)

Your entire loop needs to be in one execution group. To achieve this, use the editing command Join Execution Groups. You can find this by going to the Edit menu, then Split or Join, then Join Execution Groups. The shortcut key for this is F4. So, the entire operation can be done by placing your cursor at the top of the worksheet and repeatedly pressing F4. The cursor won't move, but if you look closely at the execution group boundaries at the left edge of the worksheet, you'll see them change.

After having merged the execution groups, add the loop header to the top and end do; to the bottom. If any of your original execution groups contained text rather than code, it doesn't matter: The parser will treat that text as comments. So, this process also gives you a way to include formatted text comments in your code.

This is the same error message as you got with your bubble sort procedure a few days ago. Like I told you, you can't ordinarily assign to a formal parameter, such as n. In other words, n shouldn't appear on the left side of :=. If you need to do something like that, first copy n to a local variable, as Markiyan has copied it to k. Also, I told you not to use the print statement as a substitute for a procedure's return value.

Here's my Collatz procedure:

Collatz:= proc(N::posint)
local n:= N, count;
     for count while n <> 1 do
          n:= `if`(n::even, iquo(n,2), 3*n+1)
     end do;
     count
end proc:

I find binary infix operators very useful. In this case, I agree with the use of the notation of that "other CAS", although I usually abhor it. The following code allows for the use of notation very similar to what you posted. You just need to change the angle symbol to &<.

phasor:= module()
option package;
export
     `+`:= proc(Z1::`&<`(realcons,realcons), Z::seq(`&<`(realcons,realcons)))
     option overload;
     local z, r:= add(op(1,z)*exp(op(2,z)*Pi/180*I), z= args);
          evalf(abs(r)) &< evalf(argument(r)*180/Pi)
     end proc,

     `-`:= proc(Z::`&<`(realcons,realcons))
     option overload;
          (-1*op(1,Z)) &< op(2,Z)
     end proc
;
end module:

with(phasor):

4&<45 + 5&<30;

This allows for sums and differences with an arbitrary number of phasors. With a little more work, I could get the output to print without the & and the extra parentheses.

After you fix the error pointed out by Robert Lopez, you should get an error about pi. In Maple, this must be capitalized (Pi) if you mean the famous mathematical constant.

Here's an example of what you want. It involves using the keyword unevaluated function typeset in the plot command (see ?plot,typesetting). I tested with 99 frames of triplets of random numbers, and there's no jitter. All that you need from this are the procedures P3 and T3; the rest is just code I wrote to test them.

P3:= (x::realcons)-> sprintf("= %5.3f   ", x):
frames:= 99:
L||(1..3):= 'RandomTools:-Generate(list(float(method= uniform), frames))' $ 3:
T3:= (f::seq(realcons))-> typeset(seq([theta[_k], P3(f[_k])][], _k= 1..nargs)):
plots:-display(
     [seq(
          plot(
               0,
               title= T3(seq(L||_k[_j], _k= 1..3)),
               titlefont= [TIMES,ROMAN,14]
          ), _j= 1..frames
     )],
     insequence
);


First 247 248 249 250 251 252 253 Last Page 249 of 395