Advanced Search

Carl Love

Carl Love

28110 Reputation

25 Badges

13 years, 119 days
Himself
Wayland, Massachusetts, United States
Contact Carl Love
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

My mental process...

Carl Love 28110
June 03 2013

@Markiyan Hirnyk 

Here I present some insight into my mental process in analyzing such a sequence and finding the patterns. It is mostly mental arithmetic and mental pattern recognition, but Maple helps a little. Here's a worksheet:

 

> 

A:=
[5, 7, 11, 13, 15, 19, 21, 29, 31, 33, 37, 39, 45, 55, 57, 63, 83, 85, 87, 91, 93, 99, 109,
 111, 117, 135, 163, 165, 171, 189, 245, 247, 249, 253, 255, 261, 271, 273, 279, 297, 325,
 327, 333, 351, 405, 487, 489, 495, 513, 567, 731, 733, 735, 739, 741, 747, 757, 759, 765,
 783, 811, 813, 819, 837, 891, 973, 975, 981, 999, 1053, 1215, 1459, 1461, 1467, 1485, 1539,
 1701, 2189, 2191, 2193, 2197, 2199, 2205, 2215, 2217, 2223, 2241, 2269, 2271, 2277, 2295,
 2349, 2431, 2433, 2439, 2457, 2511, 2673, 2917, 2919, 2925, 2943, 2997, 3159, 3645, 4375,
 4377, 4383, 4401, 4455, 4617, 5103, 6563, 6565, 6567, 6571, 6573, 6579, 6589, 6591, 6597,
 6615, 6643, 6645, 6651, 6669, 6723, 6805, 6807, 6813, 6831, 6885, 7047, 7291, 7293, 7299,
 7317, 7371, 7533, 8019, 8749, 8751, 8757, 8775, 8829, 8991, 9477, 10935, 13123, 13125,
 13131, 13149, 13203, 13365, 13851, 15309, 19685, 19687, 19689, 19693, 19695, 19701, 19711,
 19713, 19719, 19737, 19765, 19767, 19773, 19791, 19845, 19927, 19929, 19935, 19953, 20007,
 20169, 20413, 20415, 20421, 20439, 20493, 20655, 21141, 21871, 21873, 21879, 21897, 21951,
 22113, 22599, 24057, 26245, 26247, 26253, 26271, 26325, 26487, 26973, 28431, 32805, 39367,
 39369, 39375, 39393, 39447, 39609, 40095, 41553, 45927]:

I notice that the sequence is strictly increasing, so I will look at the first differences. Further, every term is odd, so I might as well divide the differences by two.

> 

A1:= [seq]((A[k]-A[k-1])/2, k= 2..nops(A));

[1, 2, 1, 1, 2, 1, 4, 1, 1, 2, 1, 3, 5, 1, 3, 10, 1, 1, 2, 1, 3, 5, 1, 3, 9, 14, 1, 3, 9, 28, 1, 1, 2, 1, 3, 5, 1, 3, 9, 14, 1, 3, 9, 27, 41, 1, 3, 9, 27, 82, 1, 1, 2, 1, 3, 5, 1, 3, 9, 14, 1, 3, 9, 27, 41, 1, 3, 9, 27, 81, 122, 1, 3, 9, 27, 81, 244, 1, 1, 2, 1, 3, 5, 1, 3, 9, 14, 1, 3, 9, 27, 41, 1, 3, 9, 27, 81, 122, 1, 3, 9, 27, 81, 243, 365, 1, 3, 9, 27, 81, 243, 730, 1, 1, 2, 1, 3, 5, 1, 3, 9, 14, 1, 3, 9, 27, 41, 1, 3, 9, 27, 81, 122, 1, 3, 9, 27, 81, 243, 365, 1, 3, 9, 27, 81, 243, 729, 1094, 1, 3, 9, 27, 81, 243, 729, 2188, 1, 1, 2, 1, 3, 5, 1, 3, 9, 14, 1, 3, 9, 27, 41, 1, 3, 9, 27, 81, 122, 1, 3, 9, 27, 81, 243, 365, 1, 3, 9, 27, 81, 243, 729, 1094, 1, 3, 9, 27, 81, 243, 729, 2187, 3281, 1, 3, 9, 27, 81, 243, 729, 2187]

Many of the terms in A1 look familiar, powers of 3. Apply ifactor.

> 

A2:= map(ifactor, A1);

