Carl Love

## 25 Badges

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

## Concentration is just amount / volume...

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?

## combine...

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     combine(%);                               0sqrt(2)*sqrt(3)*sqrt(5)*sqrt(7);                   (1/2)  (1/2)  (1/2)  (1/2)                  2      3      5      7     combine(%);                               (1/2)                            210     `

## combine...

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);                                  3                            sin(x) combine(%);                      1            3                           - - sin(3 x) + - sin(x)                      4            4       `

## Solved: Minimal subcover with Dijkstra's...

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;local    #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))    local        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                AddEdge([k1,k2])            elif Lower[k2] < Lower[k1] and Lower[k1] < Upper[k2] and Upper[k2] < Upper[k1] then                AddEdge([k2,k1])            fi        end do;        try                    #Compute shortest path through graph.             Soln:= GT:-DijkstrasAlgorithm(GT:-Graph({entries(Edges, nolist)}), "Start", "End");             Soln[2], Subcover[Soln[1][2..-2]]        catch:             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):

>

nops(M);

>

st:= time():

>

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

>

time()-st;

>

Download Subcover1.mw```

## The initial square is a different case...

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.

## RootFinding:-NextZero...

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]:= rend do:k:= k-2; # Number of solutions                               25seq(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`

## Try something more complex (not just lar...

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:

`restart;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.20srestart;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.38sDone on an i7 3630QM 4x2.4 GHz`

## Change = to =~...

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.

## map2, map[2], map[3]...

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).

## The whole code, compressed and multi-thr...

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.

## By Statistics:-Tally...

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];

1

## You forgot about phi...

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.

## Unique solution...

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}`

Download recurrence.mw

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

## n::{odd,posint}, `if`(...)...

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

## Need to count the turns and use correct ...

@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
﻿