Advanced Search

dharr

Dr. David Harrington

8762 Reputation

22 Badges

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

Error in paper...

dharr 8762
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 8762
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 8762
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 8762
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 8762
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 8762
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 8762
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 8762
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 8762
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 8762
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 8762
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 8762
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

No polynomial solutions...

dharr 8762
November 04 2025
2 0

The Heunc help page has some detailed instructions for finding polynomial solutions.

> 

restart;
> 

local gamma:
> 

Hc := HeunC(-4*b, 4*b, 2, 8*b^2, -8*b^2 - _c[2] + 1, (-r + 2*b)/(2*b));

HeunC(-4*b, 4*b, 2, 8*b^2, -8*b^2-_c[2]+1, (1/2)*(-r+2*b)/b)

> 

alpha, beta, gamma, delta, eta := op(1..5, Hc);

-4*b, 4*b, 2, 8*b^2, -8*b^2-_c[2]+1

A necessary condition for a polynomial solution of degree n is (see the HeunC help page)

> 

poly_cond := delta = -(n + (gamma + beta + 2)/2)*alpha;

8*b^2 = 4*(n+2+2*b)*b

So b has to be (for any n)

> 

solve(poly_cond, b);

0

But the last argument has b in the denominator, so b=0 does not work, and there can be no polynomial solutions

> 

eval(Hc,b=0);

Error, numeric exception: division by zero

> 

NULL

Download polysols.mw

no solution...

dharr 8762
November 02 2025
1 2

There doesn't seem to be a solution, assuming the data you give is exactly correct. The presence of 1+A*3.7653333333333333333*10^(-18)*nu needs very large A for the types of nu values you have, otherwise the A*c*nu term is too small compared to 1 to make a difference.

MLsolve.mw

For the followup examples...

dharr 8762
November 02 2025
1 5

[Edit: I give a more general answer below.]

> 

restart;
> 

local gamma;

gamma

If we know that a is in the field generated by b and we have each of them written as RootOf  a single polynomial

> 

a := RootOf(4*_Z^3 - 12*_Z^2 - 24*_Z - 25, index = 1); # which actually equals “(b + 1)^2”, but assuming we do not know this beforehand
b := RootOf(2*_Z^3 - 3, index = 1);

RootOf(4*_Z^3-12*_Z^2-24*_Z-25, index = 1)

RootOf(2*_Z^3-3, index = 1)

> 

alias(alpha=a,beta=b);

alpha, beta

Then we "just" evaluate a specifying b as a field extension. Unfortunately, many of the evala subcommands such as Normal deduce the field extension from the coefficients, and don't allow us to specify it. Factor allows us to specify it, and even though we don't need to actually factor here it will serve the purpose:

> 

evala(Factor(a,b));

beta^2+2*beta+1

We don't need Primfield, which just confirms that a can be written in terms of b, so the field extension can just be b

> 

evala(Primfield({a,b}));

[[beta = beta], [beta = beta]]

Algfield works similarly, but gives the relationship between a and b

> 

evala(Algfield({a,b}));

[[alpha = beta^2+2*beta+1], [], {beta}, true, {}]

In the second case, Maple already simplifies each of these in terms of RootOf(_Z^2-2, index = 1) and RootOf(_Z^2+1, index = 1), which doesn't help us here.

> 

c:= convert((-1)^(1/4), RootOf);
evalf(%);
d:= convert((sqrt(2) - 1)*I, RootOf); # which actually equals c*(c^2 - c + 1), but assuming we do not know this beforehand
evalf(%);

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

.7071067810+.7071067810*I

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

.4142135620*I

Convert these to forms that are RootOf a single polynomial

> 

mp:=evala(Minpoly(c,x)):
c:=select(q->evala(c-q)=0,[seq(RootOf(mp,index=i),i=1..degree(mp,x))])[];
evalf(%);
mp:=evala(Minpoly(d,x)):
d:=select(q->evala(d-q)=0,[seq(RootOf(mp,index=i),i=1..degree(mp,x))])[];
evalf(%);

RootOf(_Z^4+1, index = 1)

.7071067812+.7071067812*I

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

.4142135624*I

> 

alias(gamma=c, delta=d);

alpha, beta, gamma, delta

> 

evala(Factor(d,c));

gamma^3-gamma^2+gamma

In this example, we can write RootOf(_Z^4+6*_Z^2+1) in terms of  RootOf(_Z^4+1, index = 1) The method doesn't work the other way around - it just returns c:

> 

evala(Factor(c,d));

gamma

Same problem for the first one if we want to do it the other way around.

So this method also does not work on the example in the original question

> 

evala(Factor(b,a));

beta

> 

NULL

Download evala.mw

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