[1, ``(2), 1, 1, ``(2), 1, ``(2)^2, 1, 1, ``(2), 1, ``(3), ``(5), 1, ``(3), ``(2)*``(5), 1, 1, ``(2), 1, ``(3), ``(5), 1, ``(3), ``(3)^2, ``(2)*``(7), 1, ``(3), ``(3)^2, ``(2)^2*``(7), 1, 1, ``(2), 1, ``(3), ``(5), 1, ``(3), ``(3)^2, ``(2)*``(7), 1, ``(3), ``(3)^2, ``(3)^3, ``(41), 1, ``(3), ``(3)^2, ``(3)^3, ``(2)*``(41), 1, 1, ``(2), 1, ``(3), ``(5), 1, ``(3), ``(3)^2, ``(2)*``(7), 1, ``(3), ``(3)^2, ``(3)^3, ``(41), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(2)*``(61), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(2)^2*``(61), 1, 1, ``(2), 1, ``(3), ``(5), 1, ``(3), ``(3)^2, ``(2)*``(7), 1, ``(3), ``(3)^2, ``(3)^3, ``(41), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(2)*``(61), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(5)*``(73), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(2)*``(5)*``(73), 1, 1, ``(2), 1, ``(3), ``(5), 1, ``(3), ``(3)^2, ``(2)*``(7), 1, ``(3), ``(3)^2, ``(3)^3, ``(41), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(2)*``(61), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(5)*``(73), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(3)^6, ``(2)*``(547), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(3)^6, ``(2)^2*``(547), 1, 1, ``(2), 1, ``(3), ``(5), 1, ``(3), ``(3)^2, ``(2)*``(7), 1, ``(3), ``(3)^2, ``(3)^3, ``(41), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(2)*``(61), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(5)*``(73), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(3)^6, ``(2)*``(547), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(3)^6, ``(3)^7, ``(17)*``(193), 1, ``(3), ``(3)^2, ``(3)^3, ``(3)^4, ``(3)^5, ``(3)^6, ``(3)^7]

Now you can see why my first guess was to add 2*3^8 to the last term (recall that I divided by 2 above). But I missed a small part of the cyclic pattern: The last group of powers of 3 repeats, followed by 2 times the number before the preceding 1. That's why my second guess was the last term plus 4*17*193.

 

After studying sequence A2 for a while, the cycles and epicycles become clear to me. I did not use any plots in my analysis, but the plot below shows the cycles nicely.

> 

P1:= plot(
     [seq]([k, log[3](A1[k])], k= 1..nops(A1)),
     style= point, symbol= diamond, symbolsize= 8,
     gridlines
):
> 

P2:= plot(
     [seq]([k, log[3](A1[k])], k= 1..nops(A1)),
     thickness= 0, color= blue
):
> 

plots:-display([P1,P2]);

At this point, I did not understand the numbers like 17*193 that occasionally appeared in the pattern, but I knew exactly where they occurred. So, at this point I knew that I could generate another 54 = binomial(11,2)-1 terms of the sequence before I encountered another such number.

 

Now I analyze those special numbers, which also form an increasing sequence. I pick off those numbers starting from the end of the sequence A2, because the larger such numbers are easier to see in the pattern. It is harder to see the 1 and 2 as part of this subsequence because 1s and 2s appear in the pattern for other reasons also. But by carefully counting the position numbers within the patterns, I find them. I entered these by hand on the following line as I found them:

> 

17*193, 2*547, 5*73, 2*61, 41, 2*7, 5, 2, 1;

3281, 1094, 365, 122, 41, 14, 5, 2, 1

I immediately see that each term in this sequence is 3 times the previous minus 1. Use rsolve.

> 

rsolve({B(n)= 3*B(n-1)-1, B(0)=1}, B(n));

(1/2)*3^n+1/2

At this point, I completely understand the pattern of cycles and epicycles (and epi-epicycles), although I don't have a closed form or even a recursive form. But I have enough information to write a procedure that generates a whole number of cycles.

> 

A:= proc(n::{nonnegint,identical(-1)})
local B, R, i, j, k;
     B:= 1, 2,

          '[1,

               '[3^i $ i= 0..j, (3^(j+1)+1)/2][]'

               $ j= 0..k,
          3^i $ i= 0..k, 3^(k+1)+1

          ][]' $ k= 0..n

     ;
     R[0]:= 5;
     for k to nops([B]) do  R[k]:= R[k-1] + 2*B[k]  end do;
     convert(R, list)
end proc:

Bonus: I want a closed form for the number of terms that A(n) generates. This could be determined directly from an analysis of the code, but I think pattern finding is easier.

> 

seq(nops(A(k)), k= 0..8);

