MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
  • Latest Posts Feed
  • Error when running 'm_model.m' to generate simulink function block when using 'Group all parameters into nested structure'. No errors when using 'Do not group parameters'.

    Error using m_fullSwing (line 125)
    Invalid setting in 'MapleSim_fullSwing/MapleSim_fullSwing/MapleSimParameters' for parameter 'Value'.
    Caused by:
        Error using m_fullSwing (line 125)
        Expression '[Param.Address.ArmPlane, Param.Address.ClubRelHand,  ... ]' for
        parameter 'Value' in 'MapleSim_fullSwing/MapleSim_fullSwing/MapleSimParameters' cannot be evaluated.
            Error using m_fullSwing (line 125)
            Error: Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for
            mismatched delimiters.
     


    Mittag-Leffler function generall simulation for fractional calculus

    alpha := 1:

    beta := 1:

    f1 := sum(x^k/GAMMA(alpha*k+beta), k = 0 .. infinity)

    exp(x)

    (1)

    NULL

    NULL

    NULL

    plot(f1, x = 0 .. 1)

     

    ``


     

    Download Mitag-liffler_Function.mw

    Yesterday, I accidentally discovered a nasty bug in a fairly simple example:

    restart;
    Expr:=a*sin(x)+b*cos(x);
    maximize(Expr, x=0..2*Pi);
    minimize(Expr, x=0..2*Pi);
                                        

    I am sure the correct answers are  sqrt(a^2+b^2)  and  -sqrt(a^2+b^2)  for any real values  a  and  b .  It is easy to prove in many ways. The simplest method does not require any calculations and can be done in the mind. We will consider  Expr  as the scalar product (or the dot product) of two vectors  <a, b>  and  <sin(x), cos(x)>, one of which is a unit vector. Then it is obvious that the maximum of this scalar product is reached if the vectors are codirectional and equals to the length of the first vector, that is, sqrt(a^2+b^2).

    Bugs in these commands were noted by users and earlier (see search by keywords bug, maximize, minimize) but unfortunately are still not fixed. 

    The attached worksheet develops a procedure for extrapolating boundary data from a square to its interior.  Specifically, let's consider the square [−1,1]×[−1,1] and the continuous functions u1, u2, u3, u4 defined on its edges. The procedure constructs a continuous function u(x,y) in the interior of the square which matches the boundary data.

    The function u(x,y) is necessarily discontinuous at a corner if the boundary data of the edges meeting at that corner are inconsistent.  However, if the boundary data are consistent at all corners, then u(x,y) is continuous everywhere including the boundary.

    Here is what the extrapolating function looks like:

    proc (x, y) options operator, arrow; (1/2)*(1-y)*(u__1(x)-(1/2)*(1-x)*u__1(-1))+(1/2)*(x+1)*(u__2(y)-(1/2)*(1-y)*u__2(-1))+(1/2)*(y+1)*(u__3(x)-(1/2)*(x+1)*u__3(1))+(1/2)*(1-x)*(u__4(y)-(1/2)*(y+1)*u__4(1))+(1/4)*(u__1(-1)-u__4(-1))*(-x^2+1)*(1-y)/((x+1)^2+(y+1)^2)^(1/2)+(1/4)*(u__2(-1)-u__1(1))*(-y^2+1)*(x+1)/((y+1)^2+(1-x)^2)^(1/2)+(1/4)*(u__3(1)-u__2(1))*(-x^2+1)*(y+1)/((1-x)^2+(1-y)^2)^(1/2)+(1/4)*(u__4(1)-u__3(-1))*(-y^2+1)*(1-x)/((1-y)^2+(x+1)^2)^(1/2) end proc

     

    The worksheet Extrapolant.mw contains the details of the derivation, and an example.

     

    Maple can solve the easiest two problems of the Putnam Mathematical Competition 2018.  link

     

     

    Problem A1

     

    Find all ordered pairs (a,b) of positive integers for which  1/a + 1/b = 3/2018

     

    eq:= 1/a + 1/b = 3/2018;

    1/a+1/b = 3/2018

    (1)

    isolve(%);

    {a = 1358114, b = 673}

    (2)

    # Unfortunalely Maple fails to find all the solutions; eq must be simplified first!

    (lhs-rhs)(eq);

    1/a+1/b-3/2018

    (3)

    numer(%);

    -3*a*b+2018*a+2018*b

    (4)

    s:=isolve(%);

    {a = -678048, b = 672}, {a = 0, b = 0}, {a = 672, b = -678048}, {a = 673, b = 1358114}, {a = 674, b = 340033}, {a = 1009, b = 2018}, {a = 2018, b = 1009}, {a = 340033, b = 674}, {a = 1358114, b = 673}

    (5)

    remove(u ->(eval(a,u)<=0 or eval(b,u)<=0),[s]);

    [{a = 673, b = 1358114}, {a = 674, b = 340033}, {a = 1009, b = 2018}, {a = 2018, b = 1009}, {a = 340033, b = 674}, {a = 1358114, b = 673}]

    (6)

    # Now it's OK.

     

    Problem B1

     

    Consider the set of  vectors  P = { < a, b> :  0 ≤ a ≤ 2, 0 ≤ b ≤ 100, a, b in Z}.
    Find all v in P such that the set P \ {v} can be partitioned into two sets of equal size and equal sum.

     

    n:=100:
    P:= [seq(seq([a,b],a=0..2), b=0..n)]:

    k:=nops(P): s:=add(P):
    numsols:=0:

    for i to k do
      v:=P[i]; sv:=s-v;
      if irem(sv[1],2)=1 or irem(sv[2],2)=1 then next fi;
      cond:=simplify(add( x[j]*~P[j],j=1..k))-sv/2;
      try
        sol:=[];
        sol:=Optimization:-Minimize
             (0, {x[i]=0, (cond=~0)[], add(x[i],i=1..k)=(k-1)/2 }, assume=binary);
        catch:
      end try:
      if sol<>[] then numsols:=numsols+1;
         print(v='P'[i], select(j -> (eval(x[j],sol[2])=1), {seq(1..k)})) fi;
    od:
    'numsols'=numsols;

    [1, 0] = P[2], {5, 10, 13, 14, 16, 19, 21, 23, 24, 26, 28, 31, 32, 35, 36, 38, 40, 41, 43, 46, 47, 49, 52, 53, 55, 57, 59, 62, 64, 67, 69, 70, 73, 75, 76, 79, 81, 82, 85, 86, 88, 90, 92, 93, 94, 95, 97, 98, 99, 100, 102, 103, 106, 107, 109, 112, 113, 115, 118, 120, 121, 126, 128, 135, 138, 144, 150, 154, 159, 165, 168, 170, 171, 172, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 219, 222, 225, 228, 229, 230, 231, 235, 236, 237, 238, 240, 241, 242, 243, 246, 248, 252, 253, 254, 258, 260, 261, 264, 266, 267, 270, 273, 274, 275, 276, 279, 280, 284, 287}

     

    [1, 2] = P[8], {3, 7, 10, 13, 14, 16, 19, 21, 22, 25, 27, 28, 31, 32, 34, 36, 37, 40, 42, 44, 46, 48, 50, 52, 55, 58, 59, 61, 62, 64, 66, 69, 70, 73, 74, 76, 77, 81, 83, 85, 87, 88, 91, 93, 94, 97, 98, 100, 102, 104, 106, 108, 110, 112, 115, 117, 118, 121, 122, 124, 126, 132, 135, 139, 140, 141, 144, 151, 153, 159, 171, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 215, 216, 217, 218, 219, 220, 222, 223, 225, 229, 230, 231, 234, 236, 237, 238, 240, 241, 242, 243, 246, 249, 250, 252, 255, 257, 258, 259, 261, 263, 266, 267, 269, 273, 276, 279, 281, 282, 284}

     

    [1, 4] = P[14], {7, 10, 13, 15, 19, 21, 22, 24, 26, 28, 31, 32, 33, 37, 38, 40, 42, 44, 47, 48, 50, 52, 55, 56, 59, 60, 63, 64, 66, 68, 71, 72, 75, 76, 78, 81, 82, 87, 88, 90, 91, 94, 95, 96, 97, 100, 102, 103, 105, 106, 109, 110, 113, 114, 115, 116, 117, 118, 121, 122, 124, 127, 132, 135, 138, 142, 144, 145, 154, 156, 159, 160, 163, 165, 170, 171, 172, 176, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 216, 219, 220, 222, 225, 226, 228, 229, 231, 232, 234, 235, 237, 239, 240, 241, 244, 245, 247, 248, 250, 251, 253, 254, 255, 258, 261, 262, 264, 265, 267, 268, 270, 272, 273, 276, 284}

     

    [1, 6] = P[20], {10, 14, 15, 16, 19, 22, 24, 26, 28, 30, 31, 34, 35, 37, 40, 43, 44, 46, 47, 51, 52, 54, 56, 58, 60, 62, 65, 66, 68, 70, 73, 75, 76, 79, 82, 83, 84, 85, 86, 87, 88, 90, 93, 94, 96, 97, 99, 102, 103, 106, 107, 109, 110, 114, 116, 117, 120, 121, 123, 126, 127, 128, 129, 130, 135, 138, 142, 144, 147, 152, 153, 154, 156, 157, 162, 163, 172, 174, 177, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 223, 224, 225, 226, 228, 229, 233, 234, 235, 236, 237, 241, 242, 246, 248, 249, 250, 251, 252, 253, 255, 258, 259, 262, 263, 265, 267, 269, 270, 273, 275, 276}

     

    [1, 8] = P[26], {4, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 29, 30, 34, 36, 37, 38, 40, 43, 44, 46, 48, 51, 54, 55, 56, 58, 61, 63, 64, 66, 70, 71, 73, 74, 78, 79, 80, 82, 85, 86, 88, 90, 92, 94, 97, 100, 102, 103, 106, 108, 109, 112, 113, 115, 118, 120, 121, 123, 129, 133, 135, 141, 144, 150, 158, 159, 160, 163, 168, 169, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 219, 220, 222, 225, 226, 228, 230, 231, 232, 233, 234, 237, 240, 241, 243, 245, 246, 249, 250, 255, 257, 258, 259, 260, 261, 262, 264, 265, 267, 269, 270, 273, 279, 281, 282}

     

    [1, 10] = P[32], {3, 4, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 37, 38, 40, 43, 44, 46, 49, 51, 52, 54, 58, 61, 63, 64, 65, 67, 68, 70, 73, 74, 78, 79, 81, 82, 85, 86, 88, 91, 94, 96, 97, 100, 103, 104, 106, 108, 109, 112, 113, 115, 117, 119, 121, 126, 131, 132, 136, 138, 144, 147, 148, 153, 156, 160, 161, 168, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 217, 219, 220, 222, 224, 225, 226, 227, 228, 231, 233, 234, 236, 240, 242, 243, 245, 246, 248, 249, 251, 252, 254, 255, 257, 258, 261, 262, 264, 266, 267, 269, 270, 273, 275, 276, 279, 282}

     

    [1, 12] = P[38], {7, 10, 14, 16, 19, 22, 25, 28, 29, 31, 34, 35, 36, 37, 40, 43, 46, 47, 48, 49, 52, 54, 57, 59, 61, 63, 64, 67, 68, 70, 72, 75, 76, 79, 81, 82, 85, 87, 88, 90, 93, 94, 96, 98, 99, 100, 101, 103, 106, 108, 109, 111, 112, 114, 117, 118, 121, 122, 124, 132, 135, 141, 150, 151, 158, 159, 162, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 219, 221, 222, 223, 224, 225, 228, 229, 230, 231, 232, 234, 235, 237, 239, 240, 243, 244, 246, 249, 250, 251, 252, 255, 258, 259, 260, 261, 264, 267, 270, 273}

     

    [1, 14] = P[44], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 49, 51, 52, 54, 57, 58, 61, 62, 64, 65, 69, 71, 72, 73, 76, 78, 80, 81, 84, 85, 88, 90, 93, 94, 95, 100, 102, 103, 106, 109, 111, 112, 113, 115, 117, 120, 121, 124, 127, 132, 133, 138, 141, 145, 150, 153, 159, 162, 165, 168, 171, 173, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 224, 225, 226, 227, 228, 229, 230, 234, 235, 236, 238, 241, 242, 243, 244, 247, 248, 249, 252, 256, 258, 259, 260, 261, 264, 267, 268, 270, 272, 273, 275, 276}

     

    [1, 16] = P[50], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 53, 54, 56, 58, 61, 63, 64, 67, 68, 70, 72, 74, 76, 78, 82, 84, 85, 87, 88, 90, 92, 94, 96, 99, 103, 104, 106, 107, 109, 111, 115, 116, 118, 121, 123, 127, 132, 138, 141, 143, 147, 156, 159, 160, 168, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 215, 216, 217, 218, 219, 222, 225, 226, 227, 228, 231, 232, 233, 235, 238, 240, 242, 243, 244, 246, 247, 249, 250, 252, 253, 255, 256, 258, 259, 260, 261, 264, 267, 270, 272, 273, 275, 276, 281}

     

    [1, 18] = P[56], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 60, 62, 64, 66, 68, 70, 73, 75, 76, 78, 80, 82, 85, 86, 88, 89, 93, 94, 96, 98, 100, 102, 105, 106, 109, 110, 113, 114, 116, 118, 119, 121, 123, 129, 132, 136, 144, 148, 153, 168, 169, 170, 171, 172, 173, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 217, 218, 219, 221, 222, 223, 225, 229, 230, 234, 235, 237, 239, 240, 243, 244, 245, 246, 249, 250, 251, 252, 254, 255, 256, 257, 258, 259, 260, 261, 263, 264, 266, 267, 269, 270, 273, 276, 279}

     

    [1, 20] = P[62], {3, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 63, 64, 67, 68, 72, 73, 74, 76, 79, 80, 81, 82, 85, 87, 88, 90, 91, 93, 95, 97, 99, 101, 103, 105, 107, 109, 111, 114, 116, 117, 120, 121, 126, 127, 132, 135, 141, 147, 153, 155, 163, 169, 172, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 228, 229, 232, 234, 235, 237, 238, 240, 243, 244, 246, 248, 249, 253, 254, 255, 258, 261, 262, 264, 267, 268, 270, 273, 274, 275, 276, 279, 282}

     

    [1, 22] = P[68], {12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 73, 74, 75, 76, 78, 79, 82, 83, 84, 85, 88, 90, 91, 93, 94, 96, 98, 100, 103, 104, 109, 111, 112, 114, 115, 117, 119, 121, 124, 129, 132, 136, 138, 141, 142, 150, 153, 156, 159, 160, 168, 169, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 216, 217, 218, 219, 220, 224, 225, 226, 228, 230, 231, 232, 233, 237, 238, 240, 241, 243, 247, 249, 251, 252, 253, 254, 255, 258, 261, 262, 264, 267, 270, 273, 274, 275, 276, 279, 282}

     

    [1, 24] = P[74], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 78, 79, 81, 82, 85, 87, 88, 90, 92, 94, 97, 99, 102, 105, 106, 109, 111, 112, 113, 115, 117, 120, 122, 126, 127, 129, 132, 133, 138, 139, 144, 147, 153, 160, 165, 166, 168, 169, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 225, 227, 228, 229, 230, 231, 232, 233, 234, 237, 239, 240, 244, 245, 246, 247, 249, 251, 252, 253, 255, 258, 261, 262, 264, 267, 270, 273, 275, 276, 277}

     

    [1, 26] = P[80], {3, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 82, 83, 84, 86, 88, 91, 92, 95, 96, 99, 100, 103, 104, 106, 108, 110, 112, 115, 117, 120, 121, 124, 125, 126, 127, 128, 130, 132, 138, 139, 141, 145, 150, 153, 159, 162, 166, 174, 178, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 217, 218, 219, 221, 222, 223, 225, 226, 228, 229, 231, 233, 234, 235, 236, 237, 240, 242, 243, 245, 246, 247, 251, 252, 253, 255, 258, 261, 262, 264, 267, 268, 270, 273, 274, 275, 276, 279, 282}

     

    [1, 28] = P[86], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 90, 91, 92, 94, 97, 98, 101, 102, 104, 106, 109, 111, 112, 115, 116, 119, 120, 123, 124, 129, 132, 134, 136, 141, 147, 153, 160, 162, 173, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 218, 219, 220, 221, 224, 226, 227, 228, 233, 234, 236, 237, 240, 241, 243, 244, 245, 246, 247, 249, 250, 251, 252, 254, 255, 257, 258, 259, 260, 261, 263, 264, 266, 267, 270, 273, 276, 279}

     

    [1, 30] = P[92], {7, 10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 93, 94, 97, 98, 100, 103, 105, 107, 108, 109, 110, 111, 113, 115, 117, 119, 120, 121, 122, 123, 126, 127, 132, 135, 141, 147, 153, 157, 165, 171, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 225, 226, 228, 229, 231, 234, 235, 236, 237, 240, 241, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 262, 264, 267, 268, 270, 273, 275, 276, 279, 281}

     

    [1, 32] = P[98], {10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 102, 105, 106, 110, 111, 112, 113, 114, 118, 120, 121, 123, 124, 126, 127, 128, 129, 130, 135, 138, 144, 150, 153, 156, 157, 160, 162, 166, 167, 171, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 214, 215, 216, 217, 219, 221, 222, 223, 225, 227, 228, 229, 232, 233, 234, 235, 237, 238, 240, 243, 244, 245, 246, 247, 249, 251, 252, 253, 255, 258, 261, 262, 264, 267, 270, 273, 275, 276, 281}

     

    [1, 34] = P[104], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 95, 97, 100, 101, 103, 106, 108, 110, 113, 115, 116, 118, 120, 122, 126, 132, 134, 138, 141, 147, 152, 160, 161, 171, 172, 173, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 217, 219, 220, 222, 223, 225, 226, 228, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 248, 249, 251, 252, 254, 255, 257, 258, 260, 261, 264, 267, 270, 273, 276}

     

    [1, 36] = P[110], {10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 117, 118, 120, 121, 124, 125, 129, 134, 135, 136, 138, 139, 141, 142, 144, 145, 147, 150, 151, 154, 156, 162, 171, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 218, 219, 220, 221, 222, 223, 225, 226, 227, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 264, 267, 270, 273, 275, 276, 279, 281}

     

    [1, 38] = P[116], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 120, 121, 124, 125, 129, 135, 136, 141, 147, 150, 156, 158, 165, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 256, 258, 261, 262, 264, 267, 268, 270, 273, 274, 275, 276, 279, 282}

     

    [1, 40] = P[122], {7, 10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 124, 126, 129, 135, 136, 141, 147, 153, 156, 159, 160, 163, 166, 168, 170, 171, 172, 173, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 242, 243, 244, 245, 246, 247, 249, 251, 252, 255, 258, 261, 264, 267, 270, 273, 275, 276, 279, 281}

     

    [1, 42] = P[128], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 130, 132, 138, 143, 144, 150, 153, 154, 157, 171, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 256, 258, 261, 262, 264, 267, 270, 272, 273, 276, 278, 279, 281, 282}

     

    [1, 44] = P[134], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 138, 140, 144, 150, 159, 162, 168, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 249, 250, 251, 252, 253, 255, 256, 258, 261, 262, 264, 267, 268, 270, 273, 275, 276, 281}

     

    [1, 46] = P[140], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 159, 162, 168, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 249, 250, 251, 252, 253, 255, 257, 258, 261, 262, 264, 267, 268, 270, 273, 275, 276, 277}

     

    [1, 48] = P[146], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 155, 156, 159, 162, 165, 166, 168, 172, 176, 177, 178, 179, 180, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 229, 231, 233, 234, 235, 237, 238, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 262, 264, 267, 269, 270, 273, 275, 276, 281}

     

    [1, 50] = P[152], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 156, 159, 160, 162, 171, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 231, 232, 233, 234, 235, 237, 238, 240, 241, 243, 244, 245, 246, 248, 251, 252, 254, 255, 257, 258, 260, 261, 263, 264, 267, 270, 273, 276}

     

    [1, 52] = P[158], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 156, 157, 160, 163, 165, 171, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 231, 232, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 264, 267, 270, 273, 275, 276, 279, 281, 282}

     

    [1, 54] = P[164], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 156, 157, 160, 162, 163, 172, 173, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 229, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 250, 251, 252, 254, 255, 258, 261, 264, 267, 270, 273, 275, 276, 279, 281, 282}

     

    [1, 56] = P[170], {7, 10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 156, 157, 160, 162, 163, 165, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 251, 252, 255, 258, 261, 264, 266, 267, 270, 273, 275, 276, 279, 281}

     

    [1, 58] = P[176], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 156, 157, 160, 162, 163, 165, 171, 174, 175, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 238, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 264, 267, 270, 273, 275, 276, 279, 280, 282}

     

    [1, 60] = P[182], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 156, 157, 160, 162, 163, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 238, 240, 243, 244, 245, 246, 247, 249, 251, 252, 255, 258, 261, 264, 267, 268, 270, 273, 276, 279, 282, 285}

     

    [1, 62] = P[188], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 156, 168, 171, 172, 173, 174, 175, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 228, 229, 231, 233, 234, 235, 236, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 262, 264, 265, 267, 270, 273, 275, 276, 279, 281}

     

    [1, 64] = P[194], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 153, 156, 157, 160, 162, 163, 171, 173, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 225, 226, 228, 229, 231, 233, 234, 235, 237, 240, 242, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 264, 267, 270, 273, 275, 276, 279, 281}

     

    [1, 66] = P[200], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 132, 137, 141, 144, 150, 153, 157, 168, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 250, 251, 252, 253, 255, 258, 261, 262, 264, 267, 270, 273, 275, 276, 277, 279, 282}

     

    [1, 68] = P[206], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 128, 133, 138, 141, 144, 150, 156, 157, 159, 160, 162, 168, 172, 173, 174, 175, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 218, 219, 220, 221, 222, 223, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 262, 264, 267, 268, 270, 273, 275, 276, 279, 282}

     

    [1, 70] = P[212], {7, 10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 126, 132, 135, 138, 141, 145, 150, 153, 157, 159, 161, 165, 166, 168, 171, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 249, 250, 251, 252, 255, 258, 261, 264, 267, 270, 275, 281}

     

    [1, 72] = P[218], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 112, 114, 118, 119, 121, 122, 124, 126, 132, 133, 138, 141, 144, 150, 159, 166, 169, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 215, 216, 219, 220, 221, 222, 223, 224, 225, 226, 228, 231, 233, 234, 235, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 259, 260, 261, 263, 264, 266, 267, 269, 270, 272, 273, 276, 279}

     

    [1, 74] = P[224], {3, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 101, 103, 105, 107, 109, 113, 114, 116, 117, 121, 122, 124, 126, 129, 135, 141, 143, 147, 152, 153, 160, 165, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 217, 219, 220, 221, 222, 226, 228, 231, 233, 234, 235, 237, 240, 241, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 262, 264, 267, 268, 270, 271, 273, 275, 276, 279, 281, 282}

     

    [1, 76] = P[230], {7, 10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 92, 93, 95, 97, 100, 102, 103, 105, 107, 109, 111, 116, 117, 120, 121, 123, 125, 127, 131, 135, 141, 143, 144, 148, 153, 156, 157, 163, 168, 169, 171, 172, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 228, 229, 231, 233, 234, 235, 236, 237, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 264, 267, 270, 273, 275, 276, 290}

     

    [1, 78] = P[236], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 88, 89, 91, 93, 95, 97, 100, 103, 104, 107, 108, 110, 113, 114, 116, 118, 120, 122, 126, 128, 134, 135, 139, 141, 147, 150, 153, 156, 159, 162, 165, 167, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 214, 215, 216, 217, 218, 219, 220, 223, 225, 226, 228, 229, 231, 232, 233, 234, 235, 237, 238, 240, 243, 244, 245, 246, 247, 251, 252, 253, 255, 258, 261, 262, 264, 267, 268, 270, 273, 275, 276, 281}

     

    [1, 80] = P[242], {7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 75, 79, 81, 82, 84, 85, 87, 90, 91, 96, 97, 99, 100, 102, 104, 107, 109, 110, 112, 114, 117, 120, 121, 125, 126, 127, 130, 131, 138, 141, 144, 150, 152, 153, 154, 159, 163, 166, 174, 175, 176, 177, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 217, 218, 219, 220, 221, 222, 223, 228, 229, 231, 234, 235, 236, 237, 238, 241, 243, 244, 245, 246, 247, 249, 251, 252, 253, 255, 258, 261, 262, 264, 267, 270, 272, 273, 275, 276, 281}

     

    [1, 82] = P[248], {7, 10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 62, 63, 65, 67, 70, 71, 73, 76, 77, 79, 81, 83, 85, 88, 90, 92, 93, 96, 97, 99, 102, 103, 105, 107, 108, 109, 110, 112, 113, 117, 118, 121, 122, 126, 127, 132, 133, 136, 138, 144, 147, 150, 156, 158, 162, 168, 171, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 215, 216, 217, 218, 219, 221, 222, 226, 227, 229, 230, 232, 234, 235, 237, 238, 240, 241, 243, 244, 246, 249, 251, 252, 253, 255, 258, 261, 262, 264, 267, 270, 271, 273, 275, 276, 280, 281}

     

    [1, 84] = P[254], {7, 10, 14, 16, 19, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 51, 53, 55, 58, 59, 61, 63, 65, 68, 69, 71, 72, 73, 74, 76, 79, 81, 82, 85, 86, 91, 93, 94, 96, 97, 99, 101, 103, 108, 109, 110, 111, 113, 115, 118, 120, 124, 127, 132, 133, 135, 141, 144, 147, 148, 150, 151, 153, 154, 156, 159, 162, 165, 168, 169, 171, 172, 174, 175, 176, 177, 178, 179, 180, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 217, 219, 220, 221, 222, 223, 224, 226, 228, 230, 231, 232, 234, 235, 239, 240, 242, 246, 247, 249, 250, 252, 255, 258, 261, 264, 267, 270, 273, 275, 276, 279, 281}

     

    [1, 86] = P[260], {10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 52, 53, 55, 58, 60, 63, 64, 67, 68, 70, 71, 74, 75, 77, 79, 80, 81, 82, 85, 87, 89, 93, 94, 96, 97, 99, 101, 103, 105, 107, 109, 111, 113, 115, 117, 119, 121, 123, 126, 127, 132, 133, 135, 141, 142, 144, 150, 151, 153, 159, 162, 163, 165, 166, 172, 174, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 219, 220, 221, 222, 223, 224, 228, 229, 231, 232, 234, 235, 236, 239, 240, 241, 243, 246, 247, 248, 249, 251, 253, 254, 255, 258, 261, 263, 264, 267, 270, 272, 273, 275, 276, 281}

     

    [1, 88] = P[266], {12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 36, 37, 38, 40, 43, 46, 47, 49, 52, 53, 55, 57, 59, 61, 63, 65, 67, 68, 71, 72, 74, 76, 79, 80, 82, 84, 87, 88, 91, 93, 94, 96, 100, 102, 103, 104, 108, 109, 111, 113, 117, 119, 121, 123, 127, 128, 135, 141, 145, 147, 151, 152, 153, 157, 161, 162, 163, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 203, 204, 205, 206, 207, 208, 209, 210, 213, 214, 215, 216, 217, 220, 222, 223, 225, 228, 229, 231, 233, 237, 238, 239, 240, 244, 246, 247, 252, 253, 254, 255, 258, 261, 263, 264, 267, 268, 269, 270, 273, 275, 276, 281}

     

    [1, 90] = P[272], {10, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 35, 39, 40, 41, 43, 46, 49, 50, 51, 52, 54, 55, 58, 59, 60, 62, 65, 66, 68, 71, 73, 76, 77, 79, 80, 84, 85, 87, 90, 91, 94, 97, 99, 100, 101, 103, 106, 109, 110, 111, 112, 113, 114, 115, 117, 120, 121, 123, 129, 130, 132, 138, 144, 147, 148, 150, 153, 155, 162, 163, 168, 171, 172, 173, 174, 175, 176, 177, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 214, 216, 217, 219, 220, 222, 223, 224, 225, 228, 229, 230, 231, 232, 233, 234, 235, 237, 239, 240, 243, 244, 246, 247, 250, 251, 252, 254, 259, 261, 263, 264, 267, 268, 269, 270, 273, 276}

     

    [1, 92] = P[278], {4, 7, 10, 12, 14, 15, 16, 19, 21, 22, 25, 28, 29, 31, 34, 36, 37, 39, 41, 43, 46, 47, 49, 51, 55, 56, 58, 59, 60, 64, 65, 67, 69, 71, 73, 74, 79, 81, 82, 83, 85, 87, 89, 91, 94, 96, 97, 100, 101, 104, 105, 107, 111, 114, 115, 116, 118, 120, 122, 124, 126, 130, 135, 138, 141, 144, 154, 157, 159, 168, 169, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 214, 216, 217, 219, 221, 223, 224, 225, 226, 228, 230, 231, 233, 237, 239, 240, 241, 243, 244, 246, 247, 249, 252, 253, 254, 255, 259, 261, 263, 264, 267, 270, 273, 275, 276, 279, 281, 282}

     

    [1, 94] = P[284], {3, 4, 7, 10, 12, 14, 15, 16, 19, 22, 24, 25, 28, 30, 31, 34, 35, 36, 39, 41, 43, 46, 47, 49, 52, 53, 55, 58, 59, 61, 63, 66, 67, 70, 71, 73, 76, 78, 79, 82, 84, 85, 88, 90, 93, 95, 97, 99, 101, 103, 106, 108, 110, 112, 115, 118, 120, 121, 123, 127, 129, 135, 138, 139, 144, 145, 150, 156, 157, 159, 166, 168, 171, 174, 175, 176, 177, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 213, 214, 216, 217, 219, 220, 222, 223, 225, 226, 227, 228, 229, 231, 232, 234, 235, 237, 238, 239, 240, 241, 242, 243, 246, 247, 251, 252, 254, 255, 257, 258, 260, 261, 264, 267, 269, 270, 273, 274, 275, 279}

     

    [1, 96] = P[290], {7, 10, 12, 13, 16, 17, 20, 21, 23, 25, 27, 29, 31, 32, 35, 37, 39, 43, 44, 46, 48, 49, 52, 53, 55, 58, 60, 61, 64, 65, 66, 67, 69, 70, 73, 75, 77, 80, 81, 84, 85, 88, 89, 91, 94, 96, 98, 101, 102, 104, 106, 109, 111, 112, 114, 116, 118, 121, 124, 126, 127, 129, 133, 135, 141, 142, 144, 150, 156, 161, 163, 166, 168, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 217, 219, 220, 222, 224, 226, 229, 230, 231, 232, 234, 236, 237, 238, 239, 243, 244, 245, 246, 247, 249, 250, 251, 252, 255, 256, 258, 261, 264, 266, 267, 270, 273, 276, 279, 285}

     

    [1, 98] = P[296], {8, 10, 13, 17, 19, 20, 22, 25, 27, 29, 31, 32, 34, 37, 38, 39, 41, 42, 43, 45, 46, 49, 51, 52, 54, 57, 58, 61, 63, 64, 67, 68, 73, 74, 76, 77, 78, 82, 83, 85, 88, 89, 92, 94, 95, 97, 101, 102, 106, 108, 109, 110, 112, 115, 118, 119, 120, 121, 124, 126, 127, 132, 135, 137, 138, 144, 145, 150, 156, 162, 165, 166, 169, 171, 174, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 225, 226, 227, 228, 230, 231, 232, 233, 234, 236, 237, 240, 241, 245, 246, 247, 248, 249, 252, 253, 254, 258, 260, 264, 267, 269, 270, 273, 276, 279}

     

    [1, 100] = P[302], {6, 8, 10, 13, 15, 17, 19, 20, 22, 25, 27, 29, 31, 32, 34, 37, 38, 39, 41, 42, 43, 45, 46, 49, 51, 52, 54, 57, 58, 61, 63, 64, 67, 68, 71, 72, 74, 76, 78, 80, 82, 84, 87, 88, 91, 92, 94, 96, 98, 100, 102, 104, 106, 109, 112, 113, 115, 117, 119, 121, 126, 128, 132, 136, 144, 150, 157, 165, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 216, 217, 219, 220, 222, 225, 226, 227, 228, 229, 230, 231, 233, 234, 235, 237, 239, 241, 243, 244, 246, 248, 252, 253, 254, 258, 259, 260, 264, 265, 267, 268, 270, 273, 276, 278, 279, 282}

     

    numsols = 51

    (7)

     


    Download putnam2018.mw

     

    Edit.
    Maple can be also very useful in solving the difficult problem B6; it can be reduced to compute the (huge) coefficient of x^1842 in the polynomial

    g := (1 + x + x^2 + x^3 + x^4 + x^5 + x^9)^2018;

    The computation is very fast:

    coeff(g, x, 1842):   evalf(%);
            
    0.8048091229e1147

     

     

    While generating a 3D plot of the solution of an ODE with a parameter, I noticed that better performance could be obtained by calling the plot3d command with a procedure argument, done in a special manner.

    I don't recall this being discussed before, so I'll share it. (It it has been discussed, and this is widely known, then I apologize.)

    I tweaked the initial conditions of the original ODE system, to obtain a non-trivial solution. I don't think that the particular nature of the solution has a bearing on this note.

    restart;

    Digits := 15:

     

     

    The ODE system has two parameters. One, A, will get a fixed value. The

    other, U0, will be used like an independent variable for plotting.

     

     

    eq1:= diff(r[1](t), t, t)+(.3293064114+209.6419478*U0)*(diff(r[1](t), t))
          +569.4324330*r[1](t)-0.3123434112e-2*V(t) = -1.547206836*U0^2*q(t):
    eq2:= 2.03*10^(-8)*(diff(V(t), t))+4.065040650*10^(-11)*V(t)
          +0.3123434112e-2*(diff(r[1](t), t)) = 0:
    eq3 := diff(q(t), t, t)+1047.197552*U0*(q(t)^2-1)*(diff(q(t), t))
           +1.096622713*10^6*U0^2*q(t) = -2822.855019*A*(diff(r[1](t), t, t)):

     

    ics:=r[1](0)=0,D(r[1])(0)=1e-1,V(0)=0,q(0)=0,D(q)(0)=0:

     

    res := dsolve({eq1,eq2,eq3, ics},numeric,output=listprocedure,parameters=[A,U0]):

     

    I will call the procedure returned by dsolve, for evalutions of V(t), as the

    dsolve numeric solution-procedure in the discussion below.

     

    WV := eval(V(t), res):

     

    WV(parameters=[A=1e0]):

     

     

    The goal is to produce a 3D plot of V(t) as a function of t and the parameter U0.

     

     

    tlo,thi := 0.0, 2.0;
    U0lo,U0hi := 1e-3, 2e-1;

    0., 2.0

    0.1e-2, .2

     

    This is the grid size used for plot3d below. It is nothing special.

     

    (m,n) := 51,51;

    51, 51

     

    First, I'll demonstrate that a 3D plot can be produced quickly, by populating a
    Matrix for floating-point evaluations of V(t), depending on t in the first
    Matrix-dimension and on parameter U0 in the second Matrix-dimansion.

     

    The surfdata command is used. This is similar to how plot3d works.

     

    This  computes reasonably quickly.

     

    But generating the numeric values for U0 and t , based on the i,j positions

    in the Matrix, involves the kind of sequence generation formulas that are

    error prone for people.

     

    str := time[real]():
    M:=Matrix(m,n,datatype=float[8]):
    for j from 1 to n do
      u0 := U0lo+(j-1)*(U0hi-U0lo)/(n-1);
      WV(parameters=[U0=u0]);
      for i from 1 to m do
        T := tlo+(i-1)*(thi-tlo)/(m-1);
        try
          M[i,j] := WV(T);
        catch:
          # mostly maxfun complaint for t above some value.
          M[i,j] := Float(undefined);
        end try;
      end do:
    end do:
    plots:-surfdata(M, tlo..thi, U0lo..U0hi,
                    labels=["t","U0","V(t)"]);
    (time[real]()-str)*'seconds';

    1.686*seconds

     

    So let's try it using the plot3d command directly. A 2-parameter procedure
    is constructed, to pass to plot3d. It's not too complicated. This procedure
    will uses one of its numeric arguments to set the ODE's U0 parameter's

    value for the dsolve numeric solution-procedure, and then pass along

    the other numeric argument as a t value.


    It's much slower than the surfdata call above..

     

    VofU0 := proc(T,u0)
           WV(parameters=[U0=u0]);
           WV(T);
         end proc:

    str := time[real]():
    plot3d(VofU0, tlo..thi, U0lo..U0hi,
           grid=[m,n], labels=["t","U0","V(t)"]);
    (time[real]()-str)*'seconds';

    37.502*seconds

     

    One reason why the previous attempt is slow is that the plot3d command

    is changing values for U0 in its outer loop, and changing values of t in its

    inner loop. The consequence is that the value for U0 changes for every

    single evaluation of the plotted procedure. This makes the dsolve numeric

    solution-procedure work harder, by losing/discarding prior numeric
    solution details.

     

    The simple 3D plot below demonstrates that the plot3d command chooses
    x-y pairs by letting its second supplied independent variable be the one
    that changes in its outer loop. Each time the value for y changes the counter

    goes up by one.

     

    glob:=0:
    plot3d(proc(x,y) global glob; glob:=glob+1; end proc,
           0..3, 0..7, grid=[3,3],
           shading=zhue,  labels=["x","y","glob"]);

     

    So now let's try and be clever and call the plot3d command with the two
    independent variables reversed in position (in the call). That will make

    the outer loop change t instead of the ODE parameter U0.

     

    We can use the transform command to swap the two indepenent
    axes in the plot, if we prefer the axes roles switched. Or we could use the
    parametric calling sequence of plot3d for the same effect.

     

    The problem is that this is still much slower!

     

    VofU0rev := proc(u0,T)
           WV(parameters=[U0=u0]);
           WV(T);
         end proc:

    str := time[real]():
    Prev:=plot3d(VofU0rev, U0lo..U0hi, tlo..thi,
                 grid=[n,m], labels=["U0","t","V(t)"]):
    (time[real]()-str)*'seconds';

    plots:-display(
      plottools:-transform((x,y,z)->[y,x,z])(Prev),
      labels=["t","U0","V(t)"],
      orientation=[50,70,0]);

    34.306*seconds

     

    There is something else to adjust, to get the quick timing while using

    the plot3d command here.

     

    It turns out that setting the parameter's numeric value in the
    dsolve numeric solution-procedure causes the loss of previous details
    of the numeric solving, even if the parameter's value is the same.

     

    So calling the dsolve numeric solution-procedure to set the parameter

    value must be avoided, in the case that the new value is the same as

    the old value.

     

    One way to do that is to have the plotted procedure first call the

    dsolve numeric solution-procedure to query the current parameter

    value, so as to not reset the value if it is not changed. Another way

    is to use a local of an appliable module to store the running value

    of the parameter, and check against that. I choose the second way.

     

    And plot3d must still be called with the first independent variable-range

    as denoting the ODE's parameter (instead of the ODE's independent

    variable).

     

    And the resulting plot is fast once more.

     

    VofU0module := module()
           local ModuleApply, paramloc;
           ModuleApply := proc(par,var)
             if not (par::numeric and var::numeric) then
               return 'procname'(args);
             end if;
             if paramloc <> par then
               paramloc := par;
               WV(parameters=[U0=paramloc]);
             end if;
             WV(var);
           end proc:
    end module:

     

    For fun, this time I'll use the parameter calling sequence to flip the

    axes, instead of plots:-transform. That's just because I want t displayed

    on the first axis. But for the performance gain, what matters is that it

    is U0 which gets values from the first axis plotting-range.

     

    str := time[real]():
    plot3d([y,x,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,
           grid=[n,m], labels=["t","U0","V(t)"]);
    (time[real]()-str)*'seconds';

    1.625*seconds

     

    And, naturally, I could also use the parametric form to get a fast plot

    with the axes roles switched.

     

    str := time[real]():
    plot3d([x,y,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,
           grid=[n,m], labels=["U0","t","V(t)"]);
    (time[real]()-str)*'seconds';

    1.533*seconds

     

    Download ode_param_plot.mw

    We have just released updates to Maple and MapleSim, which provide corrections and improvements to product functionality.

    As usual, the Maple update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2018.2.1 download page, where you can also find more details.

    For MapleSim users, use Help>Check for Updates or visit the MapleSim 2018.2.1 download page.

    Just a simple use of DataFrames.  Looking for a new flashlight, I decided to compile a small selection of flashlights for comparison.
     

    interface(rtablesize = 20)

    Type := `<,>`(LED, LED, LED, Incandescent, Incandescent)``

    BType := `<,>`(AAA, AA, AAA, AA, AAA)

    Make := `<,>`(Maglite, Maglite, Maglite, Maglite, Maglite)

    Batteries := `<,>`(1, 2, 3, 2, 1)

    Lumens := `<,>`(47, 245, 200, 14, 2)

    Model := `<,>`(Solitaire, `Pro+`, XL50, mini, Solitaire)

    Minutes := `<,>`(105, 135, 405, 315, 225)

    Cost := `<,>`(17.49, 41, 46.99, 16.50, 10)

    NULL

    Note:The values for cost used is in Canadian dollars, and the stores for the prices taken were from the following - MEC, Lowes, Canadian Tire and Amazon.ca

    NULL

    NULL

    NULL

    Flashlight := DataFrame(`<|>`(Make, Model, Type, BType, Batteries, Lumens, Minutes, Cost), columns = `<,>`("Make", "Model", "Type", "BType", "Batteries", "Lumens", "Minutes", "Cost"))

    _m652460160

    (1)

    NULL

    numelems(Flashlight)

    40

    (2)

    sort(Flashlight, "Cost")

    _m653259008

    (3)

    sort(Flashlight, "Lumens")

    _m629071232

    (4)

    NULL

    No surprise that the cheapest and lowest light outputs were incandescent flashlights.  The list isn't very large so lets add a few more flashlights to our list

    NULL

    new1 := DataFrame(`<,>`(`<|>`(Thrunite, Ti3, LED, AAA, 1, 130, 30, 26.95), `<|>`(Police*Security, Stealth, LED, AA, 1, 80, 60, 7.99), `<|>`(Police*Security, Shield, LED, AA, 1, 120, 90, 15.99), `<|>`(Fenix, E12, LED, AA, 1, 130, 90, 37), `<|>`(Fenix, E20, LED, AA, 2, 265, 30, 50.75), `<|>`(Fenix, E05, LED, AAA, 1, 85, 45, 31.99), `<|>`(Maglite, mini, LED, AAA, 2, 84, 345, 22)), 'columns' = ["Make", "Model", "Type", "BType", "Batteries", "Lumens", "Minutes", "Cost"], 'rows' = [6, 7, 8, 9, 10, 11, 12])

    _m624778752

    (5)

    F1 := Append(Flashlight, new1)

    _m647087136

    (6)

    sort(F1, "Cost")

    _m627084608

    (7)

    F2 := convert(sort(F1, "Cost", `<`), DataFrame)

    _m650479744

    (8)

    Interesting to note that Police Security generally seems to be the cheapest.

     

     

    sort(F2, "BType")

    _m629097504

    (9)

     

     

    Adding yet more flashlights

     

    new2 := DataFrame(`<,>`(`<|>`(Pelican, 1910*Gen3, LED, AAA, 1, 106, 150, 38), `<|>`(Pelican, 1920*Gen3, LED, AAA, 2, 224, 135, 41)), 'columns' = ["Make", "Model", "Type", "BType", "Batteries", "Lumens", "Minutes", "Cost"], 'rows' = [13, 14])

    _m629911872

    (10)

    F3 := Append(F2, new2)

    _m625417184

    (11)

    sort(F3, "Cost")

    _m641386752

    (12)

    NULL

    Now lets select just the Maglite models

     

    F3[`~`[`=`](F3[() .. (), "Make"], Maglite)]

    _m650462944

    (13)

    ``

    Or select the flashlights that have a burn time longer than 100 minutes

     

    F3[`~`[`>`](F3[() .. (), "Minutes"], 100)]

    _m689998912

    (14)

    NULL


     

    Download flashlight3.mw

    From time to time, people ask me about visualizing knots in Maple. There's no formal "Knot Theory" package in Maple per se, but it is certainly possible to generate many different knots using a couple of simple commands. The following shows various examples of knots visualized using the plots:-tubeplot and algcurves:-plot_knot commands.

    The unknot can be defined by the following parametric equations:

     

    x=sin(t)

    y=cos(t)

    z=0

     

    plots:-tubeplot([cos(t),sin(t),0,t=0..2*Pi],
       radius=0.2, axes=none, color="Blue", orientation=[60,60], scaling=constrained, style=surfacecontour);

     

    plots:-tubeplot([cos(t),sin(t),0,t=0..2*Pi],    radius=0.2,axes=none,color=

     

    The trefoil knot can be defined by the following parametric equations:

     

    x = sin(t) + 2*sin(2*t)

    y = cos(t) + 2*sin(2*t)

    z = sin(3*t)

     

    plots:-tubeplot([sin(t)+2*sin(2*t),cos(t)-2*cos(2*t),-sin(3*t), t= 0..2*Pi],
       radius=0.2, axes=none, color="Green", orientation=[90,0], style=surface);

     

    plots:-tubeplot([sin(t)+2*sin(2*t),cos(t)-2*cos(2*t),-sin(3*t),t= 0..2*Pi],    radius=0.2,axes=none,color=

     

    The figure-eight can be defined by the following parametric equations:


    x = (2 + cos(2*t)) * cos(3*t)

    y = (2 + cos(2*t)) * sin(3*t)

    z = sin(4*t)

     

    plots:-tubeplot([(2+cos(2*t))*cos(3*t),(2+cos(2*t))*sin(3*t),sin(4*t),t=0..2*Pi],
       numpoints=100, radius=0.1, axes=none, color="Red", orientation=[75,30,0], style=surface);

     

    plots:-tubeplot([(2+cos(2*t))*cos(3*t),(2+cos(2*t))*sin(3*t),sin(4*t),t=0..2*Pi],    numpoints=100,radius=0.1,axes=none,color=

     

    The Lissajous knot can be defined by the following parametric equations:

     

    x = cos(t*n[x]+phi[x])

    y = cos(t*n[y]+phi[y])

    z = cos(t n[z] + phi[z])

    Where n[x], n[y], and n[z] are integers and the phase shifts phi[x], phi[y], and phi[z] are any real numbers.
    The 8 21 knot ( n[x] = 3, n[y] = 4, and n[z] = 7) appears as follows:
     

    plots:-tubeplot([cos(3*t+Pi/2),cos(4*t+Pi/2),cos(7*t),t=0..2*Pi],
       radius=0.05, axes=none, color="Brown", orientation=[90,0,0], style=surface);

     

    plots:-tubeplot([cos(3*t+Pi/2),cos(4*t+Pi/2),cos(7*t),t=0..2*Pi],    radius=0.05,axes=none,color=

     

    A star knot can be defined by using the following polynomial:
     

    f = -x^5+y^2

     

    f := -x^5+y^2
    algcurves:-plot_knot(f,x,y,epsilon=0.7,
       radius=0.25, tubepoints=10, axes=none, color="Orange", orientation=[60,0], style=surfacecontour);

     

     

    By switching x and y, different visualizations can be generated:

     

    g=(y^3-x^7)*(y^2-x^5)+y^3

     

    g:=(y^3-x^7)*(y^2-x^5)+y^3;
    plots:-display(<
    algcurves:-plot_knot(g,y,x, epsilon=0.8, radius=0.1, axes=none, color="CornflowerBlue", orientation=[75,30,0])|
    algcurves:-plot_knot(g,x,y, epsilon=0.8, radius=0.1, axes=none, color="OrangeRed", orientation=[75,0,0])>);

     

     

    f = (y^3-x^7)*(y^2-x^5)

     

    f:=(y^3-x^7)*(y^2-x^5);
    algcurves:-plot_knot(f,x,y,
      epsilon=0.8, radius=0.1, axes=none, orientation=[35,0,0]);

     

     

    h=(y^3-x^7)*(y^3-x^7+100*x^13)*(y^3-x^7-100*x^13)

     

    h:=(y^3-x^7)*(y^3-x^7+100*x^13)*(y^3-x^7-100*x^13);

    algcurves:-plot_knot(h,x,y,
       epsilon=0.8, numpoints=400, radius=0.03, axes=none, color=["Blue","Red","Green"], orientation=[60,0,0]);

     

    Please feel free to add more of your favourite knot visualizations in the comments below!

    You can interact with the examples or download a copy of these examples from the MapleCloud here: https://maple.cloud/app/5654426890010624/Examples+of+Knots

    In this post an interesting geometric problem is solved: for an arbitrary convex polygon, find a straight line that divides both the area and the perimeter in half. The following results on this problem are well known:
    1. For any convex polygon there is such a straight line.
    2. For any convex polygon into which a circle can be inscribed, in particular for any triangle, the desired line must pass through the center of the inscribed circle.
    3. For a triangle, the number of solutions can be 1, 2, or 3.
    4. If the polygon has symmetry with respect to a point, then any straight line passing through this point is a solution.

    The procedure called  InHalf  (the code below) symbolically solves this problem. The formal parameter of the procedure is the list of coordinates of the vertices of a convex polygon (vertices must be passed opposite or clockwise). The procedure returns all solutions in the form of a list of pairs of points, lying on the perimeter of the polygon, that are the ends of segments that implement the desired dividing.


     

    restart;

    InHalf:=proc(V::listlist)
    local L, n, a, b, M, N, i, j, P, Q, L1, L2, Area, Area1, Area2, Perimeter, Perimeter1, Perimeter2, sol, m, k, Sol;
    uses LinearAlgebra, ListTools;
    L:=map(convert,[V[],V[1]],rational); n:=nops(L)-1;
    a:=<(V[2]-V[1])[1],(V[2]-V[1])[2],0>; b:=<(V[n]-V[1])[1],(V[n]-V[1])[2],0>;
    if is(CrossProduct(a,b)[3]<0) then L:=Reverse(L) fi;
    M:=[seq([L[i],L[i+1]], i=1..n)]:
    N:=0;
    for i from 1 to n-1 do
    for j from i+1 to n do
    P:=map(t->t*(1-s),M[i,1])+map(t->t*s,M[i,2]); Q:=map(s->s*(1-t),M[j,1])+map(s->s*t,M[j,2]);
    L1:=[P,L[i+1..j][],Q,P];
    L2:=[Q,L[j+1..-1][],L[1..i][],P,Q];
    Area:=L->(1/2)*add(L[k, 1]*L[k+1, 2]-L[k, 2]*L[k+1, 1], k = 1 .. nops(L)-1);
    Area1:=Area(L1);
    Area2:=Area(L2);
    Perimeter:=L->add(sqrt((L[k,1]-L[k+1,1])^2+(L[k,2]-L[k+1,2])^2), k=1..nops(L)-2);
    Perimeter1:=Perimeter(L1);
    Perimeter2:=Perimeter(L2);
    sol:=[solve({Area1=Area2,Perimeter1=Perimeter2,s>=0,s<1,t>=0,t<1}, {s,t}, explicit)] assuming real;
    if sol<>[] then m:=nops(sol);
    for k from 1 to m do
    N:=N+1; if nops(sol[k])=2 then Sol[N]:=simplify(eval([P,Q],sol[k])) else Sol[N]:=simplify(eval([P,Q],s=t)) fi;
    od; fi;
    od; od;
    Sol:=convert(Sol, list);
    `if`(indets(Sol)={},Sol,op([Sol,t>=0 and t<1]));
    end proc:  


    Examples of use

    # For the Pythagorean triangle with sides 3, 4, 5, we have a unique solution

    L:=[[4,3],[4,0],[0,0]]:
    P:=InHalf(L);
    plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(P,  color=red), scaling=constrained);

    [[[8/5-(2/5)*6^(1/2), 6/5-(3/10)*6^(1/2)], [4, (1/2)*6^(1/2)]]]

     

     

    # For an isosceles right triangle, there are 3 solutions. We see that all the cuts pass through the center of the inscribed circle

    L:=[[0,0],[4,0],[4,4]]:
    InHalf(L);
    P:=InHalf(L);
    r:=(4+4-4*sqrt(2))/2: a:=4-r: b:=r:
    plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(P,  color=red), plot([r*cos(t)+a,r*sin(t)+b, t=0..2*Pi], color=blue), scaling=constrained);

    [[[2*2^(1/2), 0], [2*2^(1/2), 2*2^(1/2)]], [[4, 0], [2, 2]], [[4, -2*2^(1/2)+4], [-2*2^(1/2)+4, -2*2^(1/2)+4]]]

     

    [[[2*2^(1/2), 0], [2*2^(1/2), 2*2^(1/2)]], [[4, 0], [2, 2]], [[4, -2*2^(1/2)+4], [-2*2^(1/2)+4, -2*2^(1/2)+4]]]

     

     

    # There are 3 solutions for the quadrilateral below

    L:=[[0,0],[4.5,0],[4,3],[0,2]]:
    P:=InHalf(L);
    plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(P,  color=red), scaling=constrained);

    [[[(1/44844)*6^(1/2)*(17^(1/2)-13/2)*37^(3/4)*(2*17^(1/2)+13)^(1/2)*((1836*37^(1/2)-6956)*17^(1/2)+7995*37^(1/2)-56425)^(1/2)-(1/4)*17^(1/2)-(1/8)*37^(1/2)+23/8, 0], [(6^(1/2)*37^(1/4)*((1836*37^(1/2)-6956)*17^(1/2)+7995*37^(1/2)-56425)^(1/2)*(2*17^(1/2)+13)^(1/2)+(-156*37^(1/2)+7770)*17^(1/2)-711*37^(1/2)+50505)/(1776*17^(1/2)+11544), (-6^(1/2)*37^(1/4)*((1836*37^(1/2)-6956)*17^(1/2)+7995*37^(1/2)-56425)^(1/2)*(2*17^(1/2)+13)^(1/2)+(156*37^(1/2)+222)*17^(1/2)+711*37^(1/2)+1443)/(296*17^(1/2)+1924)]], [[(1/90576)*17^(3/4)*(37^(1/2)+37)^(1/2)*(37^(1/2)-37)*((1461*37^(1/2)+29415)*17^(1/2)-986*37^(1/2)-149702)^(1/2)+(5/4)*17^(1/2)+(1/8)*37^(1/2)-27/8, 0], [(17^(1/4)*((1461*37^(1/2)+29415)*17^(1/2)-986*37^(1/2)-149702)^(1/2)*(37^(1/2)+37)^(1/2)+(37*17^(1/2)-51)*37^(1/2)+703*17^(1/2)-1887)/(17*37^(1/2)+629), (17^(1/4)*((1461*37^(1/2)+29415)*17^(1/2)-986*37^(1/2)-149702)^(1/2)*(37^(1/2)+37)^(1/2)+(37*17^(1/2)+85)*37^(1/2)+703*17^(1/2)+3145)/(68*37^(1/2)+2516)]], [[-(1/90576)*17^(3/4)*(37^(1/2)+37)^(1/2)*(37^(1/2)-37)*((1461*37^(1/2)+29415)*17^(1/2)-986*37^(1/2)-149702)^(1/2)+(5/4)*17^(1/2)+(1/8)*37^(1/2)-27/8, 0], [(-17^(1/4)*((1461*37^(1/2)+29415)*17^(1/2)-986*37^(1/2)-149702)^(1/2)*(37^(1/2)+37)^(1/2)+(37*17^(1/2)-51)*37^(1/2)+703*17^(1/2)-1887)/(17*37^(1/2)+629), (-17^(1/4)*((1461*37^(1/2)+29415)*17^(1/2)-986*37^(1/2)-149702)^(1/2)*(37^(1/2)+37)^(1/2)+(37*17^(1/2)+85)*37^(1/2)+703*17^(1/2)+3145)/(68*37^(1/2)+2516)]]]

     

     

    # There are infinitely many solutions for a polygon with a center of symmetry. Any cut through the center solves the problem. The picture shows 2 solutions.

    L:=[[1,0],[1+2*sqrt(3),2],[2*sqrt(3),sqrt(3)+2],[0,sqrt(3)]]:
    P:=InHalf(L);
    plots:-display(plot([L[],L[1]], color=green, thickness=3), plot(eval(P[1],t=1/3),  color=red), scaling=constrained);

    [[[2*3^(1/2)*t+1, 2*t], [-2*3^(1/2)*(t-1), 3^(1/2)-2*t+2]], [[2*3^(1/2)-t+1, 3^(1/2)*t+2], [t, -3^(1/2)*(t-1)]]], 0 <= t and t < 1

     

     

     


     

    Download In_Half.mw

    We have just released an update to Maple, Maple 2018.2. This release includes improvements in a variety of areas, including code edit regions, Workbooks, and Physics, as well as support for macOS 10.14.

    This update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2018.2 download page, where you can also find more details.

    For MapleSim users, the update includes optimizations for handling large models, improvements to model import and export, updates to the hydraulics and pneumatics libraries, and more. For more details and download instructions, visit the MapleSim 2018.2 download page.

    The procedure presented here does independence tests of a contingency table by four methods:

    1. Pearson's chi-squared (equivalent to Statistics:-ChiSquareIndependenceTest),
    2. Yates's continuity correction to Pearson's,
    3. G-chi-squared,
    4. Fisher's exact.

    (All of these have Wikipedia pages. Links are given in the code below.) All computations are done in exact arithmetic. The coup de grace is Fisher's. The first three tests are relatively easy computations and give approximations to the p-value (the probability that the categories are independent), but Fisher's exact test, as its name says, computes it exactly. This requires the generation of all matrices of nonnegative integers that have the same row and column sums as the input matrix, and for each of these matrices computing the product of the factorials of its entries. So, there are relatively few implementations of it, and perhaps none that do it exactly. (Could some with access check Mathematica please?)

    Our own Joe Riel's amazing and fast Iterator package makes this computation considerably easier and faster than it otherwise would've been, and I also found inspiration in his example of recursively counting contingency tables found at ?Iterator,BoundedComposition

    ContingencyTableTests:= proc(
       O::Matrix(nonnegint), #contingency table of observed counts 
       {method::identical(Pearson, Yates, G, Fisher):= 'Pearson'}
    )
    description 
       "Returns p-value for Pearson's (w/ or w/o Yates's continuity correction)" 
       " or G chi-squared or Fisher's exact test."
       " All computations are done in exact arithmetic."
    ;
    option
       author= "Carl Love <carl.j.love@gmail.com>, 27-Oct-2018",
       references= (                                                           #Ref #s:
          "https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test",         #*1
          "https://en.wikipedia.org/wiki/Yates%27s_correction_for_continuity",  #*2
          "https://en.wikipedia.org/wiki/G-test",                               #*3
          "https://en.wikipedia.org/wiki/Fisher%27s_exact_test",                #*4
          "Eric W Weisstein \"Fisher's Exact Test\" _MathWorld_--A Wolfram web resource:"
          " http://mathworld.wolfram.com/FishersExactTest.html"                 #*5
       )
    ;
    uses AT= ArrayTools, St= Statistics, It= Iterator;
    local
       #column & row sums: 
       C:= AT:-AddAlongDimension(O,1), R:= AT:-AddAlongDimension(O,2),
       r:= numelems(R), c:= numelems(C), #counts of rows & columns
       n:= add(R), #number of observations
       #matrix of expected values under null hypothesis (independence):
       #(A 0 entry would mean a 0 row or column in the original, which is not allowed.)
       E:= Matrix((r,c), (i,j)-> R[i]*C[j], datatype= 'positive') / n,
       #Pearson's, Yates's, and G all use a chi-sq statistic, each computed by 
       #slightly different formulae.
       Chi2:= add@~table([
           'Pearson'= (O-> (O-E)^~2 /~ E),                     #see *1
           'Yates'= (O-> (abs~(O - E) -~ 1/2)^~2 /~ E),        #see *2
           'G'= (O-> 2*O*~map(x-> `if`(x=0, 0, ln(x)), O/~E))  #see *3
       ]), 
       row, #alternative rows generated for Fisher's
       Cutoff:= mul(O!~), #cut-off value for less likely matrices
       #Generate recursively all contingency tables whose row and column sums match O.
       #Compute their probabilities under independence. Sum probabilities of all those
       #at most as probable as O. (see *5, *4)
       #Parameters: 
       #   C = column sums remaining to be filled; 
       #   F = product of factorials of entries of contingency table being built;
       #   i = row to be chosen this iteration
       AllCTs:= (C, F, i)->
          if i = r then #Recursion ends; last row is C, the unused portion of column sums. 
             (P-> `if`(P >= Cutoff, 1/P, 0))(F*mul(C!~))
          else
             add(
                thisproc(C - row[], F*mul(row[]!~), i+1), 
                row= It:-BoundedComposition(C, R[i])
             )
          fi      
    ;
       userinfo(1, ContingencyTableTests, "Table of expected values:", print(E));
       if method = 'Fisher' then AllCTs(C, 1, 1)*mul(R!~)*mul(C!~)/n!
       else 1 - St:-CDF(ChiSquare((r-1)*(c-1)), Chi2[method](O)) 
       fi   
    end proc:
    

    The worksheet below contains the code above and one problem solved by the 4 methods


     

     

    DrugTrial:= <
       20, 11, 19;
       4,  4,  17
    >:

    infolevel[ContingencyTableTests]:= 1:

    ContingencyTableTests(DrugTrial, method= Pearson):  % = evalf(%);

    ContingencyTableTests: Table of expected values:

    Matrix(2, 3, {(1, 1) = 16, (1, 2) = 10, (1, 3) = 24, (2, 1) = 8, (2, 2) = 5, (2, 3) = 12})

    exp(-257/80) = 0.4025584775e-1

    #Compare with:
    Statistics:-ChiSquareIndependenceTest(DrugTrial);

    hypothesis = false, criticalvalue = HFloat(5.991464547107979), distribution = ChiSquare(2), pvalue = HFloat(0.04025584774823787), statistic = 6.425000000

    infolevel[ContingencyTableTests]:= 0:
    ContingencyTableTests(DrugTrial, method= Yates):  % = evalf(%);

    exp(-1569/640) = 0.8615885805e-1

    ContingencyTableTests(DrugTrial, method= G):  % = evalf(%);

    exp(-20*ln(5/4)+4*ln(2)-11*ln(11/10)-4*ln(4/5)-19*ln(19/24)-17*ln(17/12)) = 0.3584139703e-1

    CodeTools:-Usage(ContingencyTableTests(DrugTrial, method= Fisher)):  % = evalf(%);

    memory used=0.82MiB, alloc change=0 bytes, cpu time=0ns, real time=5.00ms, gc time=0ns

    747139720973921/15707451356376611 = 0.4756594205e-1

     


     

    Download FishersExact.mw

    Problem:

    Suppose you have a bunch of 2D data points which:

    1. May include points with the same x-value but different y-values; and
    2. May be unevenly-spaced with respect to the x-values.

    How do you clean up the data so that, for instance, you are free to construct a connected data plot, or perform a Discrete Fourier Transform? Please note that Curve Fitting and the Lomb–Scargle Method, respectively, are more effective techniques for these particular applications. Let's start with a simple example for illustration. Consider this Matrix:

    A := < 2, 5; 5, 8; 2, 1; 7, 8; 10, 10; 5, 7 >;

    Consolidate:

    First, sort the rows of the Matrix by the first column, and extract the sorted columns separately:

    P := sort( A[..,1], output=permutation ); # permutation to sort rows by the values in the first column
    U := A[P,1]; # sorted column 1
    V := A[P,2]; # sorted column 2

    We can regard the sorted bunches of distinct values in U as a step in a stair case, and the goal is replace each step with the average of the y-values in V located on each step.

    Second, determine the indices for the first occurrences of values in U, by selecting the indices which give a jump in x-value:

    m := numelems( U );
    K := [ 1, op( select( i -> ( U[i] > U[i-1] ), [ seq( j, j=2..m ) ] ) ), m+1 ];
    n := numelems( K );

    The element m+1 is appended for later convenience. Here, we can quickly define the first column of the consolidated Matrix:

    X1 := U[K[1..-2]];

    Finally, to define the second column of the consolidated Matrix, we take the average of the values in each step, using the indices in K to tell us the ranges of values to consider:

    Y1 := Vector[column]( n-1, i -> add( V[ K[i]..K[i+1]-1 ] ) / ( K[i+1] - K[i] ) );

    Thus, the consolidated Matrix is given by:

    B := < X1 | Y1 >;

    Spread Evenly:

    To spread-out the x-values, we can use a sequence with fixed step size:

    X2 := evalf( Vector[column]( [ seq( X1[1]..X1[-1], (X1[-1]-X1[1])/(m-1) ) ] ) );

    For the y-values, we will interpolate:

    Y2 := CurveFitting:-ArrayInterpolation( X1, Y1, X2, method=linear );

    This gives us a new Matrix, which has both evenly-spaced x-values and consolidated y-values:

    C := < X2 | Y2 >;

    Plot:

    plots:-display( Array( [
            plots:-pointplot( A, view=[0..10,0..10], color=green, symbol=solidcircle, symbolsize=15, title="Original Data", font=[Verdana,15] ),
            plots:-pointplot( B, view=[0..10,0..10], color=red, symbol=solidcircle, symbolsize=15, title="Consolidated Data", font=[Verdana,15] ),
            plots:-pointplot( C, view=[0..10,0..10], color=blue, symbol=solidcircle, symbolsize=15, title="Spread-Out Data", font=[Verdana,15] )
    ] ) );

    Sample Data with Noise:

    For another example, let’s take data points from a logistic curve, and add some noise:

    # Noise generators
    f := 0.5 * rand( -1..1 ):
    g := ( 100 - rand( -15..15 ) ) / 100:
    
    # Actual x-values
    X := [ seq( i/2, i=1..20 ) ];
    
    # Actual y-values
    Y := evalf( map( x -> 4 / ( 1 + 3 * exp(-x) ), X ) );
    
    # Matrix of points with noise
    A := Matrix( zip( (x,y) -> [x,y], map( x -> x + f(), X ), map( y -> g() * y, Y ) ) );

    Using the method outlined above, and the general procedures defined below, define:

    B := ConsolidatedMatrix( A );
    C := EquallySpaced( B, 21, method=linear );

    Visually:

    plots:-display( Array( [
        plots:-pointplot( A, view=[0..10,0..5], symbol=solidcircle, symbolsize=15, color=green, title="Original Data", font=[Verdana,15] ),
        plots:-pointplot( B, view=[0..10,0..5], symbol=solidcircle, symbolsize=15, color=red, title="Consolidated Data", font=[Verdana,15]  ),
        plots:-pointplot( C, view=[0..10,0..5], symbol=solidcircle, symbolsize=15, color=blue, title="Spread-Out Data", font=[Verdana,15] )
    ] ) );

      

    Generalization:

    Below are more generalized custom procedures, which are used in the above example. These also account for special cases.

    # Takes a matrix with two columns, and returns a new matrix where the new x-values are unique and sorted,
    # and each new y-value is the average of the old y-values corresponding to the x-value.
    ConsolidatedMatrix := proc( A :: 'Matrix'(..,2), $ )
    
            local i, j, K, m, n, P, r, U, V, X, Y:
      
            # The number of rows in the original matrix.
            r := LinearAlgebra:-RowDimension( A ):
    
            # Return the original matrix should it only have one row.
            if r = 1 then
                   return A:
            end if:
    
            # Permutation to sort first column of A.
            P := sort( A[..,1], ':-output'=permutation ):       
    
            # Sorted first column of A.
            U := A[P,1]:
    
            # Corresponding new second column of A.
            V := A[P,2]:
    
            # Return the sorted matrix should all the x-values be distinct.
            if numelems( convert( U, ':-set' ) ) = r then
                   return < U | V >:
            end if:
    
            # Indices of first occurrences for values in U. The element m+1 is appended for convenience.
            m := numelems( U ):
            K := [ 1, op( select( i -> ( U[i] > U[i-1] ), [ seq( j, j=2..m ) ] ) ), m+1 ]:
            n := numelems( K ):
    
            # Consolidated first column.
            X := U[K[1..-2]]:
    
            # Determine the consolidated second column, using the average y-value.
            Y := Vector[':-column']( n-1, i -> add( V[ K[i]..K[i+1]-1 ] ) / ( K[i+1] - K[i] ) ):
    
            return < X | Y >:
    
    end proc:
    
    # Procedure which takes a matrix with two columns, and returns a new matrix of specified number of rows
    # with equally-spaced x-values, and interpolated y-values.
    # It accepts options that can be passed to ArrayInterpolation().
    EquallySpaced := proc( M :: 'Matrix'(..,2), m :: posint )
    
            local A, i, r, U, V, X, Y:
    
            # Consolidated matrix, the corresponding number of rows, and the columns.
            A := ConsolidatedMatrix( M ):
            r := LinearAlgebra:-RowDimension( A ):
            U, V := evalf( A[..,1] ), evalf( A[..,2] ):
    
            # If the consolidated matrix has only one row, return it.
            if r = 1 then
                   return A:
            end if:
    
            # If m = 1, i.e. only one equally-spaced point is requested, then return a matrix of the averages.
            if m = 1 then
                   return 1/r * Matrix( [ [ add( U ), add( V ) ] ] ):
            end if:
    
            # Equally-spaced x-values.
            X := Vector[':-column']( [ seq( U[1]..U[-1], (U[-1]-U[1])/(m-1), i=1..m ) ] ):
    
            # Interpolated y-values.
            Y := CurveFitting:-ArrayInterpolation( U, V, X, _rest ):    
    
            return < X | Y >:
    
    end proc:

    Worth Checking Out:

     

    The attached worksheet shows how to evaluate and graphically analyze an autonomous first-order nonlinear recurrence with two dependent variables and multiple symbolic parameters. 

    This worksheet shows how a small module that simply encapsulates the given information of a problem combined with some use statements can greatly facilitate the organization of one's work, can encapsulate the setting of parameter values, and can allow one to work with symbolic parameters.

    Edit: In the first version of this Post, I forgot to include the qualifier "autonomous".  The system being autonomous substantially simplifies its treatment.
     

    Autonomous first-order nonlinear recurrences with parameters and multiple dependent variables

    Author: Carl Love <carl.j.love@gmail.com> 20-Oct-2018

     

    The techniques used in this worksheet can be applied to most autonomous first-order nonlinear recurrences with multiple dependent variables and parameters.

     

    This worksheet shows how a small module that simply encapsulates the given information of a problem combined with some use statements

    • 

    can greatly facilitate the organization of one's work,

    • 

    can encapsulate the setting of parameter values,

    • 

    can allow one to work with symbolic parameters.

     

    A Problem from MaplePrimes: A discrete Lottka-Volterra population model is applied to an isolated island with a population of predators (foxes), R, and prey (rabbits), K. [Note that R is the foxes, not the rabbits! Perhaps this problem statement originated in another language.] The change over one time period is given by

    K[n+1]:= K[n]*(-b*R[n]+a+1);  R[n+1]:= R[n]*(b*e*K[n]-c+1),

    where a, b, c, e are parameters of the model. In this problem we will use a= 0.15, b= 0.01, c= 0.02, e= 0.01, when numeric values are needed.

     

    a) Show that there exists an equilibrium (values of K[n] and R[n] such that K[n+1] = K[n] and R[n+1] = R[n]).

     

    b) Write Maple code that solves the recurrence numerically. Assume that if any population is less than 0.5 then it has gone extinct and set the value to 0. Check that your program is idempotent at the equilibrium.

     

    restart:

    We begin by collecting all the given information (except for specific numeric values) into a module. The ModuleApply lets the user set the numeric values later.

     

    For all two-element vectors used in this worksheet, K is the first value and R is the second value.

    KandR:= module()
    local
       a, b, c, e, #parameters

       #procedure that lets user set parameter values:
       ModuleApply:= proc({
           a::algebraic:= KandR:-a, b::algebraic:= KandR:-b,
           c::algebraic:= KandR:-c, e::algebraic:= KandR:-e
       })
       local k;
          for k to _noptions do thismodule[lhs(_options[k])]:= rhs(_options[k]) od;
          return
       end proc,

       Extinct:= (x::realcons)-> `if`(x < 0.5, 0, x) #force small, insignificant values to 0
    ;
    export
       #Procedure that does one symbolic iteration
       #(Note that this procedure uses Vector input and output.)
       iter_symb:= KR-> KR *~ <-b*KR[2]+a+1, b*e*KR[1]-c+1>, 

       #Such simple treatment as above is only possible for autonomous
       #recurrences.

      
       iter_num:= Extinct~@iter_symb #one numeric iteration
    ;
    end module:

    #The following expression is the discrete equivalent of the derivative (or gradient).
    #It represents the change over one time period.
    P:= <K,R>:  
    OneStep:= KandR:-iter_symb(P) - P

    Vector(2, {(1) = K*(-R*b+a+1)-K, (2) = R*(K*b*e-c+1)-R})

    #An equilibrium occurs when the gradient is 0.
    Eq:= <K__e, R__e>:
    Eqs:= solve({seq(eval(OneStep=~ 0, [seq(P=~ Eq)]))}, [seq(Eq)]);

    [[K__e = 0, R__e = 0], [K__e = c/(b*e), R__e = a/b]]

    #We're only interested here in nonzero solutions.
    EqSol:= remove(S-> 0 in rhs~(S), Eqs)[];

    [K__e = c/(b*e), R__e = a/b]

    #Set parameters:
    KandR(a= 0.15, b= 0.01, c= 0.02, e= 0.01);

    #Show idempotency at equilibrium:
    use KandR in Eq0:= eval(Eq, EqSol); print(Eq0 = iter_num(Eq0)) end use:

    (Vector(2, {(1) = 200.0000000, (2) = 15.00000000})) = (Vector(2, {(1) = 200.0000000, (2) = 15.00000000}))

    #procedure that fills a Matrix with computed values of a 1st-order recurrence.
    #(A more-efficient method than this can be used for linear recurrences.)
    #This procedure has no dependence on the module.
    Iterate:= proc(n::nonnegint, iter, init::Vector[column])
    local M:= Matrix((n+1, numelems(init)), init^+, datatype= hfloat), i;
       for i to n do M[i+1,..]:= iter(M[i,..]) od;
       M
    end proc:

    We want to see what happens if the initial conditions deviate slightly from the equilibrium. It turns out that any deviation (as long as the
    initial values are still nonnegative!) will cause the same effect. I simply chose the deviation <7,2> because it was the smallest for which

    the plot clearly showed what happens using the scale that I wanted to show the plot at. By using a finer scale, it is possible to see the

    "outward spiral" efffect from even the tiniest deviation.

    dev:= <7,2>:
    use KandR in KR:= Iterate(1000, iter_num, Eq0 + dev) end use:

    plot(
       [
           KR, #trajectory of population
           KR[[1,1],..], #1st point
           KR[-[1,1],..], #last point,
           <Eq0|Eq0>^+, #equilibrium
           #every 100th point (helps show time scale):
           KR[100*[$1..iquo(numelems(KR[..,1]), 100)-1], ..]
       ],
       #This group of options are all lists, each element of which corresponds
       #to one of the above components of the plot:
       style= [line, point$4],
       symbol= [solidcircle$4, soliddiamond],
       symbolsize= [18$4, 12],
       color= [black, green, red, brown, blue],
       thickness= [0$5],
       legend= [`pop.`, init, final, equilibrium, `100 periods`],

       #This group of options are lists, each element of which corresponds to one
       #coordinate axis (horizontal, then vertical).
       view= [0..max(KR[..,1]), 0..max(KR(..,2))],
       labels= [rabbits, foxes],
       labeldirections= [horizontal,vertical],
       size= [700,700], #measured in pixels

       #options applied to whole plot:
       labelfont= [TIMES, BOLDITALIC, 14],
       title= "Population of foxes and rabbits over time" "\n", titlefont= [TIMES,16],
       caption=
          "\n" "Choosing an initial point near the equilibrium causes"
          "\n" "outward spiraling divergence." "\n",
       gridlines= false
    );
     

    A fieldplot helps show what happens for any starting values. An arrow is drawn from each of a 2-D grid of point. The magnitude and direction of the arrow show the gradient (as a vector) in this case.

    plots:-fieldplot(
       rtable_eval(OneStep),
       K= 0..max(KR[..,1]),  R= 0..max(KR[..,2]), grid= [16,16],

       #arrow-specific options:
       anchor= tail, fieldstrength= log, arrows= slim, color= "DarkGreen",

       #other options (same as any 2D plot):
       labels= [rabbits, foxes], labeldirections= [horizontal,vertical],
       labelfont= [TIMES, BOLDITALIC, 14],
       title= "One-step population changes from any point" "\n", titlefont= [TIMES,16],
       caption= "\n" "All trajectories spiral outward from the equilibrium." "\n",
       size= [700,700],
       gridlines= false
    );

    The above plot is computed only from the symbolic discrete gradient expression OneStep; it does not use the computed population values from the first plot. It only uses the maxima of those computed values to determine the length of the axes.

     

    Conclusion: While this is interesting stuff mathematically, and makes for great plots, divergence from the equilibrium doesn't seem realistic to me.

     


     

    Download FoxesAndRabbits.mw

    With this application, you visualize the DNA chain using the position vector as a fundamental tool to describe its curvature and radius of curvature. I also show the native maple syntax for the graphics. You can download the maple center app to show different DNA trajectories based on the position vector. Developed for students of health sciences.
     

    k_adn.zip  (In spanish)

    https://www.maplesoft.com/applications/download.aspx?SF=154486/Plot_of_curvature_and_radius_of_curvature.mw

    Lenin Araujo Castillo

    Ambassador of Maple

    1 2 3 4 5 6 7 Last Page 1 of 263