Advanced Search

dharr

Dr. David Harrington

9097 Reputation

22 Badges

21 years, 229 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

Collecting with symbolic exponents...

dharr 9097
December 07 2025
4 6

Edit: handling of terms like ln(r) improved.

https://www.mapleprimes.com/questions/242062-Collect-Using-Pattern

> 

restart;

We want to collect A with respect to r, including symbolic exponents, in order to get B. I've added a constant term to the OP's example, to further test the method.
It is not clear from the question if the order matters, or what it would be for more complicated exponents.

For simplicity, assume that the order of exponents is just the default sort order, which would work here.

> 

A:=r^(2*a)+r^2*(1+a+r^(2*a)) + r + a*r + a;
B:=a+(1+a)*r+(1+a)*r^2+r^(2*a)+r^(2+2*a);

r^(2*a)+r^2*(1+a+r^(2*a))+r+a*r+a

a+(1+a)*r+(1+a)*r^2+r^(2*a)+r^(2+2*a)

Notice that Maple's automatic simplification may frustrate any attempt to get a desired order, e.g., the following is already collected, but does not retain the input order, and is perhaps not in the expected order

> 

r^3 + a*r^2 + r;

a*r^2+r^3+r

The reason collect does not work is that it is for polynomials, which have integer exponents but not symbolic ones.
The following finds the exponent wrt r for any exponent not containing r

> 

ex := (y,r) -> r*diff(y,r)/y;

proc (y, r) options operator, arrow; r*(diff(y, r))/y end proc

Examples. The presence of r in the second case indicates we should probably just leave this alone (this is done in the procedure below).

> 

ex(7*r^(2*tan(a)), r);
ex(a*ln(r), r);

2*tan(a)

1/ln(r)

Expand, find the terms and their corresponding exponents and coefficients. (assumes it is type `+`)

> 

terms := [op(expand(A))];
exps := expand(ex~(terms, r));
cfs := zip((t,e)->expand(t/r^e), terms, exps);

[(r^a)^2, r^2, a*r^2, r^2*(r^a)^2, r, a*r, a]

[2*a, 2, 2, 2+2*a, 1, 1, 0]

[1, 1, a, 1, 1, a, a]

Sort them

> 

exps2, p := sort(exps, output = [sorted, permutation]);
cfs2 := cfs[p];

[0, 1, 1, 2, 2, 2*a, 2+2*a], [7, 5, 6, 2, 3, 1, 4]

[a, 1, a, 1, a, 1, 1]

Assemble the results

> 

s:=0: cf:=cfs2[1]: expt:=exps2[1]:
for i from 2 to nops(exps2) do
  if exps2[i] = expt then
    cf += cfs2[i]
  else
    s += cf*r^expt;
    expt := exps2[i];
    cf := cfs2[i]
  end if;
end do;
s += cf*r^expt;

a+(1+a)*r+(1+a)*r^2+r^(2*a)+r^(2+2*a)

As a procedure.

> 

Collect := proc(expr, x::name)
local zz, ex, term, cfs, exps, p, s, cf, i;
if expr::function then return expr end if;
for i,term in expand(expr)+zz do # zz forces type `+`
  ex := expand(x*diff(term,x)/term);
  if has(ex, x) then
     cfs[i] := term;
     exps[i] := 0;
  else
     exps[i] := ex;
     cfs[i] := expand(term/x^ex)
  end if;
end do;
exps := convert(exps, list);
cfs := convert(cfs, list);
exps, p := sort(exps, 'output' = ['sorted', 'permutation']);
cfs := cfs[p];
s := 0;
cf := cfs[1];
ex := exps[1]:
for i from 2 to nops(exps) do
  if exps[i] = ex then
    cf += cfs[i]
  else
    s += cf*x^ex;
    ex := exps[i];
    cf := cfs[i]
  end if;
end do;
s += cf*x^ex - zz;
end proc:

