Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

For longer sections of code within a worksheet, use Code edit Regions. These are available from the Insert menu. They do allow the tab key. They offer syntax-based color coding and highighting and bracket matching.

One drawback: In M17, Code Edit Regions do not place the cursor near a syntax error when you make one; you're on your own for finding syntax errors. This is a major drawback that M16 did not have.

To excecute a Code Edit Region, use Ctrl+E. The option is also accessible from the context menu.

Here is my solution to the 4-door Monty Hall problem. I managed to simplify the code a lot compared to my 3-door solution.

# Proc to make a random selection from a set.
Rand:= (S::set)-> S[(rand() mod nops(S)) + 1]:

# Simulation for the 4-door Monty Hall problem. N is the number
# of trials. Returns the number of trials the player wins.
MontyHall4:= proc(
     N::posint,
     {switch1st::truefalse:= false},
     {switch2nd::truefalse:= false},
     $
)
local
     k, wins:= 0, Doors,
     PrizeDoor, #The door with the prize behind it.
     PlayersDoor,  # The door the player picks.
     AllDoors:= {$1..4}   
;
     to N do
          Doors:= AllDoors;
          PrizeDoor:= Rand(Doors);
          PlayersDoor:= Rand(Doors);
          # Monty opens a door, which is thus removed from the set
          # of available doors.
          Doors:= Doors minus {Rand(Doors minus {PrizeDoor,PlayersDoor})};       
          if switch1st then
               PlayersDoor:= Rand(Doors minus {PlayersDoor})
          end if;
          # Monty opens another door:
          Doors:= Doors minus {Rand(Doors minus {PrizeDoor,PlayersDoor})};
          if switch2nd then
               PlayersDoor:= Rand(Doors minus {PlayersDoor})
          end if;
          if PlayersDoor = PrizeDoor then  wins:= wins+1  end if
     end do;
     wins
end proc:     

MontyHall4(100000);
                             25092
MontyHall4(100000, switch1st);
                             37562
MontyHall4(100000, switch1st, switch2nd);
                             62565
MontyHall4(100000, switch2nd);
                             74751

If the blocks of a partition cannot be made to have equal products, then we seek to minimize the "variation" (defined below). I will use Brian's example set

S:= {3, 4, 5, 6, 8, 9, 28, 30, 35}.

