Carl Love

Carl Love

26386 Reputation

25 Badges

11 years, 142 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

You have three eval statements of the form

eval(f, x= 0 and diff(x,x)= 0);

Those should be of the form

eval(f, [x= 0, diff(f(x), x)= 0]);

Note that in addition to getting rid of the and, I changed diff(x,x) to diff(f(x), x). You may have intended that to be diff(f1(x), x). You also may have intended for the initial f to be ODE, as in eval(ODE, [x= ...]).

The basic technique is

  1. Assign fsolve's result to a variable:
     Sol:= fsolve(...);
  2. Use eval to evaluate an expression with respect to the solution:
     eval(X[CO2], Sol);

You assigned fsolve's solution to l, but I don't like that because it is difficult to distiguish from 1 or I in many fonts. So I used Sol instead. Make the last three lines of your for loop this:

     Sol:= fsolve(sys3,{a,b,c,d,e,k1,k2,k3}):
     print('T'= T, map(c-> 'X[c]'= eval(X[c], Sol), [CH4, CO2, H2O, H2, CO]))
end do:

and you will get one line of output for each loop iteration, each line like this (but it will all fit on one line):

I know that that is somewhat cryptic, so please let me know if it works for you and if you understand.

From your questions of a few days ago, I know that you are getting to this point from a group of matrices. In particular, you think that there is some correspondence between

1  1  2
1  2  2

and the matrix

1  1
0  1.

Here, I will generate several representations of a small group of matrices, the 2x2 matrices with determinant 1 over the finite field Z2. The above matrix is in this group. Since I think that you don't have Maple 17, I won't use the new GroupTheory package.

restart:

Generate all 2x2 matrices of determinant 1 in Z2.

G:= table():

for i to 2 do for j to 2 do for k to 2 do for m to 2 do
     A:= < < i, j > | < k, m > > mod 2;
     if Det(A) mod 2 = 1 then G[A]:= [][] fi
od od od od:

G:= [{indices}(G, nolist)[]];

G := [Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (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) = 0, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0})]

Convert matrices to listlist form because member (like other normal equality tests) doesn't work on matrices.

G1:= map(convert, G, listlist);

[[[1, 1], [1, 0]], [[1, 0], [1, 1]], [[1, 1], [0, 1]], [[1, 0], [0, 1]], [[0, 1], [1, 1]], [[0, 1], [1, 0]]]

This is a multplication operator for G where each element is represented by its index in the list.

Gmul:= proc(i,j)
     local z,p;
     z:= convert(G[i].G[j] mod 2, listlist);
     member(z,G1,'p');
     p
end proc:   

Use that to generate a multiplication table (also called a Cayley table) for the group.

C:= Matrix(6, 6, Gmul);

C := Matrix(6, 6, {(1, 1) = 5, (1, 2) = 6, (1, 3) = 2, (1, 4) = 1, (1, 5) = 4, (1, 6) = 3, (2, 1) = 3, (2, 2) = 4, (2, 3) = 1, (2, 4) = 2, (2, 5) = 6, (2, 6) = 5, (3, 1) = 6, (3, 2) = 5, (3, 3) = 4, (3, 4) = 3, (3, 5) = 2, (3, 6) = 1, (4, 1) = 1, (4, 2) = 2, (4, 3) = 3, (4, 4) = 4, (4, 5) = 5, (4, 6) = 6, (5, 1) = 4, (5, 2) = 3, (5, 3) = 6, (5, 4) = 5, (5, 5) = 1, (5, 6) = 2, (6, 1) = 2, (6, 2) = 1, (6, 3) = 5, (6, 4) = 6, (6, 5) = 3, (6, 6) = 4})

So, C[i,j] is the index of  G[i].G[j].

 

Each row (and each column, for that matter) of the table is a permutation of the six elements of the group.

 

The first row is a permutation of order 3.

x1:= convert(convert(C[1,..], list), disjcyc);

[[1, 5, 4], [2, 6, 3]]

That corresponds to G[1] having order 3.

G[1]^2 mod 2, G[1]^3 mod 2;

Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1})

The second row has order 2.

y1:= convert(convert(C[2,..], list), disjcyc);

[[1, 3], [2, 4], [5, 6]]

G[2]^2 mod 2;

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1})

