Joe Riel

9660 Reputation

23 Badges

20 years, 3 days

MaplePrimes Activity


These are answers submitted by Joe Riel

Events generally occur when conditionals in the model change.  For example, a step is modeled as

y = if time < 1 then 0 else 1;

At time = 1 an event occurs. This snippet is written in Modelica, which MapleSim uses as the model language.  MapleSim parses the Modelica language, generates Maple code with events that can be handled by dsolve (it's significantly more complicated than that, but that's the basic idea).

 

Use simplify with side-relations, see ?simplify[siderels]

p := x^10*y^2+x^9*y^8+x^7*y^3+x^5*y+x*y+y:

eqs := {NULL
        , x^3=a*x+y
        , x^2*y=b*x+c*y
        , x*y^2=x+d*y
        , y^3=e*y
       }:

simplify(p, eqs);
 4    3      4    2      3  4        2  2  3      4          4    2    4
a  d e  y + a  d e  y + a  e  y + 3 a  d  e  y + a  b x y + a  c y  + a  d e y

          2  2  2            4      4          3    2      2  2
     + 6 a  d  e  y + 3 a d e  y + a  d y + 5 a  d y  + 9 a  d  e y

              3      5      4        3           2             2
     + 3 a d e  y + e  y + a  x + 5 a  x y + 12 a  b d y + 12 a  c e y

          2    2            2            2          2        2          2
     + 4 a  e y  + 5 a b d y  + 5 a c e y  + 3 a d e  y + 3 a  b x - 9 a  c y

        2                                      2                2
     + a  x y + 5 a b x y + 3 a d e y + 2 d e y  + 3 a d y + a y  + b x y

          2        2
     + c y  + 2 d y  + 3 a x + 3 x y + y


You can do

A := ImportMatrix("A.txt"):
B := ImportMatrix("B.txt"):
C := ImportMatrix("C.txt"):

with(DynamicSystems):

ss := StateSpace(A,B,C):

BodePlot(ss, range = 0.1..100);

The system has 2 inputs and 14 outputs, so you get 28 superimposed plots.  You can plot just one by doing, say

BodePlot(ss, range = 0..1000, subsystem=[1,2]);

Is the source for the procedures text files or a Maple worksheet?  If the former, just use a text editor to search/replace Return with return. For a Maple worksheet you could do the same provided you used Maple input format rather than 2D format as the input region. That is, edit the *.mw file in text editor and change Return(whatever) to return whatever. If you used 2D format, select everything in the worksheet (Ctrl-A), right click, pick 2D Math -> convert -> 1D Math input.  Save the worksheet, open in an editor, change as described above.  If you prefer the 2D input, then reopen the worksheet, select all, and convert all input back to 2D input.

If you have to convert between 2D and 1D math, I'd be sure to save a backup of the worksheet before changing anything.

 

 

pde := { diff(u(x,t),t) = k*diff(u(x,t),x,x) };
bc := { u(0,t)=0, u(l,t)=0, u(x,0)=exp(-x) };

pdsolve(pde union bc) assuming l>0;
          infinity
           -----
            \
u(x, t) =    )
            /
           -----
          _Z2~ = 1

    /                                        2     2        3                \
    |                           _Z2~      -Pi  _Z2~  k t - l       Pi _Z2~ x |
    |  2 Pi _Z2~ (-exp(l) + (-1)    ) exp(-------------------) sin(---------)|
    |                                              2                   l     |
    |                                             l                          |
    |- ----------------------------------------------------------------------|
    |                                2     2    2                            |
    \                              Pi  _Z2~  + l                             /


I'm not clear on what you are doing (nor how to interpret zero.ta[2]) but you could try the following

cfs := {coeffs}(p,x):
eqs := map(`=`, cfs, zero_ta[2]*beta[3]*beta[4]*omega^2);
sol := [solve(eqs)];

That gives six different solutions, however, I don't know what are parameters and what are constants.

