Robert Israel

6522 Reputation

21 Badges

18 years, 186 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are replies submitted by Robert Israel

>  with(Statistics):
   Coin1:= Sample(Bernoulli(4/5), 20000):
   Coin2:= Sample(Bernoulli(0), 20000):
   Coin3:= Sample(Bernoulli(1/2), 20000):
   A[0]:= 0: B[0]:= 0: C[0]:= 0:
   for t from 1 to 20000 do
     d1:= 2*round(Coin1[t])-1;
     d2:= 2*round(Coin2[t])-1;
     if A[t-1] mod 3 = 0 then A[t]:= A[t-1] + d2
     else A[t]:= A[t-1] + d1
     end if;
     if B[t-1] mod 3 = 1 then B[t]:= B[t-1] + d2
     else B[t]:= B[t-1] + d1
     end if;
     if C[t-1] mod 3 = Coin3[t] then
       C[t]:= C[t-1] + d2
     else
       C[t]:= C[t-1] + d1
     end if;
   end do:
  plot([[seq([i,A[i]],i=0..2000)],[seq([i,B[i]],i=0..2000)],[seq([i,C[i]],i=0..2000)]],
    symbol=point,style=point,colour=[red,blue,green]);

Note btw that in this case A[t] = B[t] + 2 for t >= 3.  Do you see why?

>  with(Statistics):
   Coin1:= Sample(Bernoulli(4/5), 20000):
   Coin2:= Sample(Bernoulli(0), 20000):
   Coin3:= Sample(Bernoulli(1/2), 20000):
   A[0]:= 0: B[0]:= 0: C[0]:= 0:
   for t from 1 to 20000 do
     d1:= 2*round(Coin1[t])-1;
     d2:= 2*round(Coin2[t])-1;
     if A[t-1] mod 3 = 0 then A[t]:= A[t-1] + d2
     else A[t]:= A[t-1] + d1
     end if;
     if B[t-1] mod 3 = 1 then B[t]:= B[t-1] + d2
     else B[t]:= B[t-1] + d1
     end if;
     if C[t-1] mod 3 = Coin3[t] then
       C[t]:= C[t-1] + d2
     else
       C[t]:= C[t-1] + d1
     end if;
   end do:
  plot([[seq([i,A[i]],i=0..2000)],[seq([i,B[i]],i=0..2000)],[seq([i,C[i]],i=0..2000)]],
    symbol=point,style=point,colour=[red,blue,green]);

Note btw that in this case A[t] = B[t] + 2 for t >= 3.  Do you see why?

To have an animation of a surface, you need three independent variables (two space and one time).  You have only two independent variables.  You could consider them both as space variables, producing a surface, or you can consider one as time and the other as space, producing an animation of a curve (as I did in my example).  

To have an animation of a surface, you need three independent variables (two space and one time).  You have only two independent variables.  You could consider them both as space variables, producing a surface, or you can consider one as time and the other as space, producing an animation of a curve (as I did in my example).  

The only difference is that you would use a*exp(b*x) instead of a*x^b.  Not a*e^(b*x), because Maple considers e to be an ordinary variable.


 

The only difference is that you would use a*exp(b*x) instead of a*x^b.  Not a*e^(b*x), because Maple considers e to be an ordinary variable.


 

If the lists have the same length, you can simply subtract:

> A-B;

 If not, you can use

> zip(`-`, A, B, 0);

 

which will pad the shorter one with 0's. 

If the lists have the same length, you can simply subtract:

> A-B;

 If not, you can use

> zip(`-`, A, B, 0);

 

which will pad the shorter one with 0's. 

indets finds the set of all sub-expressions of an expression that have a given type.  In this case the type is specfunc(anything,CURVES), which means any function call of the form CURVES(...).  Any plot command that produces curves will contain such objects.  Then op([1,1], ...) returns the first operand of the first operand of the set, i.e. the first operand of the CURVES function call, which is generally a list (but sometimes an Array) of points.

 

indets finds the set of all sub-expressions of an expression that have a given type.  In this case the type is specfunc(anything,CURVES), which means any function call of the form CURVES(...).  Any plot command that produces curves will contain such objects.  Then op([1,1], ...) returns the first operand of the first operand of the set, i.e. the first operand of the CURVES function call, which is generally a list (but sometimes an Array) of points.

 

You might be able to increase the stacklimit, if that is below the system's hard limit.  See ?kernelopts.
Alternatively, you might try using global variables rather than local variables of the procedures for storing your large data.

 