First, I compute the ideal target product for a block---the cube root of the product of the elements of S. Call this p. Then for each block b, I measure its deviation from this ideal as abs(ln(`*`(b[]) - ln(p)). I use ln because I think that being a factor of x too large should be considered the same deviation as being a factor of x too small. For each partition P, I measure the variation as the sum of the deviations of its blocks.

The partitions are generated with Joe's Iterator package. You'll need to download this from the Maple Applications Center if you don't already have it.

restart:
S:= {3, 4, 5, 6, 8, 9, 28, 30, 35}:
AllP:= [seq(P, P= Iterator:-SetPartitions(S, [[3,3]], compile= false))]:
lnp:= evalf(ln(`*`(S[])^(1/3))):
Var:= proc(P::list(set))
local r:= evalf(`+`(map(b-> abs(ln(`*`(b[]))-lnp), P)[]));
end proc:
sort(AllP, (x,y)-> Var(x) < Var(y))[1];
              [{3, 9, 35}, {4, 8, 28}, {5, 6, 30}]
---which agrees with Brian's result.

Rather than using sort, I wanted to use Joe's trick of setting the attribute of a float and then using min. But this does not work:

Var:= proc(P::list(set))
local r:= evalf(`+`(map(b-> abs(ln(`*`(b[]))-lnp), P)[]));
     setattribute(r,P)
end proc:

The procedure returns garbage which displays as a question mark. I am hoping that Joe can tell me why.

@fluff

Forget about the identical(0,1,2) thing. It won't work, because I didn't realize that you allowed the possibility of the player not switching in the first round and then switching in the second. You need two switch parameters:

MontyHall:= proc(N::posint, {switch1st::truefalse:= false}, {switch2nd::truefalse:= false})

Then is the code you'll have

if switch1st then
. . .

if switch2nd then
. . .

Also, you still need to correct rendomize to randomize, as pointed out by Acer.

It looks like the command defining params is not being executed. Try placing the restart in its own execution group.

Your newsys is a list. In order to use it with union, it must be a set. In your dsolve command change newsys to {newsys[]}.


A:= Matrix([[.7, .3, .3], [.2, .6, .1], [.1, .1, .6]]);

A := Matrix(3, 3, {(1, 1) = .7, (1, 2) = .3, (1, 3) = .3, (2, 1) = .2, (2, 2) = .6, (2, 3) = .1, (3, 1) = .1, (3, 2) = .1, (3, 3) = .6})

macro(LA= LinearAlgebra):

The matrix I-A is singular, so it doesn't make sense to refer to 1/(I-A). Rather, we solve the system (I-A).X = 0.

LA:-LinearSolve(convert(LA:-IdentityMatrix(3)-A, rational), <0,0,0>);

Vector(3, {(1) = 2.50000000000000*_t[1], (2) = 1.50000000000000*_t[1], (3) = _t[1]})

eval(%, _t[1]= 100);

Vector(3, {(1) = 250., (2) = 150.000000000000, (3) = 100})

 


Download Leontief2.mw


restart:

params:= [
     z= 0, Omega= 2.2758, tau= 13.8,
     T2= 200, s= 1, r= 0.7071,

     s= 2.2758, Eta= 1.05457173*e-34,
     omega = 0.5, k = 1666666.667
];

[z = 0, Omega = 2.2758, tau = 13.8, T2 = 200, s = 1, r = .7071, s = 2.2758, Eta = 1.05457173*e-34, omega = .5, k = 1666666.667]

sys1 := {
     diff(u(t),t) =
          Omega*v(t)-u(t)/T2,

     diff(v(t),t) =
          -Omega*u(t) - v(t)/T2 -
          2*s*exp(-r^2/omega^2-t^2*1.177^2/tau^2)*
          cos(k*z-omega*t)*w(t),

     diff(w(t),t) =
          2*s*exp(-r^2/omega^2-t^2*1.177^2/tau^2)*
          cos(k*z-omega*t)*v(t)
};

 

ICs1 := {u(-20) = 0, v(-20) = 0, w(-20) = -1}:

{diff(u(t), t) = Omega*v(t)-u(t)/T2, diff(v(t), t) = -Omega*u(t)-v(t)/T2-2*s*exp(-r^2/omega^2-1.385329*t^2/tau^2)*cos(k*z-omega*t)*w(t), diff(w(t), t) = 2*s*exp(-r^2/omega^2-1.385329*t^2/tau^2)*cos(k*z-omega*t)*v(t)}

ans1:= dsolve(
     eval(sys1, params) union ICs1,
     numeric, output= listprocedure
):

 

plots:-odeplot(
     ans1, [[t,u(t)], [t,v(t)], [t,w(t)]],
     t= -20..20, legend= [u,v,w]
);

V:= eval(v(t), ans1):

x:= eval(2*V(t)*cos(k*z-omega*t), params);

2*V(t)*cos(.5*t)

plot(x, t= -20..20);

 


Download odesys.mw

GraphTheory:-Graph will accept an adjacency matrix but not an incidence matrix. Here is a procedure that will convert an incidence matrix into an adjacency matrix.

IncidenceToAdjacency:= proc(J::Matrix)
uses LA= LinearAlgebra;
local
     r:= LA:-RowDimension(J),
     A:= Matrix(r,r, shape= symmetric),
     c:= LA:-ColumnDimension(J),
     j, ones
;
     for j to c do
          ones:= ArrayTools:-SearchArray(J[..,j]);
          A[ones[1], ones[2]]:= 1
     end do;
     A
end proc:

(This procedure only handles the case of an unweighted, undirected graph.)

Using your example incidence Matrix:

M:=Matrix([[1,1,0,1,0,0,0],[0,1,1,0,1,0,0],[1,0,1,0,0,1,1],[0,0,0,1,1,1,0],[0,0,0,0,0,0,1]]):
A:= IncidenceToAdjacency(M);

G:= GraphTheory:-Graph(A):
GraphTheory:-DrawGraph(G);

Here is a complete simulation for the three-door Monty Hall problem (the classic). The code is very straightforward to read. See if you can extend this to the four-door problem.

MontyHall:= proc(N::posint, {switch::truefalse:= false})
# Simulation for the 3-door Monty Hall problem. N is the number
# of trials. Returns the number of trials the player wins.
local
     wins:= 0,
     PrizeDoor, #The door with the prize behind it (1..3).
     MontysDoor,  # The door that Monty reveals (1..3).
     Players1stDoor,  # The door the player picks first.
     Players2ndDoor,  # The door the player picks second.
     Doors:= {1,2,3},
     rand3:= rand(1..3), rand2:= rand(1..2)
;
     to N do
          PrizeDoor:= rand3();
          Players1stDoor:= rand3();
          if Players1stDoor = PrizeDoor then
               MontysDoor:= (Doors minus {PrizeDoor})[rand2()]
          else
               MontysDoor:= (Doors minus {PrizeDoor,Players1stDoor})[]
          end if;
          if switch then
               Players2ndDoor:= (Doors minus {MontysDoor,Players1stDoor})[]
          else
               Players2ndDoor:= Players1stDoor
          end if;
          if Players2ndDoor = PrizeDoor then
               wins:= wins+1
          end if
     end do;
     wins
end proc:     

MontyHall(100000);
                             33282
MontyHall(100000, switch);
                             66915

You asked:

How would i label the axis and change its font and size?

By using plot options labels and labelfont. (Example shown below.)

Also, how would i move the horizontal axis down, so that the whole diagram can be seen?

By using plot options axes= box and view. (Example shown below.)

Finally, there are some odd points between 0 and 1 on the vertical axis. How do you get rid of these?

Those are very interesting and seem to be an integral result of the computation for small values of r and some x, although it may also be an effect of round-off error. I will investigate that further. To eliminate them from the plot, slightly increase the lower bound for the view of the horizontal axis. I used 0.1 as the increase amount in the example show below.

Bonus: Here's how you can generate this plot in well under 1 second by using evalhf. It's at least a factor-of-50 improvement timewise over the old way. I set this up such that it will be easy to modify for other bifurcation diagrams.

restart:
Rloop:= proc(Pts::Matrix, r_min::realcons, r_max::realcons, N::posint, M::posint)
local n, r, x, r_range:= r_max - r_min;
     for n to N do
          r:= r_min+n/N*r_range;
          x:= Pts[n,2];
          to M do  x:= x*exp(r*(1-x))  end do;
          Pts[n,1]:= r;  Pts[n,2]:= x
     end do;
     NULL
end proc:     

(N,M):= 10^~(4,2):  (x_min,x_max):= (0,1):  (r_min,r_max):= (0,4):
bifpoint:= Matrix(N,2, datatype= float[8]):
bifpoint[.., 2]:=
     < RandomTools:-Generate(list(float(range= x_min..x_max, method= uniform), N))[] >:
evalhf(Rloop(bifpoint, r_min, r_max, N, M)):
pitchf:= plots:-pointplot(bifpoint, symbol=point):
plots:-display(
     pitchf,
     axes= box, view= [r_min+0.1..r_max, -1..5],
     labels= ['r','x'], labelfont= [TIMES,BOLDITALIC,16],
     axesfont= [HELVETICA,BOLD,12]
);

 

You can use the events option to dsolve to detect when the derivative is 0. If you post the equations, I'll help to set it up. Otherwise, you can read about it at ?dsolve,events . (Warning: This help page is quite esoteric.)

While it's true that you can't use solve with the output of dsolve(..., numeric), you can use fsolve.

I believe that your error "Empty script base" is related to the 2D input. I have never seen it before. When I convert your expression to 1D input with Convert tool on the Format menu, it looks very weird to me:

E~:=(k__x,k__y,k__z)~->~&+-Ɣsqrt((1)+(4~cos(3/(2)k__ya__cc)cos((((sqrt(3)))/(2))k__xa__cc))+((4)((cos^(2)(((sqrt(3))/(2))k__xa__cc))^()))6";

The tildes (~) and quotation marks make no sense to me.

Here is how I would type your expression in 1D input:

E:= (k__x, k__y, k__z)->
    &+- gamma*sqrt(1+4*cos(3/2*k__y*a__cc)*cos(sqrt(3)/2*k__x*a__cc) +
               4*cos(sqrt(3)/2*k__x*a__cc)^2);

To substitute the values for gamma and a__cc, do

E:= subs([gamma= 3*1.6e-19, a__cc= 0.142e-9], eval(E));

Names created with parse are always global. So you'll need to make PEA global instead of local for your code to work.

If you can describe what you're actually trying to accomplish, I can probably suggest a better way to do it.

Maple allows the upper bound of a plot range to be infinity:

plots([n^3, (n+1)^2], n= 3..infinity, legend= [n^3, (n+1)^2]);

or

plot(n^3 - (n+1)^2, n= 3..infinity);

Maple can prove some simple inequalities "on its own":

is(n^3 > (n+1)^2) assuming n >= 3;
                    
              true

 

First 328 329 330 331 332 333 334 Last Page 330 of 395