Carl Love

Carl Love

28050 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Given any table, you can construct its setwise inverse table with this procedure:

TableInverse:= (T::table)->  (([lhs]~)~@ListTools:-Classify)(rhs, [entries](T, ':-pairs')):

Example of use:

T:= table():  R:= rand(0..5):  for i to 9 do for j to 9 do T[i,j]:= 'R'()$2 od od:

TI:= TableInverse(T):
T[6,7];

                              1, 2

TI[1,2];
                    {[2, 5], [6, 7], [7, 2]}

If the table is one-to-one (aka injective), or you want to ignore multiple indices mapped to the same entry, then a simpler procedure can be used:

TableInverse1to1:= (T::table)-> (table@(rhs=lhs)~@[entries])(T, ':-pairs'):

The time complexity of TableInverse is O(n*) where n is the number of indices, and the * represents unknown but small (< < n) logarithmic factors. The time complexity of TableInverse1to1 is simply O(n).

You can teach diff a recursive formula for the derivative of P(n,x). Once you do that, it'll be able to automatically do higher-order derivatives and/or derivatives of compositions (such as with cos(theta)).

Important: If you want to do this, do not use with(orthopoly). Instead, if and when you want to eliminate P from an expression ex, do eval(ex, P= orthopoly:-P).

`diff/P`:= proc(n,x,z)
    if depends(n,z) then 
        error "derivative of P wrt its 1st argument not defined"
    else
        diff(x,z)*n/(x^2-1)*(x*'P'(n,x) - 'P'(n-1, x))
    fi
end proc
:
#Optional (only affects display): Display P(n,x) with the n as a 
#subscript:
`print/P`:= (n,x)-> nprintf(`#msub(mi("P"),mo("%a"))`, n)(x)
:
dP:= diff(P(n, cos(theta)), theta$2);
simplify(eval(dP, n=2));
simplify(eval(%, P= orthopoly:-P));

 

I don't know what version of Maple you are using. I am using Maple 2020. So, I don't know if this Answer will work for you. It's useful to make assumptions on x and alpha. If need be, the system can internally derive an appropriate assumption for t based on those. It may be counter-productive to make an assumption on t

This works instantaneously for me:

b:= GAMMA(2-alpha)/(1-alpha)/GAMMA(1-alpha):
J:= Int((x-t)^(-alpha)*a*(t-b*ln(1+t/b)), t= 0..x);
value(J) assuming x > 0, alpha > 0, alpha < 1;

If you need to resort to numerical integration, there are appropriate easy ways to do that in Maple. The method that you show is not very good.
 

Here's another way:

numboccur(evalb~(A >~ B), true)

@mthkvv As far as I know, it's not possible to programmatically set the fps (frames per second) with any of the animation commands: plots:-animateplots:-display, or Explore. You must set it from the toolbar (or context menu) animation controls, which only appear when the animation is the current on-screen "context".

If you save the animation to a .GIF file, then you can change the fps using numerous freely downloadable GIF-editing apps.

An often-used workaround is to repeat frames. So, to change from the default 10 fps to x fps, repeat each frame round(10/x) times. However, if you do too much of this, the GUI can become overwhelmed by the amount of data and become unbearably sluggish. It largely depends on the amount of numerical data in each frame.

Change to

f:= unapply(convert(g(x), surd), x);

I'm not 100% sure that this will work with in Maple 13. Let me know. If it doesn't work, I can give you something else.

But regardless of whether it works with or not, you should change f as shown. Here's a rule to help decide between x-> ... and unapply(..., x): If the right side of the arrow (in this case convert(g(x), surd)) contains a symbolic operation that doesn't depend on the specific value(s) of the variable(s) on the left side of the arrow (in this case x)---it only depends on x as a symbol---then use unapply. The convert(g(x), surd) is a purely symbolic operation: It has no dependence on the numeric value of x.