8, 17, 31, 51, 78, 113, 157, 211, 276

It's an increasing sequence: Look at first differences.

> 

seq(nops(A(k+1))-nops(A(k)), k= 0..8);

9, 14, 20, 27, 35, 44, 54, 65, 77

I immediately recognize that as the following sequence:

> 

seq(binomial(n,2)-1, n= 5..13);

9, 14, 20, 27, 35, 44, 54, 65, 77

> 

rsolve({C(n)= C(n-1)+binomial(n+4,2)-1, C(0)=8}, C(n));

5+(n+1)*((1/2)*n+1)*((1/3)*n+1)+2*n+2*(n+1)*((1/2)*n+1)

> 

expand(%);

8+(1/6)*n^3+2*n^2+(41/6)*n

> 

numer(%);

n^3+12*n^2+41*n+48

> 

 

 

Download SeqAnalysis.mw

Vote up?...

Carl Love 28110
June 03 2013
0 0

@Carl Love Why doesn't that get a vote up?

Vote up?...

Carl Love 28110
June 03 2013

@Carl Love Why doesn't that get a vote up?

A procedure to generate the sequence...

Carl Love 28110
June 03 2013
0 0

@Markiyan Hirnyk 

The following procedure, when called with nonnegint argument n, will generate the first

 (n^3 + 12*n^2 + 41*n + 48)/6

terms of the sequence.

A:= proc(n::nonnegint)
local B, R, i, j, k;
     B:= [1, 2, seq([1, seq([seq(3^i, i= 0..j), (3^(j+1)+1)/2][], j= 0..k),
          seq(3^i, i= 0..k), 3^(k+1)+1][], k= 0..n)];
     R[0]:= 5;
     for k to nops(B) do  R[k]:= R[k-1] + 2*B[k]  end do;
     convert(R, list)
end proc:

A(8);
[5, 7, 11, 13, 15, 19, 21, 29, 31, 33, 37, 39, 45, 55, 57, 63, 
  83, 85, 87, 91, 93, 99, 109, 111, 117, 135, 163, 165, 171, 189, 
 245, 247, 249, 253, 255, 261, 271, 273, 279, 297, 325, 327, 
333, 351, 405, 487, 489, 495, 513, 567, 731, 733, 735, 739, 
741, 747, 757, 759, 765, 783, 811, 813, 819, 837, 891, 973, 
975, 981, 999, 1053, 1215, 1459, 1461, 1467, 1485, 1539, 1701, 
 2189, 2191, 2193, 2197, 2199, 2205, 2215, 2217, 2223, 2241, 
2269, 2271, 2277, 2295, 2349, 2431, 2433, 2439, 2457, 2511, 
2673, 2917, 2919, 2925, 2943, 2997, 3159, 3645, 4375, 4377, 
4383, 4401, 4455, 4617, 5103, 6563, 6565, 6567, 6571, 6573, 
6579, 6589, 6591, 6597, 6615, 6643, 6645, 6651, 6669, 6723, 
6805, 6807, 6813, 6831, 6885, 7047, 7291, 7293, 7299, 7317, 
7371, 7533, 8019, 8749, 8751, 8757, 8775, 8829, 8991, 9477, 
10935, 13123, 13125, 13131, 13149, 13203, 13365, 13851, 15309, 
19685, 19687, 19689, 19693, 19695, 19701, 19711, 19713, 19719, 
19737, 19765, 19767, 19773, 19791, 19845, 19927, 19929, 19935, 
19953, 20007, 20169, 20413, 20415, 20421, 20439, 20493, 20655, 
21141, 21871, 21873, 21879, 21897, 21951, 22113, 22599, 24057, 
 26245, 26247, 26253, 26271, 26325, 26487, 26973, 28431, 32805, 
 39367, 39369, 39375, 39393, 39447, 39609, 40095, 41553, 45927, 
59051, 59053, 59055, 59059, 59061, 59067, 59077, 59079, 59085, 
59103, 59131, 59133, 59139, 59157, 59211, 59293, 59295, 59301, 
59319, 59373, 59535, 59779, 59781, 59787, 59805, 59859, 60021, 
60507, 61237, 61239, 61245, 61263, 61317, 61479, 61965, 63423, 
65611, 65613, 65619, 65637, 65691, 65853, 66339, 67797, 72171, 
78733, 78735, 78741, 78759, 78813, 78975, 79461, 80919, 85293, 
98415, 118099, 118101, 118107, 118125, 118179, 118341, 118827, 
120285, 124659, 137781, 177149]

