Advanced Search

dharr

Dr. David Harrington

8205 Reputation

22 Badges

20 years, 337 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
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

For simple units...

dharr 8205
January 27 2024
0 0

This is fine for simple units but is fairly limited for compound units which are just formatted as Maple expressions with the "*" etc. And for printing formatted strings, symbols such as Angstroms are expressed as their html code, so it would take a lot of work to format unit strings for the general case.

There are lines in the worksheet that are not displayed in Mapleprmes, e.g.,

printf("The voltage is %a %a", NumberUnit(x))

restart

with(Units:-Simple)

NULL

NumberUnit := proc(x)
 local num;
 num := convert(x, unit_free);
 num, eval(x/num, Units:-Unit = (() -> _passed));
end proc:

NULL

x := 5*Unit('V')

5*Units:-Unit(V)

NumberUnit(x)

5, V

printf("The voltage is %a %a", NumberUnit(x))

The voltage is 5 V

NULL

x := 5*Unit('kg'*'m'/('s'))

5*Units:-Unit(kg*m/s)

printf("The momentum is %a %a", NumberUnit(x))

The momentum is 5 kg*m/s

NULL

Download PrintUnits.mw

Here's a second version you may prefer

PrintUnits2.mw

printf("The electric field in base units is %5.2f %s", NumberUnit(x))

The electric field in base units is  5.00 m kg / A s^3

useunits...

dharr 8205
January 27 2024
1 0

The plot option useunits means the units are interpreted correctly. (And your f2 in 1D was missing a multiplication.)

FunctionUnits.mw

one way...

dharr 8205
January 27 2024
1 1

It could probably be made more efficient but I wasn't sure I understood the descriptions and all the assumptions.

Given:
How Can I sort a list of polynomials with parametric coefficients based on the following criteria?
1. Fewer different parameters in the coefficients.
2. If the number of parameters is the same, the lower power of the parameters.
3. If the previous criteria are the same, sentences with parametric coefficients appear later.

Assume: the variables are indexed and the parameters are not.
For the last condition, you mean
x[1]  + a*x[2] is after a*x[1] + x[2] because the a is later (with x[2] not x[1])
and the result depends on the variable indices rather than the order they are wriiten, so
a*x[2  + x[1] is the same as x[1]  + a*x[2] and is after a*x[1] + x[2]
The variable order is as returned by indets, so x[1], x[2] etc, but more complicated if we have x[1] and y[2] etc

> 

restart;

Correct answer for test:

> 

req:=[(a+2)*x[1]+3*x[2]-x[3]+(a-2)*x[4],2*x[1]+2*x[2]-b*x[3]+(c+1)*x[4], (a-1)*x[1]+x[2]+(a+1)*x[3]+(C^2-1)*x[4], (a+b)*x[1]+(c-3)*x[2]+(-b-1)*x[3]+2*x[4]];

[(a+2)*x[1]+3*x[2]-x[3]+(a-2)*x[4], 2*x[1]+2*x[2]-b*x[3]+(c+1)*x[4], (a-1)*x[1]+x[2]+(a+1)*x[3]+(C^2-1)*x[4], (a+b)*x[1]+(c-3)*x[2]+(-b-1)*x[3]+2*x[4]]

> 

test:=[(a-1)*x[1]+x[2]+(a+1)*x[3]+(C^2-1)*x[4], 2*x[1]+2*x[2]-b*x[3]+(c+1)*x[4], (a+2)*x[1]+3*x[2]-x[3]+(a-2)*x[4], (a+b)*x[1]+(c-3)*x[2]+(-b-1)*x[3]+2*x[4]];

[(a-1)*x[1]+x[2]+(a+1)*x[3]+(C^2-1)*x[4], 2*x[1]+2*x[2]-b*x[3]+(c+1)*x[4], (a+2)*x[1]+3*x[2]-x[3]+(a-2)*x[4], (a+b)*x[1]+(c-3)*x[2]+(-b-1)*x[3]+2*x[4]]

Comparison function

> 

comp:=proc(p1,p2)
  local vars, params, c1, c2, n1, n2, d1, d2;
  vars, params := selectremove(type,indets([p1,p2],name),indexed);
  c1 := [coeffs(p1,vars)];
  c2 := [coeffs(p2,vars)];
  n1 := nops(indets(c1)); # number of parameters in coeffs
  n2 := nops(indets(c2));
  d1 := max(map(degree, c1)); # max degree of parameters
  d2 := max(map(degree, c2));
  member(true,map(x->evalb(nops(indets(x))>0),c1),'pos1'); # pos1 = position of first param
  member(true,map(x->evalb(nops(indets(x))>0),c2),'pos2');
  if n2<n1 then
    false
  elif d2>d1 then
    false
  elif pos2>pos1 then
    false
  else
    true
  end if;
end proc:

Check

> 

sort(test,comp);

[(a+2)*x[1]+3*x[2]-x[3]+(a-2)*x[4], 2*x[1]+2*x[2]-b*x[3]+(c+1)*x[4], (a-1)*x[1]+x[2]+(a+1)*x[3]+(C^2-1)*x[4], (a+b)*x[1]+(c-3)*x[2]+(-b-1)*x[3]+2*x[4]]

> 

req;

[(a+2)*x[1]+3*x[2]-x[3]+(a-2)*x[4], 2*x[1]+2*x[2]-b*x[3]+(c+1)*x[4], (a-1)*x[1]+x[2]+(a+1)*x[3]+(C^2-1)*x[4], (a+b)*x[1]+(c-3)*x[2]+(-b-1)*x[3]+2*x[4]]

> 

NULL

 

Download sort2.mw

Calkin Wilf sequence...

dharr 8205
January 27 2024
2 1 
seq(NumberTheory:-CalkinWilfSequence(i),i=1..10)

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

one way...

dharr 8205
January 25 2024
1 4

Here's one way.

Download TabPropQ.mw

(You may want to exclude 1 from k.)

simpler example...

dharr 8205
January 23 2024
2 29

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

another way...

dharr 8205
January 22 2024
1 4

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

fnormal and simplify/zero...

dharr 8205
January 20 2024
3 12

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

one way...

dharr 8205
January 18 2024
1 0 
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

RootOf...

dharr 8205
January 17 2024
1 3

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

new context...

dharr 8205
January 16 2024
2 0

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

one way...

dharr 8205
January 16 2024
3 3

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

save to local file...

dharr 8205
January 14 2024
1 9

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

FromNames2RV procedure...

dharr 8205
January 13 2024
1 1

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.

nesting...

dharr 8205
January 12 2024
2 0

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

First 21 22 23 24 25 26 27 Last Page 23 of 81