Carl Love

Carl Love

27423 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Your reqn has two solutions for f(n+1), and each of them specifies a different recurrence:

psol:= solve({f(n+1)=f(n)+sqrt(2*f(n+1)-f(n))}, {f(n+1)});

I don't think that there's any hope of a symbolic solution for either of these (and the rsolves return unevaluated), but you may be able to make a numeric asymptotic analysis. In the following command, note that I've specified the initial condition as 15., with a decimal point, to force floating-point evaluations. If you don't do this, the expressions become extremely large with deeply nested square roots very fast.

F:= rsolve({psol[1], f(1)=15.}, f(n), makeproc);

Now you can explore, for example,

seq(F(k), k= 1..99);

If I used psol[2] instead, the solution rapidly decays to 0.

Using package VectorCalculus, I'd do it with the commands RootedVector and PlotVector (definitely not PlotPositionVector). For example, we're given the following initial tail point, base, and vectors u and v:

base:= [1,1,1]:  u:= <1,2,3>:  v:= <4,5,6>: 

We want to represent the vector sum u + v head-to-tail in that order:

VC:= VectorCalculus:
VC:-PlotVector([VC:-RootedVector(root= base, u), VC:-RootedVector(root= base+~u, v)]);

There are numerous options to PlotVector to control the size, color, shape, etc., of the vectors; see its help page.

The same process works just as well for a 2D plot if all vectors and points have 2 components.

From a Reply, I see that you really want to generate discrete uniform random numbers from 1..120 using a 6-sided die (note that dice is strictly plural, like mice and lice). The process can be simplified because 120 can be expressed as a product of positive integers all <= 6; specifically, 120 = 4*5*6. (Number theorists would say that 120 is {2, 3, 5}-smooth, the {2, 3, 5} being the primes <= 6.) Since those are three factors, it'll often be possible with three rolls of the 6-sided die. However, sometimes a roll will need to be discarded and rerolled. The expected number of rolls is exactly 3.7 (Maple-based proof below).

I use a multi-radix representation. This is similar to an alternate-base representation, but with a different base for each digit position. For any n in {$120}, there exists a unique triple [d4, d5, d6] such that

n = 1 + 5*6*(d4 - 1) + 6*(d5 - 1) + (d6 - 1),   (formula 1)

where d4 in {$4}, d5 in {$5}, and d6 in {$6}[*1] (proof outline in footnote). I developed this formula to make the arithmetic simple enough for your son. So, how do we generate d4 and d5 uniformly with a 6-sided die? For d4, roll the die, perhaps repeatedly, until it gives a number <= 4; for d5, do the analogous thing. The expected number of rolls for d_k (i.e., d4d5, etc.) is 

sum(n*k/s*(1 - k/s)^(n-1), n= 1..infinity) assuming 0 < k, k <= s,

where s is the number of sides on the actual die, and n runs over all possible numbers of rolls. The sum simplifies to s/k, as you can check with Maple. (You can also look up the negative binomial distribution in Wikipedia.) So, in this particular case the expected number of rolls to get all three "digits" is 6/4 + 6/5 + 6/6 = 37/10 = 3.7

If n is not {2, 3, 5}-smooth, a slightly more-complicated process can be used, and the expected number of rolls will be slightly more.

Here's a Maple simulation for n = 120, which I coded to be easily modifiable to other numbers of sides of the die and appropriately smooth values of n

restart
:
Roll_120:= module()
local
    sides::And(posint, Not(1)):= 6, 
    ds:= rand(1..sides),
    roll:= (k::And(posint, Not(1)))-> local r; 
        do Roll_Ct++; if (r:= ds()) <= k then return r fi od,  
    radices:= [4,5,6],  #4*5*6 = 120
    Coeffs:= (P-> P[-1]/~P)([seq[scan= `*`]](radices))  #[30, 6, 1]
;
export
    ModuleApply:= ()-> (R-> R = 1 + inner(R-~1, Coeffs))(roll~(radices)),
    ExpectedRolls:= add(sides/~radices),  #expected # of rolls per turn
    Roll_Ct:= 0  #Count the rolls only to verify expectation.