Test

> 

Collect(A + ln(r), r);

a+ln(r)+(1+a)*r+(1+a)*r^2+r^(2*a)+r^(2+2*a)

> 

Collect(r^3 + a*r^2 + r, r);

a*r^2+r^3+r

> 

Collect(ln(x)+3*x^2, x);

ln(x)+3*x^2

x^f(x) is left alone, i.e., is treated as a constant term.

> 

Collect(3*x^2 + 2*x^x, x);

3*x^2+2*x^x

> 

 

Download CollectSymbolicExponents.mw

simplify with side relations...

dharr 9097
December 04 2025
1 3 
sqrt(1-cos(x)^2);
simplify(%,{cos(x)^2=1-sin(x)^2});

expand...

dharr 9097
December 01 2025
1 0

For some reason, the original simplify gives a numerator/denominator, with sqrt(1-t^2) in both.  If expand is used first then we don't have numerator/denominator, Then simplify can deal with the diffs. As well the denominator sqrt(1-t^2) might cancel or be simpler. So simplify(expand(B)) works here. The collect serves the same purpose as expand - giving priority to the diffs means we get a polynomial in diffs and not numerator/denominator, which then simplifies. The collect can be just wrt diff, so simplify(collect(B,diff))works.

Error in paper...

dharr 9097
November 23 2025
0 1

As far as I can see, there is some error in the paper. Maybe one of the other case(s) works.

Download pde.mw

Numbers...

dharr 9097
November 19 2025
3 1

For the numbers.

> 

restart

On a 64 bit machine

> 

kernelopts(wordsize);

64

> 

evalhf(DBL_MIN);  #$MinMachineNumber

0.100000010000000008e-306

> 

evalhf(DBL_MAX);  #$MaxMachineNumber

0.99999999000000001e308

> 

Maple_floats(MAX_FLOAT); #$MaxNumber

0.1e9223372036854775807

> 

Maple_floats(MIN_FLOAT); #$MinNumber

0.1e-9223372036854775805

> 

NULL

Download Numbers.mw

For the calling of previous commands by referring to them by the equation labels on the right-hand side, use ctrl-L to insert the number in the appropriate place.

open panel/palettes...

dharr 9097
November 18 2025
3 1

If you seach for "slow" on this site, there are mentions that having the context panel - the panel on the right - closed can increase performance; you might also try closing the palettes on the left when you are not using them.

Elliptic curve solutions...

dharr 9097
November 16 2025
4 0

[Edited]

parametrization is good for finding the first rational point (not necessarily integer) on genus 0, but doesn't work for genus 1.

I assumed at first you wanted integer solutions. My usual attacks here were not that successful. For the first equation 1 found 1 more solution, but rather trivial, and then the addition process gives the infinity point. For the second case, it seems to be working nicely in producing rational solutions.

Elliptic.mw (first equation)

Second equation:

https://www.mapleprimes.com/questions/242001-Elliptic-Curves

> 

restart
> 

with(algcurves)

Find integer solutions of the following equation of interest; we know one solution is (3,7);

> 

eq := y^2-2*y+14 = 2*x^3+11*x^2-29*x-17; f := (rhs-lhs)(eq)

y^2-2*y+14 = 2*x^3+11*x^2-29*x-17

2*x^3+11*x^2-y^2-29*x+2*y-31

> 

sol1 := {x = 3, y = 7}; eval(f, sol1)

{x = 3, y = 7}

0

Check genus is 1

> 

genus(f, x, y)

1

Define addition of points in the group algebra - see Wikipedia

A is the coefficient in y^2=x^3+A*x+B

> 