This algorithm returns duplicates and also misses some paths.  To see that execute Routes(4,2) and inspect the first two results (T[1] and T[2]).  With duplicates removed, it returns only 82 paths for Route(4,1), the correct answer is easily computed to be 3*4 + 5*8 + 8*4 = 84.

I believe the correct result, for Routes(4,6), is 68272. I computed that using a simpler approach, an extension to the Iterator package I've recently described.

Addendum

The original code for Routes give the correct answer for Routes(4,1), 84.  It is my modification of it that is flawed for that case. However, the original is incorrect with Routes(4,2), it has the duplicates I've mentioned. 

Further Add

Actually, while the original code does return 84 for Routes(4,1), it returns only 82 unique routes.

Correction

The bug is in the points used for the left side; there is a duplicate.  With that fixed I get the same results as with the Iterator method. Here is the revised procedure

Routes:=proc(N::posint, n::nonnegint)

global T;
local i,j,L, Rule;

    if n>=N^2 then
        T:=[]:
        0;
    else
        Rule := proc(K)  # Continuation of the route by 1 step
        local S, k, r, rk, p, pts, j;

            S := table();
            j := 0;
            k := nops(K[1]);

            for r in K do
                # Assign pts the points to consider
                if r[k]=[1, 1] then
                    # Bottom left corner
                    pts := [[1, 2], [2, 2], [2, 1]];
                elif r[k]=[1, N] then
                    # Top left corner
                    pts := [[1, N-1], [2, N-1], [2, N]];
                elif r[k]=[N, N] then
                    # Top right corner
                    pts := [[N-1, N], [N-1, N-1], [N, N-1]];
                elif r[k]=[N, 1] then
                    # Bottom right corner
                    pts := [[N-1, 1], [N-1, 2], [N, 2]];
                elif r[k,1]=1 and r[k,2]<>1 and r[k,2]<>N then
                    # Left side
                    pts := [[1, r[k,2]-1], [2, r[k,2]-1], [2, r[k,2]], [2, r[k,2]+1], [1, r[k,2]+1]];
                elif r[k,2]=N and r[k,1]<>1 and r[k,1]<>N then
                    # Top side
                    pts := [[r[k,1]-1, N], [r[k,1]-1, N-1], [r[k,1], N-1], [r[k,1]+1, N-1], [r[k,1]+1, N]];
                elif r[k,1]=N and r[k,2]<>1 and r[k,2]<>N then
                    # Right side
                    pts := [[N, r[k,2]+1], [N-1, r[k,2]+1], [N-1, r[k,2]], [N-1, r[k,2]-1], [N, r[k,2]-1]];
                elif r[k,2]=1 and r[k,1]<>1 and r[k,1]<>N then
                    # Bottom side
                    pts := [[r[k,1]-1, 1], [r[k,1]-1, 2], [r[k,1], 2], [r[k,1]+1, 2], [r[k,1]+1, 1]];
                elif r[k,1]<>1 and r[k,1]<>N and r[k,2]<>1 and r[k,2]<>N then
                    # Inside
                    pts := [[r[k,1]-1, r[k,2]-1], [r[k,1]-1, r[k,2]], [r[k,1]-1, r[k,2]+1], [r[k,1], r[k,2]+1], [r[k,1]+1, r[k,2]+1], [r[k,1]+1, r[k,2]], [r[k,1]+1, r[k,2]-1], [r[k,1], r[k,2]-1]];
                fi;
                rk := r[1..k-1];
                j := j+1;
                S[j] := seq(`if`(member(p,rk)
                                 , NULL
                                 , [op(r),p]
                                ), p = pts);
            od;
            [seq(S[j], j=1..j)];
        end proc;

        L:=[seq(seq([[i, j]], j=1..N), i=1..N)];
        T:=(Rule@@n)(L);
        nops(T);  # Number of all the routes
    fi;
end proc:

As Carl points out, this is a weird kernel/gui bug.  To work around it, terminate the call to dsolve with a colon (so nothing is printed). You can then use soltn as desired:

