Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

In addition to Preben's Answer: There's an extra constant, ga, in the fourth equation, for which no value is provided. If you provide a value for this, then it only expects 9 boundary conditions.

Sorry this took a few days. It was a bit more complicated than I anticipated. The main complexity was that Graphs need to have unqiue labels for each vertex, but, of course, we may want to repeat labels in the tree.

I welcome any critiques of my module. Especially, I'd like to know if there are any modern aspects modules that I could have used.

(* Written by Carl Love, 16-Mar-2013.

 A module to convert a function tree representation (such as is returned by ToInert
 or various commands in XMLtools) into a tree represented in the GraphTheory package
                                                                                     *)
FunctionToTree:= module()
local
(*____________________________________________________________________________________
Phase 1: Build a representation of the tree in tables with each vertex represented by
by an integer and each edge by an ordered pair of integers.                        
``````````````````````````````````````````````````````````````````````````````````````*)
     # Vertices is a table that maps the original functional representation of the
     # node to an integer.  VertexLabels is simply the inverse of that table.
     Vertices::table, VertexLabels::table, Vertex_index::nonnegint,
     
    # Edges is a set (stored as a table indexed by an integer for efficiency)
    # of ordered pairs of vertices in their integer representation.
    Edges::table, Edge_index::nonnegint,
    
    AddEdge:= proc(e::list(function))
        Edge_index:= Edge_index + 1;  
        Edges[Edge_index]:= [Vertices[e[1]], Vertices[e[2]]]
    end proc,

    AddVertex:= proc(x::function)
        Vertex_index:= Vertex_index + 1;
        Vertices[x]:= Vertex_index;
        VertexLabels[Vertex_index]:= x
    end proc,

    # Recursive
    AddSubTree:= proc(f::function(function))
    local x::function;
        for x in f do
            #Make functions unique.
            if assigned(Vertices[x]) then  x:= convert(op(0,x), `local`)(op(x))  fi;
            AddVertex(x);
            #Leaves are type function but not function(function).
            if x::function(function) then thisproc(x) end if;
            AddEdge([f,x])
        end do
    end proc,    