;
end module
:     
#Choose the hockey team!:
'Roll_120'() $ 5;
[2, 5, 6] = 60, [2, 3, 4] = 46, [4, 5, 3] = 117, [1, 5, 2] = 26, [3, 2, 2] = 68
#If a card number were repeated, we'd simply reroll for that card.

Roll_120:-ExpectedRolls;
                               37 / 10

#Performance test for 100,000 triples:
Roll_120:-Roll_Ct:= 0:
Rolls:= CodeTools:-Usage([to 10^5 do Roll_120() od]):
memory used=74.07MiB, alloc change=-3.01MiB, 
cpu time=766.00ms, real time=727.00ms, gc time=62.50ms

Roll_120:-Roll_Ct;  #Check expected roll count.
                             369526

#Check uniformity:
RRolls:= rhs~(Rolls):
[{Statistics:-Tally(RRolls)[]}[]];
[1 = 785, 2 = 822, 3 = 821, 4 = 827, 5 = 777, 6 = 806, 7 = 874, 8 = 857, 9 = 823, 10 = 837,
 11 = 827, 12 = 811, 13 = 856, 14 = 820, 15 = 791, 16 = 874, 17 = 861, 18 = 848, 19 = 796, 20 = 855, 
 21 = 850, 22 = 826, 23 = 856, 24 = 845, 25 = 861, 26 = 829, 27 = 806, 28 = 890, 29 = 787, 30 = 846,
 31 = 873, 32 = 841, 33 = 836, 34 = 853, 35 = 855, 36 = 854, 37 = 779, 38 = 842, 39 = 854, 40 = 827, 
 41 = 903, 42 = 814, 43 = 816, 44 = 797, 45 = 829, 46 = 835, 47 = 854, 48 = 819, 49 = 861, 50 = 868, 
 51 = 824, 52 = 802, 53 = 781, 54 = 880, 55 = 802, 56 = 858, 57 = 826, 58 = 867, 59 = 807, 60 = 839,
 61 = 851, 62 = 798, 63 = 843, 64 = 828, 65 = 832, 66 = 794, 67 = 850, 68 = 821, 69 = 821, 70 = 812,
 71 = 809, 72 = 831, 73 = 819, 74 = 805, 75 = 811, 76 = 842, 77 = 805, 78 = 850, 79 = 816, 80 = 799,
 81 = 809, 82 = 851, 83 = 805, 84 = 863, 85 = 818, 86 = 841, 87 = 863, 88 = 893, 89 = 865, 90 = 784,
 91 = 840, 92 = 808, 93 = 833, 94 = 851, 95 = 889, 96 = 797, 97 = 843, 98 = 804, 99 = 798, 100 = 856,
 101 = 862, 102 = 826, 103 = 840, 104 = 883, 105 = 825, 106 = 882, 107 = 879, 108 = 799, 109 = 804, 110 = 851,
 111 = 834, 112 = 852, 113 = 819, 114 = 807, 115 = 861, 116 = 846, 117 = 829, 118 = 837, 119 = 825, 120 = 831]

nops(%), add(rhs~(%));
                          120, 100000

Statistics:-Histogram(RRolls, discrete, frequencyscale= absolute, size= 70*[8,3]);

