Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

To load the contents of an audio file into a numeric Vector, use the command AudioTools:-Read (see ?AudioTools,Read).

If your variables were actually indexed as j[4] rather than as j_4, then the Answer of AmirHosein would suffice. To work with the actual j_4, etc., you need to use type suffixed. For example,

M:= mul(j_||k^rand(0..2)(), k= 0..15);

PrefLen:= length(j_); #Prefix length of `j_`

J:= indets(M, suffixed(j_, nonnegint));

add(parse(substring(j, PrefLen+1..-1))*degree(M,j), j= J);

     67

You can save them anywhere. I put them all in one directory (or folder) called "C:/Maple_packages". But you can put them in separate directories (or folders). To activate a package, all you have to do is append its directory name to the predefined global sequence libname:

libname:= libname, "C:/Maple_packages";

Now you can access the commands and the help pages of the packages just like you do with any regular Maple package.

Make sure that you append to libname rather than overwrite it. DON'T do

libname:= "C:/Maple_packages"; #DON'T DO THIS!!!

And you probably shouldn't do this either:

libname:= "C:/Maple_packages", libname; #Don't do this unless you really know what you're doing.

Once you're satisfied that the package is something that you'd like to use regularly, you can put the command which appends to libname into an initialization file, and then you can pretty much forget about it: Those add-on packages will be just like regular Maple packages.

