Carl Love

Carl Love

26401 Reputation

25 Badges

11 years, 147 days
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

Good work so far.

(ii) The concentration is just the amount of salt divided by the volume of the tank. The volume is constant at 100 liters. The amount is your answer to (i).

(iii) "Steady state" just means the limit as time t goes to infinity. The answer to that is probably obvious to you. But you can get it in Maple via limit(sol, t= infinity);

(iv) Part (iv) is to part (iii) as part (ii) is to part (i). Understand?

In fact, Maple's square roots are the principal values, and thus Maple's square roots of positive values are always positive. The situation that you describe arises simply because Maple chooses not to do the simplification automatically. But you tell Maple to do it with combine.

sqrt(6) - sqrt(2)*sqrt(3);
                      (1/2)    (1/2)  (1/2)
                     6      - 2      3     




                   (1/2)  (1/2)  (1/2)  (1/2)
                  2      3      5      7     



Your expression can be entered as sin(x)^3 or (sin^3)(x), but sin^3(x) will not work. Once the expression is entered, you can use the command combine. That might seem like a strange name for this operation, but this is just a special case.

ex:= (sin^3)(x);
                      1            3       
                    - - sin(3 x) + - sin(x)
                      4            4       

I got stuck on the O(n~+ h) algorithm that I was thinking of. I couldn't figure out a way to avoid comparing every pair of intervals before applying Dijkstra's algorithm. Dijkstra's itself is O(n~+ h), n being the number of vertices (intervals in this case) and h being the number of edges (pairs of intervals that overlap in this case). Since every pair of intervals is compared (to see if they overlap), this solution in O(n2). That's okay for 1000 intervals.

This problem seems closely related to other standard problems of combinatorial optimization, such as the knapsack problem. Does anyone here know if it is equivalent to one of the standard problems? Joe Riel, I have a feeling that you may know something about this.

The solution that I was first attempting tried to find the overlapping pairs of intervals by using separate sorted lists of the lower bounds and upper bounds. I tried to use the new permutation option to sort in conjunction with this. Put succinctly, the relevant subproblem is this:

Given a set of intervals, find all pairs that overlap, using sorting rather than checking every pair.

Any ideas? Joe?

So, here's my code:

MinimalSubcover:= module()
option `Written 10 May 2013 by Carl Love.`;
uses GT= GraphTheory;
    #Each interval becomes a vx in a graph. Vxs are adjacent iff their intervals overlap:
    #Let I1 = (a1,b1) and I2 = (a2,b2) be intervals. There is an edge I1 -> I2 iff
    # a1 < a2 < b1 < b2.
    Edges,  #Table of edges
    EdgeCt, #Integer index to the table
    Lower,  #List of lower bounds of intervals
    Upper,  #List of upper bounds of intervals

    #Compute an edge's weight and add it to the table. An edge's weight is the average
    #of the weights of its vxs. A vx's weight is the upper bound of its interval minus
    #the lower bound. The special vertices "Start" and "End" have weight 0.
    AddEdge:= proc(e::list)
    local w; # edge weight
        if e[1]="Start" then  w:= Upper[e[2]] - Lower[e[2]]
        elif e[2]="End" then  w:= Upper[e[1]] - Lower[e[1]]
        else  w:= Upper[e[2]]-Lower[e[2]] + Upper[e[1]]-Lower[e[1]]
        end if;
        EdgeCt:= EdgeCt+1;
        Edges[EdgeCt]:= [e, w/2];
    end proc,

    ModuleApply:= proc(Cover::{list,set}, ClosedInterval::specfunc(numeric, RealRange))
        A:= op(1, ClosedInterval),
        B:= op(2, ClosedInterval),
        # Remove duplicates.
        # Force to a list.
        # Remove any intervals that do not intersect the target interval.
        Subcover:= remove(OI-> op([1,1], OI) >= B or op([2,1], OI) <= A, [{Cover[]}[]]),
        N:= nops(Subcover),
        k, k1, k2, Soln
         Edges:= table();
         EdgeCt:= 0;
         Lower:= map(x-> op([1,1], x), Subcover);
        Upper:= map(x-> op([2,1], x), Subcover);

        #Intervals with lower bounds less than the lower bound of the target interval get
        #connected to special vx "Start". Those with upper bounds greater than the upper
        #bound of the target interval get connected to "End".
        for k to N do
            if Lower[k] < A then  AddEdge(["Start", k])  fi;
            if Upper[k] > B then  AddEdge([k, "End"])  fi
        end do;

#Check each pair of intervals for overlap.
        for k in combinat:-choose(N,2) do
            (k1,k2):= k[];
            if Lower[k1] < Lower[k2] and Lower[k2] < Upper[k1] and Upper[k1] < Upper[k2] then
            elif Lower[k2] < Lower[k1] and Lower[k1] < Upper[k2] and Upper[k2] < Upper[k1] then
        end do;

        #Compute shortest path through graph.
        Soln:= GT:-DijkstrasAlgorithm(GT:-Graph({entries(Edges, nolist)}), "Start", "End");
        Soln[2], Subcover[Soln[1][2..-2]]
infinity, []
end try
    end proc
end module:
Sorry, I don't know how to get rid of the extra blank lines below. They do not appear
 in the editor when I am creating the post.

G:= RandomTools:-Generate(makeproc, float(range= 0..10^4, digits= 4)):

L:= ['RealRange('Open(G())' $ 2)' $ 2*10^3]:

M:= remove(`=`, L, BottomProp):



st:= time():

