dharr

Dr. David Harrington

8482 Reputation

22 Badges

21 years, 33 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

seq(NumberTheory:-CalkinWilfSequence(i),i=1..10)

1, 1/2, 2, 1/3, 3/2, 2/3, 3, 1/4, 4/3, 3/5

Here's one way.

Download TabPropQ.mw

(You may want to exclude 1 from k.)

This may be a helpful, though oversimplified example. As for visualization, once you have a solution, each lambda can be plotted in a 3-D plot against two of the parameters for fixed values of the other parameters. It may be possible by some sort of "dimensional" analysis to find combinations of the variables that reduce the number of variables/parameters.

restart;

with(plots):

Symmetry between lambda_2 and lambda__3

Eq2:=lambda__2-lambda__3^2;
Eq3:=lambda__3-lambda__2^2;

-lambda__3^2+lambda__2

-lambda__2^2+lambda__3

If we are interested only in real roots, we want the intersection points between the two plots - from the symmetry we will have lambda__2 = lambda__3 in this case

display(implicitplot(Eq2,lambda__2=-2..2,lambda__3=-2..2,color=red),
        implicitplot(Eq3,lambda__2=-2..2,lambda__3=-2..2,color=blue));

For this case we could have just solved one equation after recognizing lambda__2=lambbda__3

eval(Eq2,lambda__2=lambda__3);
solve(%,lambda__3);
# or:
solve({Eq2,lambda__2=lambda__3});

-lambda__3^2+lambda__3

0, 1

{lambda__2 = 0, lambda__3 = 0}, {lambda__2 = 1, lambda__3 = 1}

In the general case, we have cases where lambda__2 is not equal to lambda__3

ans:=solve({Eq2,Eq3},{lambda__2,lambda__3});
av:=allvalues([ans]);

{lambda__2 = 0, lambda__3 = 0}, {lambda__2 = 1, lambda__3 = 1}, {lambda__2 = -RootOf(_Z^2+_Z+1)-1, lambda__3 = RootOf(_Z^2+_Z+1)}

[{lambda__2 = 0, lambda__3 = 0}, {lambda__2 = 1, lambda__3 = 1}, {lambda__2 = -1/2-((1/2)*I)*3^(1/2), lambda__3 = -1/2+((1/2)*I)*3^(1/2)}], [{lambda__2 = 0, lambda__3 = 0}, {lambda__2 = 1, lambda__3 = 1}, {lambda__2 = -1/2+((1/2)*I)*3^(1/2), lambda__3 = -1/2-((1/2)*I)*3^(1/2)}]

simplify(eval([Eq2,Eq3],av[1][3]));
simplify(eval([Eq2,Eq3],av[2][3]));

[0, 0]

[0, 0]

Stepwise solve

lambda__3_1,lambda__3_2:=solve(Eq2,lambda__3);

lambda__2^(1/2), -lambda__2^(1/2)

eval(Eq3,lambda__3=lambda__3_1); # more complicated equation to solve in last step
solve(%,lambda__2);

-lambda__2^2+lambda__2^(1/2)

0, 1

eval(Eq3,lambda__3=lambda__3_2); # more complicated equation to solve in last step
solve(%,lambda__2);

-lambda__2^2-lambda__2^(1/2)

0, -1/2-((1/2)*I)*3^(1/2), -1/2+((1/2)*I)*3^(1/2)

NULL

Download symmetry.mw

Heres another approach for more recent Maple versions in 1-D only. It could be worked up into a procedure as @mmcdara did, but I am short of time now. I think it should be reasonably efficient but I didn't do any tests.

restart

M := Matrix(4, 3, {(1, 1) = 34, (1, 2) = 67, (1, 3) = 1, (2, 1) = 35, (2, 2) = 80, (2, 3) = 1, (3, 1) = 45, (3, 2) = 45, (3, 3) = 2, (4, 1) = 56, (4, 2) = 99, (4, 3) = 2})

Search column 3 for values 2 and return corresponding rows

M[[for i to upperbound(M)[1] do `if`(M[i,3]=2,i,NULL) end do]];