(The Reply was updated for clarity due to Preben's Reply below.)

If RR1 and RR2 are two RealRanges, then RR1 being a subset of RR2 is equivalent to

is('x' in RR2) assuming 'x' in RR1;

Of course, the quotes can be omitted if you know that x is unassigned.

Operator in can be replaced with operator ::, which saves typing four spaces:

is(x::RR2) assuming x::RR1;

(Don't take this to mean that in can usually be replaced by ::. It usually can't be.)

Yes, it's possible to get the entire history of optimizations done with the Optimization package (procedures Maximize, Minimize, LPSolve, and NLPSolve) by using the procedural form of input (i.e., the objective and constraints are specified by procedures rather than expressions). See ?Optimization,General,OperatorForm. A crude way to do this would be to print out the values of the arguments inside the procedures. A more refined way is to increment an index and save those values to a table or Array. I'll post an example later.

My first inclination was to use a try ... catch statement as suggested by Robert Lopez, but then I saw an easier way for your situation. The break statement exits a do loop. For your code, setting h[i+1] to 0 was equivalent to exiting the loop. Notice how the if statement first checks the type of h[i+1]. This makes it safe to compare it with 0.

restart;
p1[0]:=
     .989347189582843*x^2 - 0.139423979061219e-1*x -
     1.82559474469870e8*x^15 + 1.30761381361453e8*x^16 -
     6.88520063191821e7*x^17 + 2.51079317463498e7*x^18 -
     5.66094206949155e6*x^19 + 5.94129446612678e5*x^20 -
     6812.74182685426*x^5 + 59230.0931084044*x^6  -
     3.83520584559500e5*x^7 + 1.90126822307036e6*x^8  -
     7.34991883857609e6*x^9 + 2.24203561757434e7*x^10 -
     5.43284775909785e7*x^11 + 1.04806113793011e8*x^12 -
     1.60600324339222e8*x^13 + 1.94090536353833e8*x^14 +
     559.557918804679*x^4 - 30.6576714427729*x^3  -
     3.93727537464007e-15
:

h[1]:= -1;
for i to 149 do
     p1[i]:= p1[0]-i^2-i:
     h[i+1]:= fsolve(p1[i]-0.312e-1, x, -.2 .. 0);
     if h[i+1]::numeric implies h[i+1] >= 0 then  break  end if;
     print(i, h[i])
end do;

Since you defined h as a procedure (because of the arrow) rather than as an expression, you wouldn't use eval to evaluate it. Just use h(10).

There's a bug in your h. You need to change cos^4*(x) to cos(x)^4.

If you want an actual decimal evauluation, use h(10.); h(10) essentially gives just a substitution without any significant computation in this case.

Most of the plot commands will try to use hardware float arithmetic to evaluate your expression, but they will fail because of the extreme exponents. To counteract this, use

UseHardwareFloats:= false:

Now assign your expression to F (or whatever name you like), and plot like this:

plots:-complexplot(F, del= -1e-13..1e-13);

and

plot([Re,Im](F), del= -1e-13..1e-13);

There are many other options for visualizing a complex function.

Rouben: Don't worry: This isn't a homework problem. Brian is a long-time (many years) poster of optimization puzzle-type problems.

Brian/general readers: There's a trick for converting objective functions with absolute values---but which are otherwise linear---into linear objective functions by adding variables. Unfortunately I don't remember the trick. It should be implemented automatically in Optimization:-LPSolve, but it's not. If it were, I'd solve this as an Integer Linear Program (ILP). Package DirectSearch solves INLPs, but package Optimization doesn't. It turns out that the integer solution can be obtained by the non-integer method of Optimization:-NLPSolve.

 

Problem setup.

restart:

Data:= <name, x, y, pop;
        A, 3, 7, 5;  B, 6, 5, 2;  C, 5, 3, 2;
        D, 1, 2, 1;  E, 7, 1, 3              >;

Data := Matrix(6, 4, {(1, 1) = name, (1, 2) = x, (1, 3) = y, (1, 4) = pop, (2, 1) = A, (2, 2) = 3, (2, 3) = 7, (2, 4) = 5, (3, 1) = B, (3, 2) = 6, (3, 3) = 5, (3, 4) = 2, (4, 1) = C, (4, 2) = 5, (4, 3) = 3, (4, 4) = 2, (5, 1) = D, (5, 2) = 1, (5, 3) = 2, (5, 4) = 1, (6, 1) = E, (6, 2) = 7, (6, 3) = 1, (6, 4) = 3})

n:= op([1,1], Data) - 1;

5

City:= Vector(
     n,
     k-> Record(convert(Data[1,..]=~ Data[k+1,..], list)[])
);

City := Vector(5, {(1) = Record(name = A, x = 3, y = 7, pop = 5), (2) = Record(name = B, x = 6, y = 5, pop = 2), (3) = Record(name = C, x = 5, y = 3, pop = 2), (4) = Record(name = D, x = 1, y = 2, pop = 1), (5) = Record(name = E, x = 7, y = 1, pop = 3)})

Dist:= (P1,P2)-> LinearAlgebra:-Norm(P1-P2, 1): #1 is taxi metric.

P:= <x,y>:

Obj:= add(Dist(<City[k]:-x, City[k]:-y>, P)*City[k]:-pop, k= 1..n);   

5*abs(x-3)+5*abs(y-7)+2*abs(x-6)+2*abs(y-5)+2*abs(x-5)+2*abs(y-3)+abs(x-1)+abs(y-2)+3*abs(x-7)+3*abs(y-1)

 

 

Solution method #1:

infolevel[Optimization]:= 2:

Optimization:-NLPSolve(Obj, method= nonlinearsimplex, evaluationlimit= 2^13);

NLPSolve: calling NLP solver

SolveUnconstrainedNM: using method=nonlinearsimplex
SolveUnconstrainedNM: number of problem variables 2
SolveUnconstrainedNM: trying evalhf mode
E04CCA: number of function evaluations 251

[51.0000000000183604, [x = HFloat(4.999999999989731), y = HFloat(4.999999999991912)]]

 

 

Solution method #2:

#Brute-force integer method requires fewer evaluations!

Min:= infinity:
for X in $(min..max)(Data[2..,2]) do
     for Y in $(min..max)(Data[2..,3]) do
          dist:= eval(Obj, [x= X, y= Y]);
          if dist < Min then               
               (Min,Xmin,Ymin):= (dist,X,Y)
          end if
     end do
end do:

Min,Xmin,Ymin;

51, 5, 5

 

 

Solution method #3:

DirectSearch:-GlobalSearch(Obj, assume= integer);

Matrix(1, 3, {(1, 1) = 51., (1, 2) = [x = 5, y = 5], (1, 3) = 98})

Hmmm. The number of objective function evaluations, 98, is exactly twice the number of lattice points used in the brute-force technique. That's weird and worthy of further investigation.

 

 

Solution method #4?

This simple method is implied by the text in your post. Does it always work? 

 

Xmed:= Statistics:-Median(Data[2..,2], weights= Data[2..,4]);

HFloat(5.0)

Ymed:= Statistics:-Median(Data[2..,3], weights= Data[2..,4]);

HFloat(5.0)

 

 

Download Brian_taxi_metric.mw

 

 A multi-argument function call such as arctan(y,x) could appear in a denominator.

That's a frequently requested feature, but it's not available now.

Here are some alternatives:

  1. Increase the size of the plot.
  2. Decrease the font size of the tickmarks (option axesfont).
  3. Remove the leading zeroes from the tickmarks.
  4. Remove some tickmarks.

Let me know if you need help with any of these.

 

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.

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