algadd := proc(P::[rational,{rational,identical(infinity)}], Q::[rational,{rational,identical(infinity)}], A::rational)
# infinity point denoted [0, infinity]
  local s,xP,yP,xQ,yQ,xR,yR;
  xP,yP := P[];
  xQ,yQ := Q[];
  if yQ = infinity and yQ = infinity then return [0, infinity] end if;
  if yP = infinity then return [xQ, yQ] end if; # infinity point is identity
  if yQ = infinity then return [xP, yP] end if;
  if xP <> xQ then
    s := (yP - yQ)/(xP - xQ);
    xR := s^2 - xP - xQ;
    yR := yP - s*(xP - xR);
    [xR, -yR];
  elif yP = -yQ then
    [0, infinity] # infinity point
  else # P = Q,sum is tangent
    s := (3*xP^2 + A)/(2*yP);
    xR := s^2 - 2*xP;
    yR := yP - s*(xP - xR);
    [xR, -yR]
  end if;
end proc:

The j_invariant has the same value despite the various transformations (birational equivalence)

> 

j_invariant(f, x, y)

1643032000/703913

Maple finds a Weierstrass form

> 

wform := Weierstrassform(f, x, y, X__1, Y__1); w1 := wform[1]; j_invariant(w1, X__1, Y__1)

X__1^3-(295/3)*X__1-5164/27+Y__1^2

1643032000/703913

The addition routine wants it in the standard form y^2 = x^3+A*x+B. For integer solutions, take A and B integers, which we can achieve by changing the sign of X, scaling the whole thing by 9*81 and redefine X and Y by factor of 9 and 27. (For rational solutions, rational A and B are OK, so we only need to change the sign of X here.)

> 

trans := {X__1 = -(1/9)*X, Y__1 = (1/27)*Y}; w := eval((9*81)*w1, trans); j_invariant(w, X, Y)

{X__1 = -(1/9)*X, Y__1 = (1/27)*Y}

-X^3+Y^2+7965*X-139428

1643032000/703913

Conversion between f and w

> 

itrans := solve(trans, {X, Y}); xy_to_XY := eval(itrans, {X__1 = wform[2], Y__1 = wform[3]})

{X = -9*X__1, Y = 27*Y__1}

{X = 18*x+33, Y = 54*y-54}

Initial solution in X,Y

> 

sol1XY := eval(xy_to_XY, sol1); eval(w, sol1XY)

{X = 87, Y = 324}

0

Conversion between w and f

> 

XY_to_xy := eval({x = wform[4], y = wform[5]}, trans); eval(XY_to_xy, sol1XY)

{x = (1/18)*X-11/6, y = 1+(1/54)*Y}

{x = 3, y = 7}

To start the algorithm, we need A and the first solution

> 

isolate(w, Y^2); A := coeff(rhs(%), X, 1); P[1] := eval([X, Y], sol1XY)

Y^2 = X^3-7965*X+139428

-7965

[87, 324]

Add P[1] to get another solution - it's rational but not integer

> 

P[2] := algadd(P[1], P[1], A); eval(XY_to_xy, `~`[`=`]([X, Y], %)); eval(f, %)

[5497/16, -394291/64]

{x = 4969/288, y = -390835/3456}

0

> 

P[3] := algadd(P[2], P[1], A); eval(XY_to_xy, `~`[`=`]([X, Y], %)); eval(f, %)

[3510043719/16851025, 189667068200172/69173457625]

{x = 164108883/16851025, y = 3581526572443/69173457625}

0

> 

P[4] := algadd(P[3], P[1], A); eval(XY_to_xy, `~`[`=`]([X, Y], %)); eval(f, %)

[1015338131548513/9949785131584, -19581792761066942460817/31384885834539095552]

{x = 686995222206241/179096132368512, y = -17887008926001831301009/1694783835065111159808}

0

To find the first solution, we want either Y = 0, or a Y value for which Y^2 divides the disciminant, for which the corresponding X is an integer  (Nagell-Lutz) .

> 

Cubic := rhs(isolate(w, Y^2)); d := discrim(Cubic, X); ifactor(d)

X^3-7965*X+139428

1496352914532