A procedure to generate the sequence...

Carl Love 28110
June 03 2013

@Markiyan Hirnyk 

The following procedure, when called with nonnegint argument n, will generate the first

 (n^3 + 12*n^2 + 41*n + 48)/6

terms of the sequence.

A:= proc(n::nonnegint)
local B, R, i, j, k;
     B:= [1, 2, seq([1, seq([seq(3^i, i= 0..j), (3^(j+1)+1)/2][], j= 0..k),
          seq(3^i, i= 0..k), 3^(k+1)+1][], k= 0..n)];
     R[0]:= 5;
     for k to nops(B) do  R[k]:= R[k-1] + 2*B[k]  end do;
     convert(R, list)
end proc:

A(8);
[5, 7, 11, 13, 15, 19, 21, 29, 31, 33, 37, 39, 45, 55, 57, 63, 
  83, 85, 87, 91, 93, 99, 109, 111, 117, 135, 163, 165, 171, 189, 
 245, 247, 249, 253, 255, 261, 271, 273, 279, 297, 325, 327, 
333, 351, 405, 487, 489, 495, 513, 567, 731, 733, 735, 739, 
741, 747, 757, 759, 765, 783, 811, 813, 819, 837, 891, 973, 
975, 981, 999, 1053, 1215, 1459, 1461, 1467, 1485, 1539, 1701, 
 2189, 2191, 2193, 2197, 2199, 2205, 2215, 2217, 2223, 2241, 
2269, 2271, 2277, 2295, 2349, 2431, 2433, 2439, 2457, 2511, 
2673, 2917, 2919, 2925, 2943, 2997, 3159, 3645, 4375, 4377, 
4383, 4401, 4455, 4617, 5103, 6563, 6565, 6567, 6571, 6573, 
6579, 6589, 6591, 6597, 6615, 6643, 6645, 6651, 6669, 6723, 
6805, 6807, 6813, 6831, 6885, 7047, 7291, 7293, 7299, 7317, 
7371, 7533, 8019, 8749, 8751, 8757, 8775, 8829, 8991, 9477, 
10935, 13123, 13125, 13131, 13149, 13203, 13365, 13851, 15309, 
19685, 19687, 19689, 19693, 19695, 19701, 19711, 19713, 19719, 
19737, 19765, 19767, 19773, 19791, 19845, 19927, 19929, 19935, 
19953, 20007, 20169, 20413, 20415, 20421, 20439, 20493, 20655, 
21141, 21871, 21873, 21879, 21897, 21951, 22113, 22599, 24057, 
 26245, 26247, 26253, 26271, 26325, 26487, 26973, 28431, 32805, 
 39367, 39369, 39375, 39393, 39447, 39609, 40095, 41553, 45927, 
59051, 59053, 59055, 59059, 59061, 59067, 59077, 59079, 59085, 
59103, 59131, 59133, 59139, 59157, 59211, 59293, 59295, 59301, 
59319, 59373, 59535, 59779, 59781, 59787, 59805, 59859, 60021, 
60507, 61237, 61239, 61245, 61263, 61317, 61479, 61965, 63423, 
65611, 65613, 65619, 65637, 65691, 65853, 66339, 67797, 72171, 
78733, 78735, 78741, 78759, 78813, 78975, 79461, 80919, 85293, 
98415, 118099, 118101, 118107, 118125, 118179, 118341, 118827, 
120285, 124659, 137781, 177149]

&^?...

Carl Love 28110
June 03 2013
0 0

Joe, Could you explain the significance of the operator &^ that you used? I don't mean explain mathematically, I mean explain its significance in Maple. Is there a package where differentials are expressed like that?

&^?...

Carl Love 28110
June 03 2013

Joe, Could you explain the significance of the operator &^ that you used? I don't mean explain mathematically, I mean explain its significance in Maple. Is there a package where differentials are expressed like that?

map(simplify, vJ) assuming y>0, y<1...

Carl Love 28110
June 02 2013
0 0

@Alejandro Jakubi I decided against making a separate Question, because I now understand what's going on. By mapping simplify, I meant something a bit more elaborate:

map(simplify, J) assuming y>0, y<1;

Here's a worksheet with a simplified example of the situation. Basically, I wanted to know why generic f and Int behave differently with respect to a mapped simplify with assuming. But I understand now.

 

``

> 

e := f(csgn(x)+(x-1)*(x+1), x = 0 .. 1);

f(csgn(x)+(x-1)*(x+1), x = 0 .. 1)

 (1)
> 

simplify(e);