Matrix(2, 3, {(1, 1) = 45, (1, 2) = 45, (1, 3) = 2, (2, 1) = 56, (2, 2) = 99, (2, 3) = 2})

Search row 3 for values 45 and return corresponding columns

M[..,[for i to upperbound(M)[2] do `if`(M[3,i]=45,i,NULL) end do]];

Matrix(4, 2, {(1, 1) = 34, (1, 2) = 67, (2, 1) = 35, (2, 2) = 80, (3, 1) = 45, (3, 2) = 45, (4, 1) = 56, (4, 2) = 99})

NULL

Download search.mw

As @Carl Love says, it is common to get these small imaginary parts. They can be removed with fnormal and simplify(...,zero). Your case is unusual in that the digits for fnormal needs to be significantly different from the digits for the calculation - I don't recall ever having to use the digits option to fnormal before. Presumably the inaccuracies in the dilog, etc functions accumulate errors.

[Edit: the imaginary part looks suspiciously like Pi times a power of 10, so perhaps there is more going on here.]

restart;

with(Units) :

c := 299792458*Unit(m)/Unit(s) :
h := 662607015*10^(-8)*10^(-34)*Unit(J)/Unit(Hz) :
k := 1380649*10^(-6)*10^(-23)*Unit(J)/Unit(K) :
lambda_V_max := 750*Unit(nm) :
lambda_V_min := 380*Unit(nm) :
T_Sol := 5772.0*Unit(K) :
T0_Sol := 5772*Unit(K) :

M := (T, lambda_min, lambda_max) -> Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T)) - 1)*lambda^5), lambda = lambda_min .. lambda_max) :
sigma := 2*GAMMA(4)*Pi*Zeta(4)*k^4/(c^2*h^3) :

M_inf := (T) -> T^4*sigma :

M_Sol := M_inf(T_Sol) :
M0_Sol := M_inf(T0_Sol) : # 6.2938592470335950467548474587300E7*Unit(W)/Unit(m)^2

exact:=M(T0_Sol, lambda_V_min, lambda_V_max): # this is an exact value (zero uncertainty)

Automatically loading the Units[Simple] subpackage
 

approx:=evalf[16](exact);

(27551199.57168703-0.3141592653589793e-5*I)*Units:-Unit(kg/s^3)

fnormal should convert the imaginary part to 0.*I if it is just numerical roundoff - using the same number of digits does not do this

 - suggests the many complicated functions have led to more roundoff than usual.

fnormal(approx,16);

(27551199.57168703-0.3141592653589793e-5*I)*Units:-Unit(kg/s^3)

Try doing the calc to lower precision. Now the zero imaginary part may be removed with simplify/zero

fnormal(approx,7);
simplify(%,zero);

(0.2755120e8-0.*I)*Units:-Unit(kg/s^3)

0.2755120e8*Units:-Unit(kg/s^3)

approx:=evalf[32](exact);
fnormal(approx,23);
simplify(%,zero);

(27551199.571700602410253741952027+0.50265482457436691815402294132472e-21*I)*Units:-Unit(kg/s^3)

(27551199.571700602410254+0.*I)*Units:-Unit(kg/s^3)

27551199.571700602410254*Units:-Unit(kg/s^3)

NULL

Download fnormal.mw

for m from 0 to 5 do
  printf("% 5d\n",<seq(binomial(m,i),i=0..m)>);
end do;
   

    1
    1     1
    1     2     1
    1     3     3     1
    1     4     6     4     1
    1     5    10    10     5     1

As in @nm's answer, I'm not quite there yet, but have resolved the issue of which is the correct root

restart;

expected:=(1+cos(arctan(3/4)/3))/3;evalf(%);

1/3+(1/3)*cos((1/3)*arctan(3/4))

.6590276223

alias(alpha=RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1));

alpha

evalf(alpha);

.1090390091

alpha1:=simplify(convert(alpha,radical));
evalf(%);

(1/72)*((36*I)*(1/270+(1/360)*I)^(2/3)*3^(1/2)-36*(1/270+(1/360)*I)^(2/3)-I*3^(1/2)+24*(1/270+(1/360)*I)^(1/3)-1)/(1/270+(1/360)*I)^(1/3)