If you fail to follow this rule, the result will still often "work" in a sense, but if it does work, it'll always be extremely inefficient. The most-common violations of this rule involve putting diff or int on the right side of the arrow (in which case the result usually won't work with numeric values). (I'm not saying that you should never put those on right side of the arrow; you need to consider the situation from the symbolic perspective explained in the previous paragraph.)

This functional-compostion diagram illustrates the above. I copied it from help page ?operators,D:

       function application -->
 f   -----------------------------   f(x)
           <-- unapply
 |                                    |
 |                                  d | ^
 |                                  i | |
 |D()                               f |
 | |                                f | i
 | |                                  | n
 | V                                | | t
 |                                  V |
 |                                    |
       function application -->
D(f) ----------------------------   f'(x)
           <-- unapply

 

Note that you've erroneously put a space after GAMMA. It must be, for example, GAMMA(alpha+1) instead of GAMMA (alpha+1). This type of error is only possible in 2-D Input.

A vertically stacked array of three plots can be made by

restart:
F:= (c::list(numeric))-> 
    alpha-> 
        t-> add(c[i]*(t^alpha)^(i-1)/GAMMA((i-1)*alpha+1), i= 1..nops(c))
:
Alpha:= [0.2, 0.6, 0.8, 1.0]:
tr:= 0..1:
c:= Record(
    "R"= [0.46, 9.1625, 8.8318, 11.6888],
    "L"= [0.32, 0.09282, 2.1126. 3.9028],
    "T"= [0.52, 0.0569, 0.0243, 1.3102]
):
plots:-display(
    <seq(
        plot(
            F(c[x])~(Alpha), tr, legend= (alpha=~ Alpha),
            title= typeset(x, " for various ", alpha),
            labels= [t, x]
        ),
        x= exports(c)
    )>,
    view= [tr, DEFAULT]       
);

If you want all the plots to have the same scale on the vertical axis, change DEFAULT to
0..max(seq(max(F(c[x])~(Alpha)(rhs(tr))), x= exports(c)))

The plots can also be horizontally arrayed.

SumNeighborhood:= (G::GRAPHLN, v)->
    add(GraphTheory:-Degree~(G, GraphTheory:-Neighborhood(G, v)))
:

Make it a list: numelems([b]);

My procedure for this uses neither a remember table nor recursion. It does have something like a "reverse remember table". The sequence values are stored internally as hardware floats for faster computation.

This procedure is about 29 times faster than Kitonum's procedure at generating the first 2^13 entries of the sequence.

A:= proc(N::And(posint, Not({1,2})))
local 
    n, m, p, U:= table([HFloat(1)= 1]), 
    A:= rtable(1..N, [1,0], datatype= float[8])
;
    for n from 2 to N-1 do
        p:= A[n];
        m:= U[p]; 
        A[n+1]:= `if`(m::posint, n-m, min(abs~(p -~ A[..n-1])));
        U[p]:= n       
    od;
    seq(trunc(a), a= A)
end proc
:
S:= CodeTools:-Usage([A(2^13)]):
memory used=151.81MiB, alloc change=27.49MiB, 
cpu time=453.00ms, real time=460.00ms, gc time=0ns

S[..30];
[1, 0, 1, 2, 1, 2, 2, 1, 3, 1, 2, 4, 1, 3, 5, 1, 3, 3, 1, 3, 2, 
  10, 5, 8, 2, 4, 14, 4, 2, 4]

 

Simply add the option view= [0..24, 0..100].

This bug is very easy to work around by passing the color to display:

plots:-display(
    plots:-matrixplot(A, heights= histogram), 
    color= "X11 Thistle1"
);

There are numerous plotting commands (I've especially noticed them in package Statistics) with similar shortcomings, which may pertain to color or other options, that can be worked around in the same way. In some cases, it's necessary to also include overrideoption as an option to display.

I don't think that the Answers presented thus far can truly be called "mapping" in the fullest sense of that word because they don't produce a Matrix as output. Here's a solution that does. First, I define a procedure analogous to map:

DimMap:= proc(f, A::rtable, d::posint:= 1)
local i, D:= rtable_dims(A), R:= rtable(A);
    for i from lhs(D[d]) to rhs(D[d]) do
        R[D[..d-1], i, D[d+1..]]:= f(R[D[..d-1], i, D[d+1..]], _rest)
    od;
    R
end proc
:

The 3rd argument is the dimension over which to operate: 1 for rows, 2 for columns, and higher-dimensional Arrays are handled also.

Usage is just like map:

Sx:=1/sqrt(2)*Matrix([[0,1,0],[1,0,1],[0,1,0]]);
lam,v:=LinearAlgebra:-Eigenvectors(Sx);

DimMap(LinearAlgebra:-Normalize, v, 2, 2);

The 3rd argument, 2, is the dimension to map over. The 4th argument, 2, is the extra argument to Normalize.

Regarding your code: Matrices should be indexed with square brackets unless you're sure that you know what you're doing with the round brackets. I recall about two months ago that you had a surprising issue caused by using round brackets. I thought that you learned your lesson then. Round brackets aren't simply an alternative syntax; their semantics are very different (especially for Arrays).

Although VV's method is fast and exact, a further substantial time improvement can be made by integerizing the matrix of rational polynomials. This is done by multiplying it by its lowest common denominator. Also, the technique below can return the LU factorization of A, which'll be beneficial if you'd like to use different bs with the same A.

#This procedure returns a factorization of A, A = d*A1, where d 
#is rational and A1 is a matrix of polynomials with integer 
#coefficients.
IntPolyMatrix:= proc(A::rtable)
local
    x, 
    AR:= convert(A, rational),
    d:= ilcm(seq(x, x= denom~(AR))),
    V:= indets(A, name)
;
    1/d, rtable(d*AR, datatype= `if`(V={}, integer, polynom(integer, V)))
end proc
:
#This procedure returns a factorization of A, A = d*P.L.U, where d is
#rational and P, L, and U are the standard matrices returned by
#LUDecomposition.
LUIntPolyMatrix:= proc(A::Matrix)
local A_denom, A_int;
    (A_denom, A_int):= IntPolyMatrix(A);
    (    
        A_denom, 
        LinearAlgebra:-LUDecomposition(
            A_int, 
            ':-method'= ':-FractionFree', 
            ':-output'= [':-P', ':-L', ':-U']
        )
    )
end proc
:
#This procedure returns the exact matrix/vector X such that A.X = b.
#X is a matrix of rational functions with integer coefficients.
SolveIntPolyMatrix:= proc(A::Matrix, b::{Vector, Matrix})
local b_denom, b_int, dPLU;
    dPLU:= LUIntPolyMatrix(A);
    (b_denom, b_int):= IntPolyMatrix(b); 
    b_denom/dPLU[1]*LinearAlgebra:-LinearSolve([dPLU[2..]], b_int)
end proc
:
#Time test (using the A and b from your worksheet):
X:= CodeTools:-Usage(SolveIntPolyMatrix(A,b)):
memory used=277.44MiB, alloc change=0 bytes, 
cpu time=1.47s, real time=1.36s, gc time=250.00ms

#Verify correctness:
[seq](x, x= simplify(convert(A, rational).X - convert(b, rational)));
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 0, 0, 0, 0]

 

Yes, it's easy. To do it completely, there are 2-1/2 or 3 steps:

1. Use the command LibraryTools:-Save(x, y, z, ..., filename) to create a library named filename and add x, y, and z to it. The same command can be used to add new names to an existing library. Anything that can be assigned to a name can be saved; they don't need to be procedures.

2. The global variable libname determines which libraries are visible and the order in which they are accessed (very important if the same name is used in multiple libraries). libname contains the names of the directories that Maple will look in for the libraries.

3. To make it fully automatic as you describe, you must update the value of libname in your initialization file. Make sure that you update it, not just set it. For example,

libname:= libname, "the new directory":
or
libname:= "the new directory", libname:

depending on the access order you want in case of duplicate names.

 

First 70 71 72 73 74 75 76 Last Page 72 of 395