An element of order 3 and an element of order 2 must generate a group of order at least lcm(2,3) = 6, so in this case they must generate the whole group.

G2:= permgroup(6, {x= x1, y= y1}):

group:-grouporder(G2);

6

Now let's look at "words" constructed from the group generators.

convert([x], disjcyc, G2);

[[1, 5, 4], [2, 6, 3]]

convert([x,x], disjcyc, G2);

[[1, 4, 5], [2, 3, 6]]

convert([x,y], disjcyc, G2);

[[1, 6], [2, 5], [3, 4]]

convert([x,x,y], disjcyc, G2);

[[1, 2], [3, 5], [4, 6]]

convert([y,x], disjcyc, G2);

[[1, 2], [3, 5], [4, 6]]

So xxy = yx. Thus xxyx^(-1)y^(-1) = e, the identity element. This gives us enough information to represent this as a "finitely presented" group.

G3:= grelgroup({X,Y}, {[X,X,X], [Y,Y], [X,X,Y,1/X,1/Y]}):

Represent the trivial identity subgroup so that we can get its cosets.

E:= subgrel({e= []}, G3):

Thus the following is another permutation representation of the same group.

G4:= group:-permrep(E, G3);

permgroup(6, {X = [[1, 3, 4], [2, 5, 6]], Y = [[1, 2], [3, 6], [4, 5]]})

 

 

Download groups.mw

You can merge an arbitrary sequence of m column Vectors into a Matrix like this:

A:= `<|>`(seq(`&varphi;`[1,k], k= 1..m));

If making an assumption about g is appropriate for your situation, then use evalc (which assumes that g is real) and assume g > -4.

ex:= -(1/2*I)*(I*Pi*ln(2) - I*Pi*ln(g+4) - 2*ln(2)^2
      + 2*ln(2)*ln(g+4) - 2*ln(2)*ln(-1/2*I) + ln(g+4)*ln(-1/2*I)
      - ln(g+4)*ln(2*I) + ln(-1/2*I)*ln(4-g) + ln(4-g)*ln(2*I))
     /Pi:
evalc(ex) assuming g>-4;
                    Pi ln(2) - Pi ln(g + 4)
                    -----------------------
                              Pi           
combine(%);
                             /1      \
                          -ln|- g + 2|
                             \2      /

Surprisingly, Maple can get the symbolic solution to this very quickly if you give it just a tiny bit of help. It is trivial to solve ode1 for V(t) and substitute that into ode2. Then dsolve solves it in a few seconds. 

restart:

ode1:= diff(L(t),t)*(V(t))=0.682*10^(-7);

(diff(L(t), t))*V(t) = 0.682000000000000e-7

SolV:= solve(ode1, {V(t)})[];

V(t) = 0.682000000000000e-7/(diff(L(t), t))

ode2:= 137*10^6*diff(V(t),t)*V(t)=(1/(4*3.14))*1.66*10^(-10)*L(t)^2*(3.5*10^3-V(t))^2;

137000000*(diff(V(t), t))*V(t) = 0.132165605095541e-10*L(t)^2*(3500.0-V(t))^2

ode:= eval(ode2, SolV);

-0.637219880000000e-6*(diff(diff(L(t), t), t))/(diff(L(t), t))^3 = 0.132165605095541e-10*L(t)^2*(3500.0-0.682000000000000e-7/(diff(L(t), t)))^2

ode:= solve(ode, {diff(L(t),t$2)})[];

diff(diff(L(t), t), t) = -254.076922775915*(diff(L(t), t))^3*L(t)^2+0.990174064760994e-8*(diff(L(t), t))^2*L(t)^2-0.964712445952854e-19*L(t)^2*(diff(L(t), t))

ic1:=V(0)=0.01;

V(0) = 0.1e-1

ic2:=L(0)=0;

L(0) = 0

ic3:= subs(V(t)= eval(V(0), ic1), ode1);

0.1e-1*(diff(L(t), t)) = 0.682000000000000e-7

ic3:= D(L)(0) = solve(ic3, diff(L(t), t));

(D(L))(0) = 0.682000000000000e-5

Sol:= dsolve({ode, ic3, ic2}, {L(t),V(t)}, convert_to_exact= false);

L(t) = RootOf(-(Int(-508.153845551830/(0.316227766016838e-15*tan(0.527046276694730e-16*_a^3+arctan(-10959187324298.1))-0.990174064760994e-8), _a = 0 .. _Z))+t)