f(csgn(x)+x^2-1, x = 0 .. 1)

 (2)
> 

`assuming`([simplify(e)], [x > 0]);

f(x^2, x = 0 .. 1)

 (3)
> 

map(simplify, e);

f(csgn(x)+x^2-1, x = 0 .. 1)

 (4)
> 

`assuming`([map(simplify, e)], [x > 0]);

f(x^2, x = 0 .. 1)

 (5)

 

I expected all the four of the above simplify results. Below, the only thing that is changed is that f is changed to Int

 

> 

J := Int(csgn(x)+(x-1)*(x+1), x = 0 .. 1);

Int(csgn(x)+(x-1)*(x+1), x = 0 .. 1)

 (6)
> 

simplify(J);

Int(x^2, x = 0 .. 1)

 (7)

 

I understand why that one behaves differently.

> 

`assuming`([simplify(J)], [x > 0]);

Int(x^2, x = 0 .. 1)

 (8)
> 

map(simplify, J);

Int(csgn(x)+x^2-1, x = 0 .. 1)

 (9)
> 

`assuming`([map(simplify, J)], [x > 0]);

Int(csgn(x)+x^2-1, x = 0 .. 1)

 (10)
> 

`assuming`([J], [x > 0]);

Int(csgn(x)+(x-1)*(x+1), x = 0 .. 1)

 (11)
> 

``

``


Download assuming_Int_vs_f.mw

map(simplify, vJ) assuming y>0, y<1...

Carl Love 28110
June 02 2013

@Alejandro Jakubi I decided against making a separate Question, because I now understand what's going on. By mapping simplify, I meant something a bit more elaborate:

map(simplify, J) assuming y>0, y<1;

Here's a worksheet with a simplified example of the situation. Basically, I wanted to know why generic f and Int behave differently with respect to a mapped simplify with assuming. But I understand now.

 

> 

e:= f(csgn(x)+(x-1)*(x+1), x= 0..1);

f(csgn(x)+(x-1)*(x+1), x = 0 .. 1)

> 

simplify(e);

f(csgn(x)+x^2-1, x = 0 .. 1)

> 

simplify(e) assuming x>0;

f(x^2, x = 0 .. 1)

> 

map(simplify, e);

f(csgn(x)+x^2-1, x = 0 .. 1)

> 

map(simplify, e) assuming x>0;

f(x^2, x = 0 .. 1)

I expected all the four of the above simplify results. Below, the only thing that is changed is that f is changed to Int.

> 

J:= Int(csgn(x)+(x-1)*(x+1), x= 0..1);

Int(csgn(x)+(x-1)*(x+1), x = 0 .. 1)

> 

simplify(J);

Int(x^2, x = 0 .. 1)

I understand why that one behaves differently.

> 

simplify(J) assuming x>0;

Int(x^2, x = 0 .. 1)

> 

map(simplify, J);

Int(csgn(x)+x^2-1, x = 0 .. 1)

> 

map(simplify, J) assuming x>0;

Int(csgn(x)+x^2-1, x = 0 .. 1)

> 

J assuming x>0;

Int(csgn(x)+(x-1)*(x+1), x = 0 .. 1)

> 

 

 

Download assuming_Int_vs_f.mw

z >= 0?...

Carl Love 28110
June 02 2013

@only math I think that Preben's and Axel's Answers are making the assumption z >= 0, and that that cannot be derived from the problem as stated (including the Asker's correction). I don't object; it's a natural assumption. I just want to make sure that I am reading their Answers correctly. So, I ask explicitly, are you making an additional assumption that z >= 0?

59051...

Carl Love 28110
June 02 2013

45927 + 4*17*193 = 59051

59049...

Carl Love 28110
June 02 2013

45927 + 2*3^8 = 59049

Multivariate case?...

Carl Love 28110
June 02 2013

I don't know if the multiple-independent-variable case is handled at all. (Are they called partial recurrences, akin to partial differential equations?) There are no such cases in the examples at ?rsolve. Do you have reason to believe otherwise?

CDF better than int(PDF(...), ...)...

Carl Love 28110
June 02 2013
0 0

@Don_Caraota Also note that it is generally better to use CDF(dist,x) rather than evalf(Int(PDF(dist,t), t= -infinity..x)).

CDF better than int(PDF(...), ...)...

Carl Love 28110
June 02 2013

@Don_Caraota Also note that it is generally better to use CDF(dist,x) rather than evalf(Int(PDF(dist,t), t= -infinity..x)).

First 655 656 657 658 659 660 661 Last Page 657 of 710