Carl Love

Carl Love

27656 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

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

1013

st:= time():

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

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

time()-st;

6.000

 

Download Subcover1.mw

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
                               25
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:

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.20s

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

                                                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}

Download recurrence.mw

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

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

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

 

The equation can be entered like this:

w(x,y) = 4*F/Pi^4/E[2]/I/a/b*

Sum(Sum(sin(m*Pi*xi/a)*sin(n*Pi*eta/b)*sin(m*Pi*x/a)*sin(n*Pi*y/b)/((m/a)^2+(n/b)^2)^2,

   n= 1..infinity), m= 1..infinity

);

"Solving" it is another thing. It's not exactly clear what "solve" means in this case.

Download BigSum.mw

Try this:

plot(sin(x), labels=[typeset(convert(a, `local`)[0]), "some text"]);

To address your first issue: Maple's int does not supply a constant of integration. If you need one, you can add it yourself, or you can use dsolve without initial conditions:

eq:= diff(v(x), x, x) = P*x/E/Iz;
                      d  / d      \   P x
                     --- |--- v(x)| = ----
                      dx \ dx     /   E Iz
dsolve(eq);
                             3               
                          P x                
                  v(x) = ------ + _C1 x + _C2
                         6 E Iz              

I noticed that you changed the initial condition on the derivative between your solution with int and your attempted solution with dsolve. You should specify the condition like this: D(v)(0) = -P*L^2/E/Iz/2. (If you need to specify a condition for a second derivative, then use (D@@2)(v)(0).) Putting it all together, we have

bcs:= v(L)=0, D(v)(0) = -P*L^2/E/Iz/2;
                                            2
                                         P L  
                  v(L) = 0, D(v)(0) = - ------
                                        2 E Iz
dsolve({eq, bcs});
                           3       2         3
                        P x     P L  x    P L  
                v(x) = ------ - ------ + ------
                       6 E Iz   2 E Iz   3 E Iz

 

 

First 369 370 371 372 373 374 375 Last Page 371 of 392