SolL:= eval(Sol);

L(t) = RootOf(-(Int(-508.153845551830/(0.316227766016838e-15*tan(0.527046276694730e-16*_a^3-1.57079632679481)-0.990174064760994e-8), _a = 0 .. _Z))+t)

SolV:= eval(SolV, SolL);

V(t) = 0.682000000000000e-7/(-0.622307139432214e-18*tan(0.527046276694730e-16*RootOf(-(Int(-508.153845551830/(0.316227766016838e-15*tan(0.527046276694730e-16*_a^3-1.57079632679481)-0.990174064760994e-8), _a = 0 .. _Z))+t)^3-1.57079632679481)+0.194857142857143e-10)

Sol:= {SolV,SolL}:

plot(eval(L(t), Sol), t= 0..100);

plot(eval(V(t), Sol), t= 0..100);

 



Download moji_solved.mw

 

 

 

This is not a problem with LinearAlgebra:-Equal.

When a Matrix (or other rtable) contains a symbolic (the a in your case) and that symbolic is changed, the Matrix is not automatically updated. One can force the update with rtable_eval:

A := Matrix([[a, 0], [0, 0]]);
for a from 0 to 10 do
     A1:= rtable_eval(A);
     if LinearAlgebra:-Equal(A1^2, A1) then print(a, A1) end if
end do;

I consider it a bad programming practice to use a for loop's index variable as a symbolic. However, Maple does allow it. So I'd rewrite the code as

A:= < < a, 0 > | < 0, 0 > >:
for k from 0 to 10 do
     A1:= eval(A, a= k);
     if LinearAlgebra:-Equal(A1^2, A1) then print(k, A1) end if
end do: 

Why the method labelled "PROCEDURE" isn't working for you:

(I think that you already figured this out on your own.)

In this method, you've just moved the global/parameter conflict to an earlier procedure. In your procedure named qpercentagedifference, the parameter height is not the same as the global height that is "in" the variables qtotal and qpresent. Note that expressions created with the arrow (->) and with proc are essentially the same to Maple--they are both considered procedures.

In order to use this method, you have to be consistent "down through the chain" and make every expression using height a procedure with height as a parameter.

Why the method labelled "SUBS" isn't working for you:

The subs does not change the expression being substituted; rather it creates a new expression. For example, subs(x= _x, y) does not change y. There are several things that you can do to fix this.

(1)
modalmasspedestriansfail:= proc(subsheight)
     if subs(height= subsheight, qpercentagedifference) >= 5/100 then
          subs(height= subsheight, modalmasswithpedestrians)
     else
          subs(height= subsheight, modalmasswithoutpedestrians)
     end if
end proc:

(2)modalmasspedestriansfail:= proc(subsheight)
     local A, B, C;
     A,B,C:= subs(
          height= subsheight,
          [qpercentagedifference,
           modalmasswithpedestrians,
           modalmasswithoutpedestrians
          ]

    )[];
     if A >= 5/100 then B else C end if
end proc

(3)
modalmasspedestriansfail:= proc(subsheight)
     subs(
          height= subsheight,
            `if`(
                  qpercentagedifference >= 5/100
,
              modalmasswithpedestrians,
                  modalmasswithoutpedestrians
             )
     )
end proc:

 

The problem has nothing to do with the if...then. In both of your examples, there are two different variables named height, one a global and one a procedure parameter. In your first example--the one the doesn't work--you are trying to use the global as the parameter. Here's a simplified version of your first example:

y:= x^2:
f:= proc(x) y end proc:

The x that is "in" y is not the same as the parameter x. This procedure won't plot.

There are several ways around this problem. One way is to make y a procedure:

y:= x-> x^2:
f:= proc(x) y(x) end proc:

A second way is to use a different name and subs one for the other:

y:= x^2:
f:= proc(_x) subs(x= _x, y)
end proc:

A third way is unapply:

y:= x^2:
Y:= unapply(y, x):
f:= proc(x) Y(x) end proc:

And there are several more ways.

In your second example, the one that works, the height appears explicitly in the procedure. But, as you probably realize, this is not a good way to do it because you had to copy the expression into the procedure. Copying sections of code is always undesirable.