.1090390092-0.2069494425e-10*I

As written (both index=1), it is the wrong root to correspond to the expected result.

u1 := RootOf(4*_Z^2 + (4*RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1) - 4)*_Z + 4*RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1)^2 - 4*RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1) + 1);
evalf(%);

RootOf(4*_Z^2+(4*alpha-4)*_Z+4*alpha^2-4*alpha+1)

.2319333685

We want index=2 of the outer RootOf, but I copied here, which is cheating

RootOf(4*_Z^2 + (4*RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1) - 4)*_Z + 4*RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1)^2 - 4*RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1) + 1,index=2);
evalf(%);

RootOf(4*_Z^2+(4*alpha-4)*_Z+4*alpha^2-4*alpha+1, index = 2)

.6590276225

u2:=subs(alpha=a,u1);
rt1,rt2:=allvalues(u2);

RootOf(4*_Z^2+(4*a-4)*_Z+4*a^2-4*a+1)

-(1/2)*a+1/2+(1/2)*(-3*a^2+2*a)^(1/2), -(1/2)*a+1/2-(1/2)*(-3*a^2+2*a)^(1/2)

We want rt1

evalf(eval(rt1,a=alpha1));
evalf(eval(rt2,a=alpha1));

.6590276224-0.381382510e-11*I

.2319333684+0.2450876934e-10*I

u4:=simplify(eval(rt1,a=alpha1));
evalf(%);

1/3+(1/12)*cos((1/3)*arctan(3/4))+(1/12)*(cos((2/3)*arctan(3/4))+2-3^(1/2)*sin((2/3)*arctan(3/4)))^(1/2)*3^(1/2)+(1/12)*sin((1/3)*arctan(3/4))*3^(1/2)

.6590276224

Stuck on the last step...

 

Download RootOf.mw

You can work around this by making a new context for the are unit in which any SI prefix is allowed.

restart

with(Units[Standard])

A := Unit('ha')

Units:-Unit(ha)

convert(A, units, Unit(a))

100*Units:-Unit(a)

convert(A, units, Unit(Pa))

Error, (in convert/units) unable to convert `ha` to `Pa`

a and ha are defined as "one-offs" with prefix=none, i.e.,can't use prefixes.

Units:-GetUnit(a)

are, abbreviations = {}, plural = are, default = true, spellings = {are}, prefix = none, abbreviation = none, spelling = are, conversion = 100*metre[SI]^2, context = SI, symbols = {a}, symbol = a

Units:-GetUnit(ha)

hectare, abbreviations = {}, plural = hectares, default = true, spellings = {hectare, hectares}, prefix = none, abbreviation = none, spelling = hectare, conversion = 10000*metre[SI]^2, context = SI, symbols = {ha}, symbol = ha

But if we define a new context "realestate" then we can define a again with a different context in which we can use any prefix

Units:-UseContexts(realestate)

Units:-AddUnit(are, context = realestate, prefix = SI, conversion = are)

B := Unit(ha[realestate])

Units:-Unit(ha)

convert(B, units, Unit(Pa[realestate]))

(1/10000000000000)*Units:-Unit(Pa)

NULL

Download Areas.mw

I'm assuming you actually want to calculate these for specific examples. Then the group ring elements could be vectors of the ring elements, with each component associated with the corresponding group element. I left the ring multiplication generic and non-commutative, but you could define &* to be a more specific operation if you want. Likewise the group operation could instead be more generic, say ".", but then you wouldn't be able to do much. This worksheet does one of the Wikipedia page examples, but without assuming a commutative ring product.

restart

with(GroupTheory)

G := Group(a = Perm([[1, 2, 3]])); n := GroupOrder(G)

_m2088506322848

3

DrawCayleyTable(G, size = [200, 200], assignelements = 'els')

Element order for the vectors

els; elsindex := table(`~`[`=`](els, [`$`(1 .. n)]))

[_m2088475132256, _m2088506358112, _m2088527908928]

Ring multiplication is &*. Make it associative (flat), with Identity E.