The process that you show in the Question has two major mathematical errors. (This is totally independent of the Maple problem of automatic simplification / premature evaluation that's already been addressed.):

  1. A sum of independent nonconstant uniform random variables can never be uniform.
  2. If r is a discrete uniform random variable on 1..n, then r mod d can only be uniform if d divides (or, quite trivially, if abs(d) > n).


Footnote:

[*1] Outline of a proof of the existence and uniqueness of the multi-radix representation

n = 1 + 5*6*(d4 - 1) + 6*(d5 - 1) + (d6 - 1).

Let d6:= irem(n-1, 6, 'q') + 1d5:= irem(q, 5, 'q') + 1d4:= irem(q, 4) + 1. Each of these steps is defined and produces a unique value.

The Submit Software Change Request feature of MaplePrimes works. You may have simply experienced some temporary server problem or Internet problem.

My mental solution without Maple:

rel(3) is obviously true; rel(2) can be verified by elementary arithmetic. For n > 3, divide both sides by ((n-3)*2)^(1/(n-1)). The equation becomes (n-3)^1 * 2^1 = (n-1)*2^1 - 4, which reduces to 2*n-6 = 2*n-6.

Here's a workaround that's as close to the "spirit" of the original command as I think is possible. It uses subs to perform simultaneous substitutions by making the first argument to subs a list.

subs(op(expr)=~ [a,b,c,d], expr)

(As always, my posting of a "workaround" is unrelated to whether I think there's a bug. In this case, I do think that there's a bug.)

Here's my workaround to find the index of the maximum entry in a table:

T:= table(['rand()'$2^18] =~ ['rand()'$(2^18-1), 2^64]):
P:= [indices](T, pairs):
j:= lhs(P[max[index](rhs~(P))]);

                        j := 65671522612
T[j];
                      18446744073709551616
2^64;
                      18446744073709551616

This works exactly the same if the table has a multi-index.

You need

map(myf@op, mylist)

The @op converts the inner lists to argument sequences before passing them to myf.

The transformation that you propose has already been implemented in the reverse direction, the right side of your equation becoming the left side. In other words, the following happens automatically:

diff(inttrans:-fourier(U(x,t), x, s), t);  lprint(%);
fourier(diff(U(x,t),t),x,s)

I'm not writing that to say that what you want can't be done; I'm just saying that we need to take it into account when implementing your proposal. If it's not taken into account, your transformed result (i.e., diff@fourier) will automatically revert to its original (i.e., fourier@diff). It is for this reason that my implementation below uses the inert Diff rather than the usual diff.

map2(
    inttrans:-addtable,
    fourier,
    [diff, Diff](U(x,t),t),
(* => *) Diff(inttrans:-fourier(U(_x,t), _x, s), t),
    x, s, {U(x,t), t}
): 
#There's no response if this command is accepted.

Test it:

inttrans:-fourier(diff(U(x,t),t), x, s);  lprint(%);
Diff(fourier(U(_x,t),_x,s),t)

I changed x to _x in the result to show that it's not really the same x. This is akin to @sand15's xi. Since they're bound variables, I don't think it makes any practical difference.

Unfortunately, this only implements your transform for the specific function name U. I may be able to figure out something more general later..

There are eight groups of three elimination variables that can be used, all of which yield the desired polynomial. The process shown in the following worksheet finds all of them quickly and automatically.

restart:

Peng-Robinson_eos

In the following computation, the variables in each numerator are never separated, so we can reduce the number of variables by considering each numerator as a single variable.

eq4:= P = RT/(v-b) - aT/(v*(v+b)+b*(v-b));

P = RT/(v-b)-aT/(v*(v+b)+b*(v-b))

elims:= {(A,B,Z) =~ P/RT *~ (aT/RT, b, v)};

{A = P*aT/RT^2, B = P*b/RT, Z = P*v/RT}

So there are potentially five variables that we can eliminate:

Vs:= indets(rhs~(elims));

{P, RT, aT, b, v}

But since there are only three elimination equations, we need to pick three variables. The possible groups of three are

V3:= [combinat:-choose(Vs, nops(elims))[]];

[{P, RT, aT}, {P, RT, b}, {P, RT, v}, {P, aT, b}, {P, aT, v}, {P, b, v}, {RT, aT, b}, {RT, aT, v}, {RT, b, v}, {aT, b, v}]

Use solve on each group of three variables:

sols:= V3 =~ map2({solve}, elims, V3);

[{P, RT, aT} = {}, {P, RT, b} = {{P = Z^2*aT/(A*v^2), RT = Z*aT/(A*v), b = B*v/Z}}, {P, RT, v} = {{P = B^2*aT/(A*b^2), RT = B*aT/(A*b), v = Z*b/B}}, {P, aT, b} = {{P = Z*RT/v, aT = A*RT*v/Z, b = B*v/Z}}, {P, aT, v} = {{P = B*RT/b, aT = A*RT*b/B, v = Z*b/B}}, {P, b, v} = {{P = A*RT^2/aT, b = B*aT/(A*RT), v = Z*aT/(A*RT)}}, {RT, aT, b} = {{RT = P*v/Z, aT = A*P*v^2/Z^2, b = B*v/Z}}, {RT, aT, v} = {{RT = P*b/B, aT = A*P*b^2/B^2, v = Z*b/B}}, {RT, b, v} = {{RT = RootOf(A*_Z^2-P*aT), b = B*RootOf(A*_Z^2-P*aT)/P, v = Z*RootOf(A*_Z^2-P*aT)/P}}, {aT, b, v} = {{aT = A*RT^2/P, b = B*RT/P, v = Z*RT/P}}]

Obviously, we can't use the one whose solution set is empty. The one with RootOf could also be problematic, so we'll exclude it also, at least for now.

ok_sols:= ((lhs=op@rhs)~@remove)(
    (r-> has(r, RootOf) or nops(r) <> 1)@rhs,
    sols
);

[{P, RT, b} = {P = Z^2*aT/(A*v^2), RT = Z*aT/(A*v), b = B*v/Z}, {P, RT, v} = {P = B^2*aT/(A*b^2), RT = B*aT/(A*b), v = Z*b/B}, {P, aT, b} = {P = Z*RT/v, aT = A*RT*v/Z, b = B*v/Z}, {P, aT, v} = {P = B*RT/b, aT = A*RT*b/B, v = Z*b/B}, {P, b, v} = {P = A*RT^2/aT, b = B*aT/(A*RT), v = Z*aT/(A*RT)}, {RT, aT, b} = {RT = P*v/Z, aT = A*P*v^2/Z^2, b = B*v/Z}, {RT, aT, v} = {RT = P*b/B, aT = A*P*b^2/B^2, v = Z*b/B}, {aT, b, v} = {aT = A*RT^2/P, b = B*RT/P, v = Z*RT/P}]

All eight of the remaining combinations yield the desired residual polynomial:

for V in lhs~(ok_sols) do
    sol[V[]]:= eliminate({eq4,elims[]}, V)[2][]
od;

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

(collect(sol[aT,b,v], Z) = 0) &where subs(RT= R*T, elims);

`&where`(Z^3+(B-1)*Z^2+(-3*B^2+A-2*B)*Z+B^3-A*B+B^2 = 0, {A = P*aT/(R^2*T^2), B = P*b/(R*T), Z = P*v/(R*T)})

Verify that the substitutions return the original eq4:

solve(eval(op(1,%), op(2,%)), P);

0, 0, -(R*T*b^2-2*R*T*b*v-R*T*v^2-aT*b+aT*v)/(b^3-3*b^2*v+b*v^2+v^3)

convert(%[3], parfrac, v);

-aT/(-b^2+2*b*v+v^2)+R*T/(v-b)

 

Download Peng_Robinson.mw

AFAIK, pdsolve(..., numericcan only handle PDE systems with 2 differentiation variables. You have X, R, and t.

Here are two ways:

subsindets(expr, specindex(a), x-> b[op(x)])

Or, if you know that all the a's that you want have exactly 1 index:

subsindets(expr, a[anything], x-> b[op(x)])

How do you expect it to react when you give it a y-range that's out-of-bounds (into complex numbers) for the function? Narrow the y-range to -1.53..1.53.

The command FunctionAdvisor will give you a huge amount of information about a special function. Try

FunctionAdvisor(Ei[1](x)).

The notations Ei[a](x) and Ei(a, x) seem to be used interchangeably, and

 Ei(x) = Ei(0, x) = Ei[0](x) = exp(-x)/x.

Here's another way, just to show the possibilities. I often prefer index-free methods, which let me treat matrices and vectors as whole mathematical objects in their own right rather than "poking inside" them.

I construct a matrix of products of powers of X and by taking the outer product of a vector of powers of and a vector of powers of Y. Then I take the elementwise product of that with m. Then use add, which can be used without indices.

(X,Y):= (2,3):  (r,c):= op(1, m) -~ 1:
add(<X^~($0..r)>.<Y^~($0..c)>^%T *~ m);

                           
-11.518

I will readily admit that the above, done on a large scale, will increase your garbage collection time. 

First 12 13 14 15 16 17 18 Last Page 14 of 391