(**) eqn1 := diff(R1(r),r) + I*R2(r) = 0:
(**) eqn2 := diff(R2(r),r) + R1(r) = 0:
(**) st := 1.3:
(**) fn := 4.3:
(**) soltn := dsolve([eqn1,eqn2,R1(st)=0.1,R2(st)=0.2], [R1(r),R2(r)], numeric, range=st..fn): # terminate with a colon
(**) soltn(2);
   [r = 2., R1(r) = 0.110429803108025 - 0.115236213462622 I, R2(r) = 0.128139507027451 + 0.0432522993325482 I]

Correction:

My work-around doesn't really do anything other than prevent the message from being displayed. Originally I hadn't realized that soltn was being assigned correctly, regardless. I've submitted this as an SCR. The bug lies in the display of soltn. To see that, execute

print(soltn);
                       Non-fatal error while reading data from kernel.

By avoiding n^2 list building and using a few other refinements, I reduced the run time of Routes from about 20 seconds to 1/2 second. Here is the modified code

Routes:=proc(N::posint, n::nonnegint)

global T;
local i,j,L, Rule;

    if n>=N^2 then
        T:=[]:
        0;
    else
        Rule := proc(K)  # Continuation of the route by 1 step
        local S, k, r, rk, p, pts, j;

            S := table();
            j := 0;
            k := nops(K[1]);

            for r in K do
                # Assign pts the points to consider
                if r[k]=[1, 1] then
                    # Bottom left corner
                    pts := [[1, 2], [2, 2], [2, 1]];
                elif r[k]=[1, N] then
                    # Top left corner
                    pts := [[1, N-1], [2, N-1], [2, N]];
                elif r[k]=[N, N] then
                    # Top right corner
                    pts := [[N-1, N], [N-1, N-1], [N, N-1]];
                elif r[k]=[N, 1] then
                    # Bottom right corner
                    pts := [[N-1, 1], [N-1, 2], [N, 2]];
                elif r[k,1]=1 and r[k,2]<>1 and r[k,2]<>N then
                    # Left side
                    pts := [[1, r[k,2]-1], [2, r[k,2]-1], [2, r[k,2]-1], [2, r[k,2]+1], [1, r[k,2]+1]];
                elif r[k,2]=N and r[k,1]<>1 and r[k,1]<>N then
                    # Top side
                    pts := [[r[k,1]-1, N], [r[k,1]-1, N-1], [r[k,1], N-1], [r[k,1]+1, N-1], [r[k,1]+1, N]];
                elif r[k,1]=N and r[k,2]<>1 and r[k,2]<>N then
                    # Right side
                    pts := [[N, r[k,2]+1], [N-1, r[k,2]+1], [N-1, r[k,2]], [N-1, r[k,2]-1], [N, r[k,2]-1]];
                elif r[k,2]=1 and r[k,1]<>1 and r[k,1]<>N then
                    # Bottom side
                    pts := [[r[k,1]-1, 1], [r[k,1]-1, 2], [r[k,1], 2], [r[k,1]+1, 2], [r[k,1]+1, 1]];
                elif r[k,1]<>1 and r[k,1]<>N and r[k,2]<>1 and r[k,2]<>N then
                    # Inside
                    pts := [[r[k,1]-1, r[k,2]-1], [r[k,1]-1, r[k,2]], [r[k,1]-1, r[k,2]+1], [r[k,1], r[k,2]+1], [r[k,1]+1, r[k,2]+1], [r[k,1]+1, r[k,2]], [r[k,1]+1, r[k,2]-1], [r[k,1], r[k,2]-1]];
                fi;
                rk := r[1..k-1];
                j := j+1;
                S[j] := seq(`if`(member(p,rk)
                                 , NULL
                                 , [op(r),p]
                                ), p = pts);
            od;
            [seq(S[j], j=1..j)];
        end proc;

        L:=[seq(seq([[i, j]], j=1..N), i=1..N)];
        T:=(Rule@@n)(L);  # List of all the routes
        nops(T);  # Number of all the routes
    fi;