define('`&*`', 'flat', 'identity' = E)

Associative and distributive

expand(`&*`(a+e, b+c)); `&*`(a, `&*`(b, c))-`&*`(`&*`(a, b), c)

`&*`(a, b)+`&*`(a, c)+`&*`(e, b)+`&*`(e, c)

0

Define group-ring product that operates on n-vectors of ring elements. First component goes with first group element etc.

`&otimes;` := proc (A, B) options operator, arrow; `<,>`(seq(add(`&*`(A[i], B[elsindex[PermInverse(els[i]).els[j]]]), i = 1 .. n), j = 1 .. n)) end proc

Warning, (in &otimes;) `j` is implicitly declared local

Warning, (in &otimes;) `i` is implicitly declared local

r := `<,>`(z__0, z__1, z__2); s := `<,>`(w__0, w__1, w__2)

Vector(3, {(1) = z__0, (2) = z__1, (3) = z__2})

Vector[column](%id = 36893490236068004316)

`&otimes;`(r, s)

Vector[column](%id = 36893490236068007076)

r+s

Vector[column](%id = 36893490236067996964)

Download group-ring.mw

I don't see a way to do this directly, but a workaround is just to save the string to a local file and then read it.

restart

filnam := cat(currentdir(), "/test.mpl")

"C:\Users\dharr\Desktop/test.mpl"

FileTools[Text]:-WriteFile(filnam, Import("https://www.nicolesharp.net/testbox/test.mpl"))

227

read filnam

Error, (in alias) equations expected

logpi(3)

ln(3)/ln(pi)

Download MPL.mw

By making FromNames2RV a procedure, you can achieved the result you want.

FromNames2RV := () -> [anames(RandomVariable)] =~ eval([anames(RandomVariable)]);

Then using FromNames2RV() multiple times gives the same (expected) result.

FromNames2RV := () -> map(x -> x = eval(x), [anames(RandomVariable)]);

is an alternative that doesn't call anames twice.

If you do the factoring step by step, each step gets slower, presumably because of the deeper nesting of the RootOfs.So Afactor and Split presumably do it in this same way. Since the polynomials are irreducible at each stage, this is the slowest case. I think galois is fast because Maple just has a database of cases for low-degree polynomials. Perhaps Maple could make more use of the symmetry, but I don't know much about this.

Download galois.mw

SubgroupLattice reports 1431 subgroups, but gives a warning that perfect subgroups may be missing

restart;

with(GroupTheory):

G:=SymmetricGroup(6);

_m2603493461024

L:=SubgroupLattice(G);

Warning, group is not solvable and soluble residual has order >= 360; perfect subgroups may be missing from lattice.

"GroupTheory:-SubgroupLattice(module() ... end module)"

List:=convert(L,list):nops(%);

1431

Choose a group - generator is given as shorthand

G2:=List[321];

_m2603594640096

Elements(G2);

{_m2603495175360, _m2603495644000, _m2603599152640, _m2603599172544}

 

Download S6.mw

Aside from maplemint, in more recent Maple versions if you write code in the startup code region, you get similar information in the diagnostics panel. For example running the code

Q:=proc(x) dsolve(y,method=rkf45) end;

gives:
Info: Line 1: These parameters were never used: x
Info: Line 1: These names were used as global names but were not declared: y
Info: Line 1: These names were used as global names but were not declared: method
Info: Line 1: These names were used as global names but were not declared: rkf45

For future reference, please upload your worksheet using the green up-arrow in the Mapleprimes editor. Maple does give that complicated solution, so you can put in numerical parameters and make various plots, even if it is not so easy to view the solution. If you want a better view of the structure of the solution, here is one way to analyze that.

Download solve.mw

Edit: the solution may look simpler depending on the order in which Maple goes about solving the system. If manipulating this to the simplest form is important, you can experiment by specifying the variable order that PolynomialSystem uses by putting the variables in a list (in a set as here Maple tries to find a nice order). But Maple can handle the complicated solution so long as you don't need to print it out, so a lot depends on what you want to do with the answer.

First 23 24 25 26 27 28 29 Last Page 25 of 83