``(2)^2*``(3)^12*``(7)*``(100559)

So test 0 and all the divisors of 2*3^6

> 

divs := `union`({0}, NumberTheory:-Divisors(mul(map(proc (q) options operator, arrow; ifelse((q[2])::even, q[1]^((1/2)*q[2]), 1) end proc, ifactors(d)[2]))))

{0, 1, 2, 3, 6, 9, 18, 27, 54, 81, 162, 243, 486, 729, 1458}

> 

j := 0; for i in divs do facs := factors(eval(w, Y = i))[2]; if nops(facs) > 1 then j := j+1; sol[j] := {X = solve(select(proc (z) options operator, arrow; type(z[1], linear) end proc, facs)[1][1], X), Y = i}; print(sol[j]) end if end do

No luck here.

> 

NULL

Download Elliptic2.mw

timelimit...

dharr 9097
November 15 2025
2 10

A few of points:

1. Sometimes a calulation can be fast but displaying the complicated expression can be slow, so it is better then to store the answer without displaying it, e.g. use ans := pdetest(eqv2, pde): then if it completes you can attempt to display it (perhaps check its length first).

2. You can use ans := timelimit(20, pdetest(eqv2, pde)): to limit the attempt to 20 s.

3. If you suspect a complicated expression might simplify to zero you can set random numerical values for parameters and see if the numerical value is nearly zero. In your case you have functions beta(t) etc, so you might give simple functions to these, e.g.,

params := {beta=(t->t^2), k[1] = 0.35, ...};
ans := eval(eval(pde, params), eval(eqv2, params)):

 

Converting to a first order system...

dharr 9097
November 14 2025
0 3

I didn't undertand what you said. It looks as though you were trying to manually make a first order system. I fixed this and get the same result as convertsys.

systems.mw

bad parameter values...

dharr 9097
November 13 2025
2 5
> 

restart
> 

eq1 := S__1 = sqrt((-beta[1]+sqrt(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2))/beta[2]); eq2 := S__2 = (1/2)*sqrt(-(2*(beta[1]+sqrt(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)))/beta[2]); eq3 := S__3 = sqrt(2)*sqrt(chi*(4*chi^2*k^2*p+4*chi^2*w+2*chi*p*beta[1]-p^2*beta[2]))/(4*chi^2); eq4 := T__1 = sqrt(-2*chi*p)/(2*chi); eqs := {eq1, eq2, eq3, eq4}