You might be able to increase the stacklimit, if that is below the system's hard limit.  See ?kernelopts.
Alternatively, you might try using global variables rather than local variables of the procedures for storing your large data.

 

The following procedure will generate an arrow of length L starting at the point on the solution Sol to the system Sys for t=T.  It is assumed that Sys is of the form {D(x)(t) = ..., D(y)(t) = ...} and Sol is the result of a call to dsolve({op(Sys), ICs}, numeric, output=procedurelist).

> arrowonsol:= proc(Sys, Sol, T, L)
         uses plots;
         local S, P, V;
         S:= Sol(T);
         P:= subs(S, [x(t), y(t)]);
         V:= evalf(subs(S, subs(Sys, [D(x)(t), D(y)(t)])));
         V:= L*V/sqrt(V[1]^2 + V[2]^2);
         plots[arrow](P,V,args[5..-1]);
    end proc;
  

For example:

>  with(plots):
   Sys:= [(D(x))(t) = -.1+x(t)^2-x(t)*y(t), (D(y))(t) = y(t)^2-x(t)^2-1];   
   Sol1:= dsolve([op(Sys), x(0)=0,y(0)=-1],numeric,output=procedurelist);
   Sol2:= dsolve([op(Sys), x(0) = 0.2, y(0) = -1.02], numeric, output=procedurelist);
   Curves:= odeplot(Sol1, [x(t),y(t)], t=-5..5),odeplot(Sol2, [x(t),y(t)], t=-5..5):
   display([Curves,
        seq(arrowonsol(Sys,Sol1,i,0.3,width=0.1),i=[-4,3]),
        seq(arrowonsol(Sys,Sol2,i,0.3,width=0.1),i=[-2.5,1.5])], 
     scaling=constrained);

The following procedure will generate an arrow of length L starting at the point on the solution Sol to the system Sys for t=T.  It is assumed that Sys is of the form {D(x)(t) = ..., D(y)(t) = ...} and Sol is the result of a call to dsolve({op(Sys), ICs}, numeric, output=procedurelist).

> arrowonsol:= proc(Sys, Sol, T, L)
         uses plots;
         local S, P, V;
         S:= Sol(T);
         P:= subs(S, [x(t), y(t)]);
         V:= evalf(subs(S, subs(Sys, [D(x)(t), D(y)(t)])));
         V:= L*V/sqrt(V[1]^2 + V[2]^2);
         plots[arrow](P,V,args[5..-1]);
    end proc;
  

For example:

>  with(plots):
   Sys:= [(D(x))(t) = -.1+x(t)^2-x(t)*y(t), (D(y))(t) = y(t)^2-x(t)^2-1];   
   Sol1:= dsolve([op(Sys), x(0)=0,y(0)=-1],numeric,output=procedurelist);
   Sol2:= dsolve([op(Sys), x(0) = 0.2, y(0) = -1.02], numeric, output=procedurelist);
   Curves:= odeplot(Sol1, [x(t),y(t)], t=-5..5),odeplot(Sol2, [x(t),y(t)], t=-5..5):
   display([Curves,
        seq(arrowonsol(Sys,Sol1,i,0.3,width=0.1),i=[-4,3]),
        seq(arrowonsol(Sys,Sol2,i,0.3,width=0.1),i=[-2.5,1.5])], 
     scaling=constrained);

Here's another version.  This time the turn number doesn't matter: which game you play is decided randomly. 

You have three coins, Coin 1 which wins with probability 4/5, Coin 2 which always loses, and Coin 3 (used only in Game C to decide whether to play a turn of A or B) which is fair: probability 1/2 for heads and 1/2 for tails.

Game A: if fortune divisible by 3, flip Coin 2, otherwise flip Coin 1.

Game B: if fortune congruent to 1 mod 3, flip Coin 2, otherwise flip Coin 1.

Game C: Flip Coin 3.  If heads play a turn of Game A, if tails play a turn of Game B.  .

In Game A, the three states (fortune == 0 mod 3, 1 mod 3 and 2 mod 3) have long-run probabilities  21/51, 5/51 and 25/51 respectively. 
Your expected gain per turn is (-1)*21/51 + 3/5*5/51 + 3/5*25/51 = -1/17.  Similarly in Game B.

In Game C, you always flip Coin 1 if your fortune == 2 mod 3, otherwise effectively you're flipping a coin with winning probability 2/5.  The long-run
probabilities of the three states are then 13/55, 19/55 and 23/55, and the expected net gain per turn is (-1/10)*13/55 + (-1/10)*19/55 + 3/5*23/55 = 53/275.

 

First 48 49 50 51 52 53 54 Last Page 50 of 187