Advanced Search

dharr

Dr. David Harrington

8697 Reputation

22 Badges

21 years, 80 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 replies submitted by dharr

Not a second order DE...

dharr 8697
4 hours ago
0 0

Your ode is only first order in Omega(xi), so converting it to a first-order system doesn't do anything. If you change it to second order, it works.

make_system.mw

extraneous solutions...

dharr 8697
Yesterday at 1:27 PM
0 0

@Math-dashti Maple has a habit of squaring both sides of an equation, which can lead to invalid solutions. It seems that there might not be a solution for the 6 variables you chose

> 

restart

Eqs written as expressions equal to zero.

> 

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); eq4 := T__1-sqrt(-2*chi*p)/(2*chi)

Equations - note we can't have beta[2] or chi zero

> 

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, T__1-(1/2)*(-2*chi*p)^(1/2)/chi}

> 

indets(eqs, name)

{S__1, S__2, S__3, T__1, chi, k, p, w, beta[1], beta[2]}

See if there is a solution for 6 chosen variables in terms of the others - yes, one solution with beta[2] and k as free parameters.

> 

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

{chi = -4*S__3^2/(beta[2]*(S__1^2*S__2^2-S__1^2*T__1^2-2*S__2^2*T__1^2+2*T__1^4)), k = k, p = 8*T__1^2*S__3^2/(beta[2]*(S__1^2*S__2^2-S__1^2*T__1^2-2*S__2^2*T__1^2+2*T__1^4)), w = -(1/2)*(S__1^4*S__2^4*beta[2]^2-S__1^4*S__2^2*T__1^2*beta[2]^2-2*S__1^2*S__2^4*T__1^2*beta[2]^2+2*S__1^2*S__2^2*T__1^4*beta[2]^2+16*S__3^2*T__1^2*k^2)/(beta[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] = -(1/2)*S__1^2*beta[2]-S__2^2*beta[2], beta[2] = beta[2]}

Check that they are solutions - not immediately zero unless S3 and T1 are zero. But S3 = 0 makes chi=0, which is a problem

> 

simplify(eval(eqs, sol), symbolic)

{0, 2*S__3, 2*T__1}

Can we put over common denominator and take the numerator - now chi is not in the denominator

> 

eqs2 := map(`@`(numer, normal), eqs)

{S__1-((-beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2), 2*S__2-(-2*(beta[1]+(4*k^2*p*beta[2]+4*w*beta[2]+beta[1]^2)^(1/2))/beta[2])^(1/2), -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)+4*S__3*chi, 2*T__1*chi-(-2*chi*p)^(1/2)}

We get the same solution as before and a new one, which solves the new equations but not the originals

> 

sol2 := [solve(eqs2, {chi, k, p, w, beta[1], beta[2]})]

[{chi = 0, k = k, p = p, w = -(1/2)*S__1^2*S__2^2*beta[2]-p*k^2, beta[1] = -(1/2)*S__1^2*beta[2]-S__2^2*beta[2], beta[2] = beta[2]}, {chi = -4*S__3^2/(beta[2]*(S__1^2*S__2^2-S__1^2*T__1^2-2*S__2^2*T__1^2+2*T__1^4)), k = k, p = 8*T__1^2*S__3^2/(beta[2]*(S__1^2*S__2^2-S__1^2*T__1^2-2*S__2^2*T__1^2+2*T__1^4)), w = -(1/2)*(S__1^4*S__2^4*beta[2]^2-S__1^4*S__2^2*T__1^2*beta[2]^2-2*S__1^2*S__2^4*T__1^2*beta[2]^2+2*S__1^2*S__2^2*T__1^4*beta[2]^2+16*S__3^2*T__1^2*k^2)/(beta[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] = -(1/2)*S__1^2*beta[2]-S__2^2*beta[2], beta[2] = beta[2]}]

> 

sol2[1]

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

Put values in

> 

params := {S__1 = 1, S__2 = -1, S__3 = 0, T__1 = 0}; A := eval(sol, params)

{S__1 = 1, S__2 = -1, S__3 = 0, T__1 = 0}

{chi = 0, k = k, p = 0, w = -(1/2)*beta[2], beta[1] = -(3/2)*beta[2], beta[2] = beta[2]}

Choose values for the free parameters

> 

free := {k = 1, beta[2] = 1}

{k = 1, beta[2] = 1}

> 

A11 := `union`(eval(A, free), free)

{1 = 1, chi = 0, k = 1, p = 0, w = -1/2, beta[1] = -3/2, beta[2] = 1}

Check solution to eqs2 - no division by zero but still not valid

> 

simplify(eval(eqs2, `union`(A11, params)))

{0, -2-2^(1/2), 1-2^(1/2)}

> 

NULL

 

NULL

Download find-p1.mw

adaptive stepsizes...

dharr 8697
Yesterday at 3:11 PM
0 0

@awass I think you are worried about unneccesary steps in the algorithm. However, both methods are efficient in this way.

For the output=Array method, dsolve knows where it needs answers and may use intermediate step between these as necessary, but will traverse the range only once. 

For the parameters method (and the usual operation with listprocedure or procedurelist), dsolve remembers internally where it is up to. So if you do ans(2.0) it will step out to 2.0 (with necessary intermediate steps). If you then ask for ans(2.1) it will start at 2.0 and step to 2.1; it does not redo 0 to 2.0. So in @acer's parameters method, Y contains the procedure to evaluate f(x) and it is successively applied to the x values in the Array A, each time just going the extra distance in x.

So both methods work essentially the same way.

The parameters method is more efficient and preferred when you have many parameters, since it does not duplicate the overhead of converting the differential equation to the numerical method each time. It is a bit more complicated to use.

Regarding your misconception about parameters being deprecated, I'm guessing you are thinking about that warning message when you don't use the parameters= option, but just leave a parameter unspecified. That warning message is saying don't leave them unspecified (deprecated), use the "new" parameters= option.

size of Array...

dharr 8697
November 12 2025
0 0

@awass "How does Maple know the size of the Array xpts?" I'm not sure I understand the question. If you mean how does dsolve or plot know, then intenally they just detect its size, e.g., by numelems(xpts), which you could also use. 

If you are asking about how I made xpts have a certain size, I made it from a list, in turn made by the sequence: seq(0..3.0, numelems=20) generates 20 points evenly spaced including the endpoints, so xpts has 20 entries.

Constants...

dharr 8697
November 09 2025
0 0

@tedh You say you left out your constant, so if we put it back you found (1/4)*x^4 + x^2 + C. Maple left out its constant and if we put it back Maple found (1/4)*x^4 + x^2 + 1 + W. So if the two constants are related by C = 1 + W, then the results are identical.

The value of the constant doesn't really matter, it will cancel out when a definite integral is involved. For example integrating from 0 to 1 from your result we get ((1/4)*1^4 + 1^2) - ((1/4)*0^4 + 0^2) = 5/4 - 0 = 5/4; from Maple's result we get ((1/4)*1^4 + 1^2 + 1) - ((1/4)*0^4 + 0^2 + 1) = 9/4 - 1 = 5/4.

Ahhh...

dharr 8697
November 08 2025
0 0

@acer OK. I've edited.

maybe...

dharr 8697
November 08 2025
0 0

@acer The OP wanted: "Then,  `-\\infty` to `- + \\infty`. "  Seemed like a typo, but I was in a literal mood this morning. Perhaps the OP can clarify.

solve for w...

dharr 8697
November 07 2025
0 0

@Andiguys Not following what tou are doing here, but you didn't include w in the variables to solve for, so maybe that is it?

native isolve...

dharr 8697
November 06 2025

@Alfred_F Maple's isolve does this directly and correctly.

Ellipse.mw

simplified cubic equation...

dharr 8697
November 06 2025
0 0

@Andiguys Solve for w separately, then Pn is found by solving a cubic polynomial, which can be given in a simple way.

Inside_root.mw

Inversion in Maple....

dharr 8697
November 05 2025
0 0

@dharr I realized that for the cases where Algfield gives one number as a polynomial of the other, inversion can be done by simply solving that polynomial. 

> 

restart;
> 

local gamma:

Example 1.

> 

a := RootOf(4*_Z^3 - 12*_Z^2 - 24*_Z - 25, index=1);
b := RootOf(2*_Z^3 - 3, index=1);
alias(alpha = a, beta = b);

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

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

alpha, beta

Algfield finds a in terms of b as a polynomial, so they a is in the same field defined by n

> 

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

eqa := %[1][];

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

alpha = beta^2+2*beta+1

So to invert, we just treat beta as an unknown x and solve
roots and evala(Roots()) work in the same way

> 

eq := alpha-subs(beta = x, rhs(eqa));
roots(eq);

alpha-x^2-2*x-1

[[(2/13)*alpha^2-(12/13)*alpha-28/13, 1], [-(2/13)*alpha^2+(12/13)*alpha+2/13, 1]]

Example 2

> 

c := RootOf(_Z^4 + 1, index=1);
d := RootOf(_Z^4 + 6*_Z^2 + 1, index=1);
alias(gamma = c, delta = d);

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

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

alpha, beta, gamma, delta

> 

evala(Algfield({c,d}));
> 

eqd := %[1][];
eq := delta-subs(gamma = x, rhs(eqd));
roots(eq);

[[delta = gamma^3-gamma^2+gamma], [], {gamma}, true, {}]

delta = gamma^3-gamma^2+gamma

delta-x^3+x^2-x

[[(1/4)*delta^3+(1/4)*delta^2+(7/4)*delta+3/4, 1]]

> 

NULL

Download evala2.mw

General procedure...

dharr 8697
November 04 2025

@Alfred_F Here is a general procedure, tested on your two examples. If it isn't clear or needs more features, please let me know. Choose more _Z1 values for more solutions, but they get very large very quickly.

Procedure QuadraticIsolve is in the startup code edit region.

It accepts a quadratic in two variables and returns integer solutions.
If there is no second argument it returns all solutions if there are a finite number of solutions

or general solutions parameterized by variable _Z1 if there are an infinite number of solutions.

If a set or range of integers is given as a second argument, then solutions corresponding only

to these values of _Z1 are given.
Use infolevel[QuadraticIsolve]:=2; to generate qualitative information about the solving process

and infolevel[QuadraticIsolve]:=3; for more quantitative information.

> 

restart
> 

infolevel[QuadraticIsolve] := 3

3

First case  -  there are 128 general solutions involving an arbitrary integer parameter _Z1

> 

P := 14*x^2-5*x*y-17*y^2+2*x+11*y-12; sols := QuadraticIsolve(P); nops(%)

14*x^2-5*x*y-17*y^2+2*x+11*y-12

QuadraticIsolve: disc = 977, transformed equation: X^2-977*Y^2 = -559328

QuadraticIsolve: hyperbola, solved by transformation

QuadraticIsolve: returning general solutions containing parameter _Z1

QuadraticIsolve: solutions may be non-integer for some values of the parameter

128

Look at one of the general solutions

> 

sols[92]

{x = (151471/1954)*977^(1/2)*((108832847723078562849+3481871275306470280*977^(1/2))^_Z1-(108832847723078562849-3481871275306470280*977^(1/2))^_Z1)+(4734529/1954)*(108832847723078562849+3481871275306470280*977^(1/2))^_Z1+(4734529/1954)*(108832847723078562849-3481871275306470280*977^(1/2))^_Z1-13/977, y = 318/977-(2524409/977)*(108832847723078562849+3481871275306470280*977^(1/2))^_Z1-(2524409/977)*(108832847723078562849-3481871275306470280*977^(1/2))^_Z1-(80763/977)*977^(1/2)*((108832847723078562849+3481871275306470280*977^(1/2))^_Z1-(108832847723078562849-3481871275306470280*977^(1/2))^_Z1)}

Choose a specific _Z1. (Some solutions may not be integer)

> 

evala(eval(sols[92], _Z1 = 1))

{x = 1054805055873880264681484, y = -1124825473059920639559012}

Specify chosen _Z1 values as the second argument (set or range); here 0 and 1. Only integer solutions are returned.

> 

sols := QuadraticIsolve(P, {0, 1}); nops(%)

QuadraticIsolve: disc = 977, transformed equation: X^2-977*Y^2 = -559328

QuadraticIsolve: hyperbola, solved by transformation

36

> 

sols

{{x = -4354274635632925848630, y = 4643321530916357188489}, {x = -842369127639866026312480, y = -650531927181484066810521}, {x = -86052394265244921080883408, y = -66455224963899021648732681}, {x = -1834717862750735136472687617, y = -1416887808357334226606681712}, {x = -1935541499837998630, y = -1494750342474165031}, {x = -10004971671472080, y = 10669125001482999}, {x = -97938928148608, y = 104440342385799}, {x = -4593556006657, y = 4898486956848}, {x = -710258662, y = -548507680}, {x = -225037758, y = 239976289}, {x = -10554777, y = 11255428}, {x = -103321, y = 110180}, {x = -4846, y = 5168}, {x = -1632, y = -1260}, {x = -16, y = -12}, {x = -1, y = 0}, {x = 20, y = -21}, {x = 3870, y = 2989}, {x = 395342, y = 305309}, {x = 8429063, y = 6509468}, {x = 54514467908, y = -58133265300}, {x = 172057228892, y = 132873721299}, {x = 3668421709343, y = 2832992527848}, {x = 374748386313759, y = 289404943697160}, {x = 7989987545117724, y = 6170385197335548}, {x = 23725301048449838, y = -25300241809328680}, {x = 2423663086517669214, y = -2584551489063081992}, {x = 51674773213364371719, y = -55105065055676535380}, {x = 163094654023171420679, y = 125952241260259919180}, {x = 3477331195842375008654, y = 2685420073047531974168}, {x = 355227494974617759355998, y = 274329648738604371091960}, {x = 1054805055873880264681484, y = -1124825473059920639559012}, {x = 22489434422549442376595999, y = -23982335477370792469719000}, {x = 2297412873087832090660316703, y = -2449920492316976033389932792}, {x = 48983000093817740901123542892, y = -52234605764925810264202912221}, {x = 154598945607990784875891494868, y = 119391305695605386194795009740}}

Check

> 

map2(eval, P, sols)

{0}

Second case

> 

QuadraticIsolve(3*x^2-5*x*y+7*y^2+2*x-11*y+19)

QuadraticIsolve: disc = -59, (a+c)*delta = 2100

QuadraticIsolve: imaginary ellipse; has no solutions

{}

 

``

Download DiophantineConics3.mw

syntax fixed...

dharr 8697
November 04 2025
0 0

@Andiguys Data needs to be in a set, not list if you use union. 

new_Q_solve_sand15.mw

remark...

dharr 8697
November 04 2025
0 0

@sursumCorda I took this idea from a much longer algorithm for finding when one field is a subfield of another, so I think it is much harder to properly distinguish the two cases. Glad you like it.

Using IntegerRelations:-LinearDependency...

dharr 8697
November 03 2025
0 0

@sursumCorda Here's an algorithm using IntegerRelations:-LinearDependency, based on the LLL algorithm. It uses floats as an intermediate step but checks the obtained result is exact. It does the original examples, but not their reverse (expressing degree 6 roots in terms of degree 2 or 3 roots).

> 

restart;
> 

local gamma:
> 

algconvert := proc(a::algnumext, b::algnumext)
description "writes algebraic number a as a polynomial in algebraic number b";
local Bs, i, Banum, cfs;
Bs := seq(b^i, i = 0..degree(op(1,b), ':-_Z')-1);
Banum:=evalf([Bs,a]);
cfs := IntegerRelations:-LinearDependency(Banum);
if evala(add(cfs*~[Bs,a])) <> 0 then return FAIL end if;
-evala(add(cfs[1..-2]*~[Bs]))/cfs[-1];
end proc:

Example 1.

> 

a := RootOf(4*_Z^3 - 12*_Z^2 - 24*_Z - 25, index=1);
b := RootOf(2*_Z^3 - 3, index=1);
alias(alpha = a, beta = b);

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

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

alpha, beta

> 

algconvert(a, b);
evala(a - %);
algconvert(b, a);
evala(b - %);

beta^2+2*beta+1

0

-(2/13)*alpha^2+(12/13)*alpha+2/13

0

Example 2

> 

c := RootOf(_Z^4 + 1, index=1);
d := RootOf(_Z^4 + 6*_Z^2 + 1, index=1);
alias(gamma = c, delta = d);

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

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

alpha, beta, gamma, delta

> 

algconvert(c, d);
evala(c - %);
algconvert(d, c);
evala(d - %);

(1/4)*delta^3+(1/4)*delta^2+(7/4)*delta+3/4

0

gamma^3-gamma^2+gamma

0

Example 3

> 

a1 := convert(3^(1/3), RootOf);
a2 := convert(4^(1/4), RootOf);
3^(1/3) + 2^(1/2) + 1^(1/1);
B := RootOf(evala(Minpoly(%, x)), index=2);
evala(%% - B);
alias(alpha__1 = a1, alpha__2 = a2, beta = B);
 

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

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

3^(1/3)+2^(1/2)+1

RootOf(_Z^6-6*_Z^5+9*_Z^4-2*_Z^3+9*_Z^2-60*_Z+50, index = 2)

0

alpha, beta, gamma, delta, alpha__1, alpha__2

> 

algconvert(a1, B);

FAIL

Succeeds with higher digits

> 

Digits := 50;

50

> 

algconvert(a1, B);
evala(a1 - %);
algconvert(a2, B);
evala(a2 - %);

277/151-(232/755)*beta-(174/755)*beta^2-(52/755)*beta^3+(213/755)*beta^4-(48/755)*beta^5

0

-428/151+(987/755)*beta+(174/755)*beta^2+(52/755)*beta^3-(213/755)*beta^4+(48/755)*beta^5

0

These fail

> 

algconvert(B, a1);
algconvert(B, a2);

FAIL

FAIL

> 

NULL

Download evala.mw

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