It looks like Maple ran out of memory trying to solve your equation symbolically. Get a numeric solution instead.

 

Sol:= dsolve({ode1,ode2,ic1,ic2},{V(t),L(t)}, numeric);

You can plot the solutions with plots:-odeplot. See ?dsolve,numeric and ?plots,odeplot .

Example: Plot t vs. L(t):

plots:-odeplot(Sol, [[t,L(t)]], t= 0..100)

For some reason I can't compile on my computer. (And I'd appeciate some advice on how to correct that.)

Anyway, Acer's evalhf code runs in 2.55 secs on my computer. By parallelizing that and making some small tweaks to the innermost loop (removing the extra trunc), I got that down to 0.815 secs (modulo my eight processors). The parallelization should also work with the compiled version, but I can't test that. (But I left in the code for the compiled version, of course.) Here's the code: Simu_para_evalhf.mw

Like this:

BVP:= {diff(y(x),x$2)-2*y(x)=0,y(0)=1.2,y(1)=0.9}:
Sol_num:= dsolve(BVP, numeric, method=bvp[midrich], output= listprocedure):
Sol_exact:= dsolve(BVP):
Res:= eval(y(x), Sol_num) - unapply(eval(y(x), Sol_exact), x):
plot(Res, 0..1);

The remaining time can be cut in half (!!!) simply by replacing frac(x) with x - trunc(x) since trunc(x) has already been computed. frac is symbolic, and trunc is kernel. It's amazing how much of a difference this makes.

Like most code that generates a lot of random numbers, this one is an obvious candidate for parallelization. It just required a trivial change to your outermost procedure and converting writable globals to locals (which is good coding practice anyway). The benefit that you derive from this will be proportional to the number of processors that you have. I have eight, most modern home computers have four. Windows Resource Monitor showed that I had 100% utilization on all processors, so parallelization can't get any better than that.

Making these two changes, I got the time down to 6 seconds. Here is the modified code: Simulation_Para.mw

There should be no loop for i from 1 to Size in procedure SampleQ. The loop in Qsim is already for i from 1 to Size. So SampleQ should begin

SampleQ:=proc()
    global Q,obs,forv,A,rSum,cSum;
    local number,i,j,k,n,m,broekdel;    
   
     obs:=Matrix(r,c):   #Creating a r times c matrix with zeroes - to hold the observed values
     forv:=Matrix(r,c):   #Creating a r times c matrix with zeroes - to hold the expected values
     A:=Sample(X,1..Total):   

...

This reduces the run time to 51 seconds for me.

Since I think that you're using a less-than-most-recent version of Maple, see ?grelgroup . If you've upgraded to Maple 17 since we've last talked about it, then see ?GroupTheory,Group . In either case, "words" are represented as lists; so instead of a*b*a^(-1), you'd use [a,b,a^(-1)] or [a,b,1/a].

To represent the group that you used in your example:

G:= grelgroup({a,b}, {[a,a], [b,b], [a,b,1/a,1/b]});

(You can use a^(-1) instead of 1/a if you prefer.)

In this representation, a and b are called the generators.

G is a group, but not a permutation group because its elements are not permutations (a is just a letter, obviously; so it's not a permutation). However, every group is isomorphic to some permutation group. G has four elements: [e, a, b, ab] (where e is the identity element). Make that list correspond to [1,2,3,4]. Now we wish to find the permutation that corresponds to a. Multiply each element of the list by a (this group is abelian (commutative), so order doesn't matter) and you get [ae, aa, ab, aab] = [a, e, ab, b], which corresponds to [2,1,4,3], which is [[1,2],[3,4]] in disjcyc notation. So, in a sense of "is", a "is" the permutation [[1,2],[3,4]], but only because of the order that we listed the elements. Likewise, if we multiply the original list by b, we find that b corresponds to the permutation [[1,3], [2,4]]. Those two elements are enough to generate the group. We can now say

GP:= permgroup(4, {a= [[1,2],[3,4]], b=[[1,3],[2,4]]});

Now, it turns out that this is the minimal-degree permutation representation of this particular group, but that was just luck. The technique outlined here will always give a representation of degree the same as the order of the group, but there is usually a representation of lesser degree.

I hope that that increased your understanding and didn't cause more confusion.

First 347 348 349 350 351 352 353 Last Page 349 of 381