S__1 = ((-beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2)

S__2 = (1/2)*(-2*(beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2)

S__3 = (1/4)*2^(1/2)*(chi*(4*chi^2*k^2*p+4*chi^2*w+2*chi*p*beta[1]-p^2*beta[2]))^(1/2)/chi^2

T__1 = (1/2)*(-2*chi*p)^(1/2)/chi

{S__1 = ((-beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2), S__2 = (1/2)*(-2*(beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2), S__3 = (1/4)*2^(1/2)*(chi*(4*chi^2*k^2*p+4*chi^2*w+2*chi*p*beta[1]-p^2*beta[2]))^(1/2)/chi^2, T__1 = (1/2)*(-2*chi*p)^(1/2)/chi}

> 

indets(eqs)

{S__1, S__2, S__3, T__1, chi, k, p, w, beta[1], beta[2], (chi*(4*chi^2*k^2*p+4*chi^2*w+2*chi*p*beta[1]-p^2*beta[2]))^(1/2), ((-beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2), (-2*chi*p)^(1/2), (-2*(beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2), (4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2)}

We have 4 eqns, so choose 4  variables to solve for

> 

sol := solve(eqs, {p, w, beta[1], beta[2]})

{p = -2*T__1^2*chi, w = 2*chi*(S__1^2*S__2^2*T__1^2*k^2-S__1^2*T__1^4*k^2-2*S__2^2*T__1^4*k^2+2*T__1^6*k^2+S__1^2*S__2^2*S__3^2)/(S__1^2*S__2^2-S__1^2*T__1^2-2*S__2^2*T__1^2+2*T__1^4), beta[1] = 2*S__3^2*chi*(S__1^2+2*S__2^2)/(S__1^2*S__2^2-S__1^2*T__1^2-2*S__2^2*T__1^2+2*T__1^4), beta[2] = -4*S__3^2*chi/(S__1^2*S__2^2-S__1^2*T__1^2-2*S__2^2*T__1^2+2*T__1^4)}

Put values in - we see a problem

> 

eval(sol, {S__1 = 1, S__2 = 1, S__3 = 1, T__1 = 1, k = 1})

Error, numeric exception: division by zero

For some other values it should work

> 

eval(sol, {S__1 = 1, S__2 = 1, S__3 = 1, T__1 = 2, k = 1})

{p = -8*chi, w = (170/21)*chi, beta[1] = (2/7)*chi, beta[2] = -(4/21)*chi}

NULL

Download find-parameter.mw

Bug...

dharr 9097
November 12 2025
1 0

Someone else found a similar problem here and I reported this bug to Maple. The same workaround might work here - put simplify(...) around the expression with the square root.

output=Array...

dharr 9097
November 11 2025
1 4

You can use the output=Array option of dsolve to do this.

> 

restart

Specify required x values

> 

xpts := Array([seq(0 .. 3.0, numelems = 20)])

Array(%id = 36893490082029227172)

> 

ode := {diff(y(x), x) = -a*y(x), y(0) = 1}

{diff(y(x), x) = -a*y(x), y(0) = 1}

> 

NULL
> 

ans20 := dsolve(eval(ode, a = 2.0), numeric, output = xpts); ans21 := dsolve(eval(ode, a = 2.1), numeric, output = xpts)

extract vectors of x and y values

> 

xx := ans20[2, 1][() .. (), 1]; y20 := ans20[2, 1][() .. (), 2]; y21 := ans21[2, 1][() .. (), 2]

Plot the difference

> 

plot(xx, y21-y20)

Plot the ratio (enter \/~)

> 

plot(xx, `~`[`/`](y21, y20))

Same thing using the new (Maple 2024) elementwise

> 

plot(xx, elementwise(y21/y20))

> 

NULL

Download dsolveArray.mw

Integration constant...

dharr 9097
November 09 2025
4 3

If you expand Maple's result you get (1/4)*x^4 + x^2 + 1, differing from your expected value by a constant. Since you and Maple are both omitting arbitrary constants, both answers are correct.

Number of BCs...

dharr 9097
November 08 2025
0 1

@KIRAN SAJJAN I would write the 6th bc as (D@@2)(f)(1) = 0. But the equations have a third order derivative of f and a second order derivative of theta, so Maple is right in only expecting 5 boundary conditions. So the paper has an error, perhaps in the odes or perhaps in the bcs (redundant one perhaps?, but I didn't check.)

I recommend putting the parameters in a set so you can look at the odes without the numbers in to check them. Then it is easier to find typos. I found two. The second lambda_val just before Bstar in the second equation should be omitted, and the first equation has an incorrect term. Corrected here:

thin_film_base_paper_comparision.mw

Successive substitutions...

dharr 9097
November 08 2025
2 6

[Updated] If I understand you correctly, then this classic method should work.

> 

restart;
> 

with(StringTools):
> 

s:="from -\\infty to \\infty and then from \\infty to -\\infty";

"from -\infty to \infty and then from \infty to -\infty"

> 

SubstituteAll(s, "-\\infty", "#@#"):
SubstituteAll(%, "\\infty", "+\\infty"):
s2:=SubstituteAll(%, "#@#","-\\infty");

"from -\infty to +\infty and then from +\infty to -\infty"

> 

NULL

Download Substitute.mw

1 2 3 4 5 6 7 Last Page 3 of 87