(*__________________________________________________________________________________
Phase 2: Shorten the function labels to something that can be shown on a plot of the
tree but which is still meaningful.
`````````````````````````````````````````````````````````````````````````````````````*)
    # Prefix is the number of chars at the beginning of each function name that
     # can be trimmed off.  For example, if every function name begins _Inert_, then
     # Prefix should be set to 7.  This value is passed in by the user.
     Prefix::nonnegint,

     # Labels is a table mapping the string form of the shortened function names to
     # alternative (usually abbreviated) representations.  For example, "SUM" can be
     # mapped to `+`.  This is passed in by the user.
     Labels::table,

     #Shorten the function names
     StripName:= proc(f::function)
       local f0:= sprintf("%a", op(0,f))[Prefix+1 .. -1];
           `if`(assigned(Labels[f0]), Labels[f0], nprintf("%s", f0))
    end proc,

    ShortenLabels:= proc()
    local k::nonnegint, f::function, f0::symbol;
        VertexMap:= table();
        for k to Vertex_index do
            f:= VertexLabels[k];
            f0:= StripName(f);
            if not f::function(function) then
                #For leaves, display the function operands also.
                f0:= nprintf("%Q", `if`(f0 = ``, op(f), f0(op(f))))
            end if;
            if assigned(VertexMap[f0]) then f0:= convert(f0, `local`) fi;
            VertexMap[f0]:= f;
            VertexLabels[k]:= f0
        end do             
    end proc,
(*____________________________________________________________________________________
                                         Main                                               
``````````````````````````````````````````````````````````````````````````````````````*)
    ModuleInits:= proc()
        Edges:= table();
        Edge_index:= 0;
        Vertices:= table();
        VertexLabels:= table();
        Vertex_index:= 0
    end proc,
    
    ModuleApply:= proc(func::function(function), {Labels::table:= table(), Prefix::nonnegint:= 0})
        ModuleInits();
        thismodule:-Prefix:= Prefix;
        thismodule:-Labels:= Labels;
        
        AddVertex(func);  #Root the tree
        AddSubTree(func); #Phase 1
        ShortenLabels();  #Phase 2
        #Subs the short labels into the integer labels, then build and return the graph
        GraphTheory:-Graph(subs(op(eval(VertexLabels)), convert(Edges, set)))
    end proc       
;
# VertexMap is a table mapping the short labels to their original functions
export VertexMap::table;
end module;

A module to convert a function tree representation (such as is returned by ToInert or by various commands in XMLtools) into a tree represented in the GraphTheory package.

 

(* Written by Carl Love, 16-Mar-2013.

module () local Vertices::table, VertexLabels::table, Vertex_index::nonnegint, Edges::table, Edge_index::nonnegint, AddEdge, AddVertex, AddSubTree, Prefix::nonnegint, Labels::table, StripName, ShortenLabels, ModuleInits, ModuleApply; export VertexMap::table; end module

Func:= ToInert(x+2*y*z^2);

_Inert_SUM(_Inert_PROD(_Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), _Inert_INTPOS(2)), _Inert_NAME("x"))

I came up with a scheme to abbreviate the vertex labels. In the line below, you need to map the second part of the _Inert_ function names to display labels. Note that the display label can be the empty name ``. Also note that this table must map strings (things in double quotes) to names (things in single backquotes).

Abbrs:= table(["SUM"= `+`, "PROD"= `*`, "INTPOS"= ``, "NAME"= ``, "POWER"= `^`]):

The Prefix= 7 argument tells it to strip off the first 7 characters of each function name.

G:= FunctionToTree(Func, Labels= Abbrs, Prefix= 7);

GRAPHLN(directed, unweighted, [`*`, `*`, `+`, `2`, `2`, `^`, x, y, z], Array(%id = 18446744073961315374), `GRAPHLN/table/41`, 0)

For any given tree, it's hard to say which command of DrawNetwork of DrawGraph will display it better. Neither is great.

GraphTheory:-DrawNetwork(G);

GraphTheory:-DrawGraph(G);

The module exports a table VertexMap that shows the correspondence between the vertex labels and the original functional representation.

eval(FunctionToTree:-VertexMap);

table( [( `2` ) = _Inert_INTPOS(2), ( y ) = _Inert_NAME("y"), ( `^` ) = _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2)), ( z ) = _Inert_NAME("z"), ( `*` ) = _Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), ( `*` ) = _Inert_PROD(_Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), _Inert_INTPOS(2)), ( `+` ) = _Inert_SUM(_Inert_PROD(_Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), _Inert_INTPOS(2)), _Inert_NAME("x")), ( `2` ) = _Inert_INTPOS(2), ( x ) = _Inert_NAME("x") ] )

 

Func:= ToInert([1,2,3,4]);

_Inert_LIST(_Inert_EXPSEQ(_Inert_INTPOS(1), _Inert_INTPOS(2), _Inert_INTPOS(3), _Inert_INTPOS(4)))

G2:= FunctionToTree(Func, Labels= Abbrs, Prefix= 7);

GRAPHLN(directed, unweighted, [`1`, `2`, `3`, `4`, EXPSEQ, LIST], Array(%id = 18446744073961304662), `GRAPHLN/table/46`, 0)

GraphTheory:-DrawNetwork(G2);

 


Download FunctionTree.mw

You can see the code for Maple's rifsimp:

showstat(`DEtools/Rif/rifsimp`);

Once you do the first problem, you know that c = 3. That allows you to set the Support option of the Statistics:-Distribution command.

with(Statistics):

MyDist:= Distribution(PDF= (x-> x^2/9), Support= 0..3):

[Mean, Variance, StandardDeviation](MyDist);

[9/4, 27/80, (3/20)*15^(1/2)]

You should also learn the more-traditional ways:

f:= x-> x^2/9:

mean = Int(x*'f'(x), x= 0..3);

mean = Int(x*f(x), x = 0 .. 3)

value(%);

mean = 9/4

mu:= eval(mean,%);

9/4

variance = Int((x-'mu')^2*'f'(x), x= 0..3);

variance = Int((x-mu)^2*f(x), x = 0 .. 3)

value(%);

variance = 27/80

stdev:= sqrt(eval(variance,%));

(3/20)*15^(1/2)

 


Download Stats.mw

The links to your uploaded worksheets are not working. But the idea of the problem is that the integral of f(x) from 0 to c must equal 1.

eqn:= int(x^2/9, x= 0..c) = 1;
solve({eqn, c>0}, c);

See ?intsolve.

One problem with using index numbers with op is that you don't know that the integral will be the, say, 3rd factor in each term. Indeed, the ordering of the factors within each term might change on different runs. If you just want to extract all the integrals from lhs(eq10_1), then do

indets(lhs(eq10_1), specfunc(anything, Int));

If you want to extract the integrals and maintain the order that they have in the equation (and you know that there is more than one term and at most one integral per term), then do

map(op@indets, [op](lhs(eq10_1)), specfunc(anything, Int));

I believe that the effect that you're seeing is just due to the pixel size of your screen, and has nothing to do with the number of points. Curves that are close to horizontal or close to vertical will appear jagged. We have to live with it until they start making screens at 300 dots per inch, or more. Try printing the plot on paper and seeing if the effect is still present.

In three short lines:

R1:= %cos ~ (< 1|1|1 >*omega*t + 2*Pi/3*<0|-1|1>):

T:= sqrt(2/3)*value(< R1, eval(R1, %cos= -%sin), < 1|1|1 >/sqrt(2) >):

(simplify@expand) ~ (T.map(diff, T, t)^%T);

Matrix(3, 3, {(1, 1) = 0, (1, 2) = -omega, (1, 3) = 0, (2, 1) = omega, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0})

 


Download trigmatrix.mw

I think that a properly constructed histogram would be a better-than-the-original substitute for a stem-and-leaf plot. What do think of this (the gridlines don't appear in the worksheet; that's an artifact of MaplePrimes)?

restart;

Data:= [130.8, 129.9, 131.5, 131.2, 129.5, 132.2, 133.7, 127.8,
131.5, 132.7, 134.8, 131.7, 133.9, 129.8, 131.4, 132.8,  132.7,
128.8, 130.6, 131.1, 133.8, 130.5, 131.4, 131.3, 129.5
]:

Statistics:-Histogram(
   Data
  ,binbounds= [$127..135]
  ,frequencyscale= absolute
  ,axis[1]= [tickmarks= [.5 +~ [$127..134] =~ [$127..134]]]
);


Quartiles = map2(Statistics:-Quartile, Data, [1,2,3]);

Quartiles = [HFloat(130.3), HFloat(131.4), HFloat(132.7)]

 


Download stem_and_leaf.mw

You need to separate the parameter definitions from the setting of the initial conditions. A clean way to do this is

eq1 := {
    diff(S_h(t),t)= mu_h*(N_h-S_h(t))-(bi_h*b)/(N_h+m)*S_h(t)*I_v(t),
    diff(I_h(t),t)=(bi_h*b)/(N_h+m)*S_h(t)*I_v(t)-(mu_h+ga)*I_h(t),
    diff(R_h(t),t)=ga*I_h(t)-mu_h*R_h(t),
    diff(S_v(t),t)=A-(bi_v*b)/(N_h+m)*S_h(t)*I_h(t)-mu_v*S_v(t),
    diff(I_v(t),t)=(bi_v*b)/(N_h+m)*S_v(t)*I_h(t)-mu_h*I_h(t)
};
params:= {
    mu_h=0.0000457,mu_v=0.25,b=1,ga=0.167,
    bi_h=0.4, bi_v=0.4,m=0,N_h=10000,A=5000
};
ICs:= {
    S_h(0)=10000, I_h(0)=1,R_h(0)=0,S_v(0)=10000,I_v(0)=1
};

sol1:=dsolve(
    eval(eq1, params) union ICs,
    {S_h(t),I_h(t),R_h(t),S_v(t),I_v(t)},
    type=numeric
);

Go to ?procedure and read about the uses clause of a procedure header. You can't use with inside a procedure, and I'd recommend against writing procedures that rely on a with performed outside the procedure. With the uses clause, you can make abbreviations:

proc()
   uses St= Statistics, LA= LinearAlgebra;
   local A,B, ...;
   ...;
   A:= St:-Rank(...);
   B:= LA:-Rank(...);
   ...
end proc;

Also, my opinion is that A:-B looks more elegant than A[B]. Maple already has enough square brackets; it is nice to balance the variety of symbols.

(A prior version of this answer was messed up due to the well-known bug of the MaplePrimes editor dropping angle brackets.)

The below assumes that the factors are represented as 2x4 Matrices, and it represents the product as a 4x4 Matrix.

`&x`:= proc(AB::Matrix(2,4), CD::Matrix(2,4))
   local C:= CD[.., 1..2], D:= CD[.., 3..4];
   < < AB[1,1]*C + AB[1,3]*D | AB[1,2]*C + AB[1,4]*D >,
     < AB[2,1]*C + AB[2,3]*D | AB[2,2]*C + AB[2,4]*D >
   >
end proc:           

Example of use:

A,B,C,E:= seq(LinearAlgebra:-RandomMatrix(2, 2, generator= rand(-1..1)), k= 1..4);

A, B, C, E := Matrix(2, 2, {(1, 1) = -1, (1, 2) = -1, (2, 1) = 0, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = -1, (1, 2) = 1, (2, 1) = 1, (2, 2) = -1}), Matrix(2, 2, {(1, 1) = -1, (1, 2) = -1, (2, 1) = 0, (2, 2) = -1}), Matrix(2, 2, {(1, 1) = -1, (1, 2) = 1, (2, 1) = -1, (2, 2) = -1})

'`&x`'('`<|>`'(A,B), '`<|>`'(C,E)) = < A|B > &x < C|E >;

`&x`(`<|>`(Matrix(2, 2, {(1, 1) = -1, (1, 2) = -1, (2, 1) = 0, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = -1, (1, 2) = 1, (2, 1) = 1, (2, 2) = -1})), `<|>`(Matrix(2, 2, {(1, 1) = -1, (1, 2) = -1, (2, 1) = 0, (2, 2) = -1}), Matrix(2, 2, {(1, 1) = -1, (1, 2) = 1, (2, 1) = -1, (2, 2) = -1}))) = (Matrix(4, 4, {(1, 1) = 2, (1, 2) = 0, (1, 3) = 0, (1, 4) = 2, (2, 1) = 1, (2, 2) = 2, (2, 3) = -1, (2, 4) = 0, (3, 1) = -1, (3, 2) = 1, (3, 3) = 1, (3, 4) = -1, (4, 1) = -1, (4, 2) = -1, (4, 3) = 1, (4, 4) = 1}))

NULL


Download tensor.mw

 

Check out ?ToInert and ?dismantle. The former returns a result that is logically identical to what you want; it is just not presented visually as a tree. The latter displays the output in tree form, but includes details that you might find irrelevant. It would be relatively easy to convert the ToInert output to an explicitly presented tree by using the GraphTheory package and its DrawNetwork command.

Example:

ToInert(x+2*y*z^2);

 

First 381 382 383 384 385 386 387 Last Page 383 of 395