MinimalSubcover(M, RealRange(1., 1000.));

1031.10940350000, [RealRange(Open(0.5965e-3), Open(99.83)), RealRange(Open(99.72), Open(1031.))]





It's the same situation as with the Koch snowflake: The initial square has no left turns: "FRFRFRFR" (or, as one may say in Maple, cat("FR"$4)). We replace every F by the recursion to get

cat(cat(Box(n), "R") $ 4);

You are correct that the U shape is the recursive "unit": "FRFLFLFRF". Note that there's no "R" on the end of that. Now replace every F in that with recursion.

For a non-polynomial equation, fsolve gives just one solution. One of the best ways to get all the roots of an analytic function in an interval of reals is to use RootFinding:-NextZero:

r:= -5:
for k while r < 5 do
   r:= RootFinding:-NextZero(x-> sin(Pi*x^2/2), r);
   Sols[k]:= r
end do:
k:= k-2; # Number of solutions
seq(Sols[i], i= 1..k);
-4.89897948556636, -4.69041575982343, -4.47213595499958,
  -4.24264068711928, -4.00000000000000, -3.74165738677394,
  -3.46410161513776, -3.16227766016838, -2.82842712474619,
  -2.44948974278318, -2.00000000000000, -1.41421356237310, 0.,
  1.41421356237310, 2.00000000000000, 2.44948974278318,
  2.82842712474619, 3.16227766016838, 3.46410161513776,
  3.74165738677394, 4.00000000000000, 4.24264068711928,
  4.47213595499958, 4.69041575982343, 4.89897948556636

I don't know the exact reason that your example doesn't show any benefit. I do duplicate your results on my machine. I also randomized the list order and still got no better results from Threads. It may be that the task (computing the 10th power of an integer less than 10^6) is too trivial. Here's an example where there is a clear benefit to using Threads:

L:= RandomTools:-Generate(list(posint, 2^16)):
CodeTools:-Usage(map(nextprime, L)):

memory used=2.25GiB, alloc change=2.00MiB, cpu time=21.17s, real time=21.20s

L:= RandomTools:-Generate(list(posint, 2^16)):
CodeTools:-Usage(Threads:-Map(nextprime, L)):

memory used=2.25GiB, alloc change=71.34MiB, cpu time=52.78s, real time=8.38s

Done on an i7 3630QM 4x2.4 GHz

Just change the = between the seqs to =~.

Or you can make a seq of equations:

seq(seq(p[j,c]= 1/(1+exp(-(mu+cat(tau,j)+cat(eta,c)+mix[j,c]))), j= 2..3), c= 1..3);

which is more efficient than the =~ option.

You wrote: map2(myfunc, phi, p1, {CH||(1..5)});

But that needs to be map[3], since you are replacing the 3rd argument of myfunc with a set/list. Note that map2 is an abbreviation for map[2], but that there is no such abbreviation for higher indices.

You wrote: Is it not by default for Maple to 'adapt' to computer and use multi thread then?

It can't be done automatically (yet); each piece of code needs to be inspected for multithreading possibility. If the computation of later iterations depends on the results from earlier iterations, then multithreading isn't possible. I believe that reviewing the whole library for multithreading possibilities is an ongoing process at Maplesoft. For example, Maple 17 introduced multithreaded garbage collection, which is done automatically (no input is required from the user). The user can control this via kernelopts(gcmaxthreads).

You wrote: ANS1:=w1*a+w2*b+(1-w1-w2)*b

I'll assume that you meant for the final b to be c; otherwise, the cs are totally unused.

With that assumption, all of the above code can be compressed to:

ANS:= < w1, w2, 1-w1-w2 > ^%T .

     Matrix([Threads:-Seq](Threads:-Map[3](myfunc, phi, p, [CH||(1..5)]), p= [p||(1..3)]))


returning a 5-vector for ANS.

First, there is no need to wrap z(.., 1) with Vector(...); it is automatically converted to a column vector by the indexing.

Since Statistics:-Tally is prepared to take vectors (or lists), you could use it, send the output to a table, and index the table for the element that you want. Sounds complicated, but all it is is

Statistics:-Tally(z(.., 1), output= table)[1];


The reason that you got an error is that you forgot to include phi, your first argument to myfunc, in your map2 call. So, it should be

map[3](myfunc, phi, p1, {CH||(1..5)});

To make that run multi-core, replace map with Threads:-Map.

Taking specific values of n, we get equations in a, b, c, d. We solve those equations.

T:= n-> add(binomial(n,j)^4, j= 0..n):
EQN:= n-> (n+1)^3*T(n+1) = (2*n+1)*(a*n^2+a*n+b)*T(n) + n*(c*n^2+d)*T(n-1):
solve({seq}(EQN(k), k= 1..4), {a, b, c, d});

                 {a = 6, b = 2, c = 64, d = -4}


So the solution is not in the ranges -30..30.

theta:= (n::{odd,posint})-> 
          2*sin((n+1)*beta/2)/(n+1)/beta +
          `if`(n=1, 1,

@jenniferchloe In order to end up with a closed curve, it needs to be a closed curve (a polygon) at the initial level. And you can't do that with just any angle. To make an initial equilateral triangle, you need to make turns of 120 degrees. If the angle is set at 60, each of those turns is RR. So, to do the spiky 80-degree thing all the way around,

  1. set the angle to 40
  2. make the initial triangle with RRR turns (3*40 = 120)
  3. replace every single turn in the recursion with a double turn (2*40 = 80).


First 358 359 360 361 362 363 364 Last Page 360 of 381