end proc:
m := 5:
LinearAlgebra:-BandMatrix([[seq(a(k),k=2..m-1)],[seq(b(k),k=1..m-1)],[seq(a(k),k=1..m-1)]]);
                                      [b(1)    a(1)     0       0       0  ]
                                      [                                    ]
                                      [a(2)    b(2)    a(2)     0       0  ]
                                      [                                    ]
                                      [ 0      a(3)    b(3)    a(3)     0  ]
                                      [                                    ]
                                      [ 0       0      a(4)    b(4)    a(4)]


with(Statistics):
X := RandomVariable(Binomial(10,1/2));
ProbabilityFunction(X,3);
                                      15
                                      ---
                                      128


The problem lies in the first ln: (1+(XL)^(2))(1+(YL)^(2).  There should be a multiplicative operator between the two additive terms. Without that, Maple interprets this as function application.  I'll make a change to AlgEquation so that it returns the offending terms when it detects an error.

Edit -> Remove Output (pick either Worksheet or Selected Region)

You can use the Iterator package to solve this puzzle with nested parentheses.  To step through all possible groupings,  use the BinaryTrees iterator.  The following appliable module does that. Understanding precisely how it does it might take a bit of study. This isn't ideal; one really should prune the trees to prevent, for example, both the following branches (1 + 2) + 3 and 1 + (2 + 3) from appearing.  Doing so isn't hard but I'll leave that as an exercise.  Looping through all possibilities takes a while; here I illustrate the operation using n=6, which limits the tree to six internal nodes (7 leaves), so the digits are from 1 to 7 rather than 1 to 9.

DigitalPuzzle := module()
export ModuleApply;
local Expr, Format;

    # Given the L and R Arrays from the BinaryTree iterator,
    # construct an expression corresponding to the tree,
    # where o[v[i]] is the arithmetic operator of the i-th
    # internal node (the v-Array is used to change the operators
    # for a give tree).
    Expr := proc(L,R)
    local leaf,prefix;
    global o,v;
        leaf := 0;
        prefix := proc(i)
            if i=0 then leaf := leaf+1;
            else 'o'['v'[i]](prefix(L[i]),prefix(R[i]));
            end if;
        end proc;
        prefix(1);
    end proc:

    # Given the L and R arrays of the current tree, and an array, v,
    # that maps the internal nodes to the arithmetic operators (stored
    # in o), return a string corresponding to the desired equality.
    # This is only used for formatting a correct result; speed is
    # not a significant issue.
    Format := proc(L,R,v,o,targ)
    local leaf,infix;
        leaf := 0;
        infix := proc(i)
            if i=0 then leaf := leaf+1;
            else sprintf("(%A %A %A)", infix(L[i]), o[v[i]], infix(R[i]));
            end if;
        end proc;
        sprintf("%s = %a", infix(1), targ);
    end proc:

    ModuleApply := proc(n::posint, targ:=100)
    local A,Accept,cnt,ops,Op,BT,LR,s;
    uses Iterator;

        # Array of arithmetic operators
        ops := Array(0..3,[`+`,`*`,`-`,`/`]);

        # Template for the accept predicate.  The try/catch is
        # needed to handle division by zero.
        Accept := proc(v,o:=ops)
            try
                evalb(_ex=targ);
            catch:
                false;
            end try;
        end proc;

        A := ()->NULL; # dummy procedure that is replaced

        # Construct an iterator that generates all possible values of
        # arithmetic operators (as an Array with values from 0 to 3
        # corresponding to the four operations), but accepts (outputs)
        # only those that meet the criteria.
        Op := MixedRadixTuples([4$n], 'accept'=A);

        # Construct an iterator that generates all binary trees
        # with n-internal nodes (and n+1 leaves).
        BT := BinaryTrees(n);

        cnt := 0; # success counter

        # Loop through all binary trees
        for LR in BT do
            # Reassign A, which is the accept predicate in
            # the Op[erator] iterator, specializing it to the
            # selected tree.
            A := subs(_ex=Expr(LR),op(Accept));

            # Loop through all possibilities of assigning
            # the arithmentic operators to the internal nodes of the
            # tree, keeping only those that evaluate to the target.
            for s in Op do
                cnt := cnt+1;
                printf("%s\n", Format(LR,s,ops,targ));
            end do;
            reset(Op);
        end do;
        cnt;
    end proc:
end module:

CodeTools:-Usage(DigitalPuzzle(6));
(1 + (((((2 + 3) * 4) * 5) + 6) - 7)) = 100
((1 + ((((2 + 3) * 4) * 5) + 6)) - 7) = 100
(((1 + (((2 + 3) * 4) * 5)) + 6) - 7) = 100
(((1 - 2) + 3) + (((4 * 5) - 6) * 7)) = 100
(1 - ((2 * 3) - (((4 + 5) + 6) * 7))) = 100
(1 - ((2 - 3) - (((4 * 5) - 6) * 7))) = 100
((1 - (2 * 3)) + (((4 + 5) + 6) * 7)) = 100
((1 - (2 - 3)) + (((4 * 5) - 6) * 7)) = 100
((1 - 2) + (3 + (((4 * 5) - 6) * 7))) = 100
((1 * 2) * (3 + (((4 + 5) * 6) - 7))) = 100
(1 * (2 * (3 + (((4 + 5) * 6) - 7)))) = 100
(1 - (2 - (3 + (((4 * 5) - 6) * 7)))) = 100
((1 * 2) * ((3 + ((4 + 5) * 6)) - 7)) = 100
(1 * (2 * ((3 + ((4 + 5) * 6)) - 7))) = 100
((1 + ((2 * 3) * 4)) * ((5 + 6) - 7)) = 100
((1 + (2 * (3 * 4))) * ((5 + 6) - 7)) = 100
(1 - ((2 * 3) - ((4 + (5 + 6)) * 7))) = 100
((1 - (2 * 3)) + ((4 + (5 + 6)) * 7)) = 100
(1 + ((((2 + 3) * (4 * 5)) + 6) - 7)) = 100
((1 + (((2 + 3) * (4 * 5)) + 6)) - 7) = 100
(1 + ((((2 + 3) * 4) * 5) + (6 - 7))) = 100
((1 + (((2 + 3) * 4) * 5)) + (6 - 7)) = 100
(((1 + ((2 + 3) * (4 * 5))) + 6) - 7) = 100
((1 - (2 * 3)) * ((4 * 5) * (6 - 7))) = 100
((1 - (2 * 3)) * ((4 * 5) / (6 - 7))) = 100
((((1 - (2 * 3)) * 4) * 5) * (6 - 7)) = 100
((((1 - (2 * 3)) * 4) * 5) / (6 - 7)) = 100
(1 + (((2 + 3) * (4 * 5)) + (6 - 7))) = 100
((1 + ((2 + 3) * (4 * 5))) + (6 - 7)) = 100
((1 + ((2 * 3) * 4)) * (5 + (6 - 7))) = 100
(((1 - (2 * 3)) * (4 * 5)) * (6 - 7)) = 100
(((1 - (2 * 3)) * (4 * 5)) / (6 - 7)) = 100
(((1 - (2 * 3)) * 4) * (5 * (6 - 7))) = 100
(((1 - (2 * 3)) * 4) * (5 / (6 - 7))) = 100
((1 + (2 * (3 * 4))) * (5 + (6 - 7))) = 100
((1 - (2 * 3)) * (4 * (5 * (6 - 7)))) = 100
((1 - (2 * 3)) * (4 * (5 / (6 - 7)))) = 100
memory used=370.72MiB, alloc change=24.00MiB, cpu time=6.43s, real time=6.60s
                                                                    37





@JotaTR You should contact Maple Support. Did you install the MapleSim 5.02 update?  That will probably fix this.

First 34 35 36 37 38 39 40 Last Page 36 of 114