Advanced Search

dharr

Dr. David Harrington

8672 Reputation

22 Badges

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

ThueSolve...

dharr 8672
January 05 2025
2 0

Well, as is often the case, Maple needs help to go the extra mile. Actually, ThueSolve sends this case to isolve, but without the inequalities.

> 

restart;
> 

isolve({a*b=4,a>=1,b>=1});

Warning, solutions may have been lost

> 

sols:=NumberTheory:-ThueSolve(a*b=4);
goodsols:=select(x->eval((a>=1 and b>=1),x),sols);

[[a = -4, b = -1], [a = -2, b = -2], [a = -1, b = -4], [a = 1, b = 4], [a = 2, b = 2], [a = 4, b = 1]]

[[a = 1, b = 4], [a = 2, b = 2], [a = 4, b = 1]]

> 

isolve(a*b=4);

{a = -4, b = -1}, {a = -2, b = -2}, {a = -1, b = -4}, {a = 1, b = 4}, {a = 2, b = 2}, {a = 4, b = 1}

> 

NULL

Download Thue.mw

No solution by elimination of x,y,z...

dharr 8672
December 31 2024
1 4

My take on it is that the elimination process does not lead to any valid solutions.

> 

restart;
> 

expr :=[a, b, c]^~2 =~ 4*[y*z/((x + y)*(x + z)), z*x/((y + z)*(y + x)), x*y/((z + x)*(z + y))];

[a^2 = 4*y*z/((x+y)*(x+z)), b^2 = 4*z*x/((y+z)*(x+y)), c^2 = 4*x*y/((x+z)*(y+z))]

The following assumes denominators are not zero

> 

expr2:=map(w->numer(normal(lhs(w)-rhs(w))),expr);

[a^2*x^2+a^2*x*y+a^2*x*z+a^2*y*z-4*y*z, b^2*x*y+b^2*x*z+b^2*y^2+b^2*y*z-4*x*z, c^2*x*y+c^2*x*z+c^2*y*z+c^2*z^2-4*x*y]

> 

denoms1:={map(denom@rhs,expr)[]};

{(x+y)*(x+z), (x+z)*(y+z), (y+z)*(x+y)}

Interesting result - probably an error message or null result would be better

> 

eliminate(expr2,{x,y,z});

[{x = 0, y = FAIL, z = 0}, {}]

Let's do it "by hand" - eliminate x between pairs, leaving 2 eqns in y and z. Denominators in expressions for x must be nonzero

> 

e12:=eliminate(expr2[1..2],x);
e23:=eliminate(expr2[2..3],x);
expr3 := e12[2] union e23[2];

[{x = -b^2*y*(y+z)/(b^2*y+b^2*z-4*z)}, {4*y*z*(a^2*b^2*y^2-a^2*b^2*z^2-b^4*y^2-2*b^4*y*z-b^4*z^2+4*a^2*z^2+8*b^2*y*z+8*b^2*z^2-16*z^2)}]

[{x = -c^2*z*(y+z)/(c^2*y+c^2*z-4*y)}, {b^2*c^2*y^3+b^2*c^2*y^2*z-b^2*c^2*y*z^2-b^2*c^2*z^3-4*b^2*y^3-4*b^2*y^2*z+4*c^2*y*z^2+4*c^2*z^3}]

{4*y*z*(a^2*b^2*y^2-a^2*b^2*z^2-b^4*y^2-2*b^4*y*z-b^4*z^2+4*a^2*z^2+8*b^2*y*z+8*b^2*z^2-16*z^2), b^2*c^2*y^3+b^2*c^2*y^2*z-b^2*c^2*y*z^2-b^2*c^2*z^3-4*b^2*y^3-4*b^2*y^2*z+4*c^2*y*z^2+4*c^2*z^3}

> 

denoms2:={denom(rhs(e12[1][])), denom(rhs(e23[1][]))};

{b^2*y+b^2*z-4*z, c^2*y+c^2*z-4*y}

Now eliminate y between the remaining 2 equations - denom has no restrictions on z

> 

e13:=eliminate(expr3,y);

[{y = -z*c^2*(-a^2*b^4*c^2+a^4*b^2+2*a^2*b^4+5*a^2*b^2*c^2+b^6+3*b^4*c^2-4*a^4-12*a^2*b^2-4*a^2*c^2-16*b^4-12*b^2*c^2+16*a^2+48*b^2)/(-a^2*b^4*c^4+a^4*b^2*c^2+2*a^2*b^4*c^2+3*a^2*b^2*c^4+b^6*c^2+3*b^4*c^4-16*b^4*c^2-4*b^2*c^4-16*a^4-16*a^2*b^2-32*a^2*c^2-16*c^4+128*a^2+64*b^2+128*c^2-256)}, {c*z*(b-2)*(b+2)*(a-2)*(a+2)*(-a*b*c+a^2+b^2+c^2-4)*(a*b*c+a^2+b^2+c^2-4)}]

> 

denoms3:=denom(rhs(e13[1][]));

-a^2*b^4*c^4+a^4*b^2*c^2+2*a^2*b^4*c^2+3*a^2*b^2*c^4+b^6*c^2+3*b^4*c^4-16*b^4*c^2-4*b^2*c^4-16*a^4-16*a^2*b^2-32*a^2*c^2-16*c^4+128*a^2+64*b^2+128*c^2-256

> 

eqz:={isolate(e13[2][],z)};

{z = 0}

From which we conclude y=0

> 

eqy:=eval(e13[1],eqz);

{y = 0}

So we have a problem at the eliminate x stage, and with the original equations

> 

eval(denoms2, eqy union eqz);
eval(denoms1, eqy union eqz);

{0}

{0, x^2}

> 

solve(expr,{x,y,z},parametric); # correct

[]

Here the c<>0 solution doesn't solve the original equations (zero denoms). And solve is too lazy to do the c=0 case (default option - see ?solve,parametric

> 

ans := solve(expr2,{x,y,z},parametric);

ans := piecewise(c = 0, %SolveTools[Parametric]({-4*y*x, a^2*x^2+a^2*x*y+a^2*x*z+a^2*y*z-4*y*z, b^2*x*y+b^2*x*z+b^2*y^2+b^2*y*z-4*x*z}, {x, y, z}, {a, b}), c <> 0, [{x = 0, y = 0, z = 0}])

Force it. None of these solve the original equations either.

> 

op(2,ans);
value(%);

%SolveTools[Parametric]({-4*y*x, a^2*x^2+a^2*x*y+a^2*x*z+a^2*y*z-4*y*z, b^2*x*y+b^2*x*z+b^2*y^2+b^2*y*z-4*x*z}, {x, y, z}, {a, b})

piecewise(b = 0, [{x = 0, y = 0, z = z}, {x = 0, y = y, z = 0}], b <> 0, [{x = 0, y = 0, z = z}])

In one step

> 

ans := solve(expr2,{x,y,z},parametric=full);

ans := piecewise(c = 0, piecewise(b = 0, [{x = 0, y = 0, z = z}, {x = 0, y = y, z = 0}], b <> 0, [{x = 0, y = 0, z = z}]), c <> 0, [{x = 0, y = 0, z = 0}])

Same result.

> 

SolveTools:-PolynomialSystem(expr2,[x,y,z]);

{x = 0, y = 0, z = 0}

> 

 

> 

NULL

Download eliminate.mw

remove restriction...

dharr 8672
December 29 2024
1 2

Your question is not clear. Is this what you want to do?

> 

restart

Output from solve - the variables we solved for are only on the left of the =.

> 

case1 := {alpha = alpha, beta = gamma, delta = delta, gamma = gamma, k = k, lambda = 0, m = 2*n, mu = mu, n = n, sigma = 32*alpha*mu^2*n^4/a[-1]^2, w = -2*alpha*k^2*n-4*alpha*mu^2*n+delta^2, a[-1] = a[-1], a[0] = 0, a[1] = 0}

To remove the restriction on beta so it is now arbitrary (actually would remove if beta was on the right also)

> 

case2 := remove(has, case1, beta)

{alpha = alpha, delta = delta, gamma = gamma, k = k, lambda = 0, m = 2*n, mu = mu, n = n, sigma = 32*alpha*mu^2*n^4/a[-1]^2, w = -2*alpha*k^2*n-4*alpha*mu^2*n+delta^2, a[-1] = a[-1], a[0] = 0, a[1] = 0}

A little more specific for the left-hand side.

> 

case3 := remove(proc (q) options operator, arrow; has(lhs(q), beta) end proc, case1)

{alpha = alpha, delta = delta, gamma = gamma, k = k, lambda = 0, m = 2*n, mu = mu, n = n, sigma = 32*alpha*mu^2*n^4/a[-1]^2, w = -2*alpha*k^2*n-4*alpha*mu^2*n+delta^2, a[-1] = a[-1], a[0] = 0, a[1] = 0}

> 

NULL

Download removal.mw

Array...

dharr 8672
December 27 2024
2 2

Put the plots in an Array to achieve this. Attached is an example. (grid is for the fineness of the points plotted, and insequence is, as you found, for animations.)

grid.mw

slightly shorter...

dharr 8672
December 23 2024
1 0

I don't know a way that doesn't involve a string. But this is slightly shorter:

cat(``,String(x)[1]);

solve...

dharr 8672
December 21 2024
2 1

I think you just want to solve the 3 equations simultaneously - if there is a solution it lies on the line; if no solution it does not

solve({seq(l[i] = P[i], i = 1 .. 3)}, alpha)

2024-12-21_Q_3D_point_lies_on_3D_line.mw

complex conjugate...

dharr 8672
December 21 2024
2 0

It means complex conjugate, i.e., every i is replaced by -i.

Hyperbolas...

dharr 8672
December 19 2024
1 2

No artificial intelligence used.

We consider an ellipse x^2/a^2+y^2/b^2-1=0 and 2 vertices of this ellipse A(a,0) and B(0,b). We imagine a variable equilateral hyperbola passing through the points O, A and B. This curve meets the ellipse at 2 other points A1 and B1. Show that the line A1B1 passes through a fixed point.

https://www.mapleprimes.com/questions/239553-Variable-Equilateral-Hyperbole-Passing

> 

restart
> 

with(Student:-MultivariateCalculus); with(plots); _local(D)

Define ellipse and an arbitrary hyperbola (B^2-4A*C>0) (see Wikipedia)

> 

ellipse := x^2/a^2+y^2/b^2-1; hyperbola := A*x^2+B*x*y+C*y^2+D*x+E*y+F

x^2/a^2+y^2/b^2-1

A*x^2+B*x*y+C*y^2+D*x+E*y+F

The hyperbola is equilateral when A=-C

> 

hyperbola2 := eval(hyperbola, C = -A)

A*x^2-A*y^2+B*x*y+D*x+E*y+F

Intersection points between ellipse and hyperbola are at O (assume this means the origin), A, B. Use these to eliminate D,E,F

> 

eqO := eval(hyperbola2, {x = 0, y = 0}); eqA := eval(hyperbola2, {x = a, y = 0}); eqB := eval(hyperbola2, {x = 0, y = b})

F

A*a^2+D*a+F

-A*b^2+E*b+F

> 

solve({eqA, eqB, eqO}, {D, E, F}); hyperbola3 := eval(hyperbola2, %)

{D = -A*a, E = A*b, F = 0}

-A*a*x+A*b*y+A*x^2-A*y^2+B*x*y

A1 and B1 are two points on both the ellipse and the hyperbola. The 4 solutions are points B, A, A1, B1; A1 and B1 are complicated expressions in A.B,a,b

> 

ans := sort([solve({ellipse, hyperbola3}, {x, y}, explicit)]); nops(%)

4

> 

A1 := eval([x, y], ans[3]); B1 := eval([x, y], ans[4])

If the assertion is true, then the line A1-B1 intersects at the same point for different {A,B} pairs.

> 

L1 := Line(eval([A1, B1][], {A = A__1, B = B__1})); L2 := Line(eval([A1, B1][], {A = A__2, B = B__2}))

And we indeed find a value independent of A__1, A__2, B__1, B__2

> 

X__0 := GetIntersection(L1, L2)

[-b^2*a/(a^2+b^2), -a^2*b/(a^2+b^2)]

Construct a diagram for two different hyperbolas

> 

params0 := {a = 1.5, b = 1.}; params := `union`(params0, {A = 1., B = 1.}); params2 := `union`(remove(has, params, B), {B = 3.})

{A = 1., B = 1., a = 1.5, b = 1.}

{A = 1., B = 3., a = 1.5, b = 1.}

Red points are O,A,B,X; Black and blue points are A1,B1 for the blue and black hyperbolas

> 

ellipseplot := implicitplot(eval(ellipse, params), x = -1.5 .. 1.5, y = -1 .. 1, color = magenta); hyperbola1plot := implicitplot(eval(hyperbola3, params), x = -2 .. 2, y = -1.5 .. 1.5, color = black); hyperbola2plot := implicitplot(eval(hyperbola3, params2), x = -2 .. 2, y = -1.5 .. 1.5, color = blue); pts := pointplot(eval([[0, 0], [a, 0], [0, b], X__0], params), symbol = solidcircle, color = red, symbolsize = 10); pts1 := pointplot(eval([A1, B1], params), symbol = solidcircle, color = black, symbolsize = 10); pts2 := pointplot(eval([A1, B1], params2), symbol = solidcircle, color = blue, symbolsize = 10); pline1 := plottools:-line(eval([A1, B1][], params), color = black, linestyle = dash); pline2 := plottools:-line(eval([A1, B1][], params2), color = blue, linestyle = dash); t1 := textplot(eval([[A1[], A__1], [B1[], B__1]], params), color = black, align = below); t2 := textplot(eval([[A1[], A__1], [B1[], B__1]], params2), color = blue, align = right); t3 := textplot(eval([[0, 0, O], [X__0[], X], [a, 0, A]], params0), color = red, align = right); t4 := textplot(eval([[0, b, B]], params0), color = red, align = above); display(t1, t2, t3, t4, pline1, pline2, ellipseplot, hyperbola1plot, hyperbola2plot, pts, pts1, pts2, axes = none, scaling = constrained)

> 

NULL

Download Ellipse.mw

sort...

dharr 8672
December 17 2024
1 3

I'm not sure exactly what you mean by labelled sets. You can do something like this

ans := sort([solve({x^2 - x + 1}, x)])

sort will sort in a consistent way from session to session, though what exact the order is is nontrivial.

missing multiplication...

dharr 8672
December 17 2024
1 0

You are missing the multiplication between the (x+7) and (x-1), so l45 := (x + 7)*(x - 1) = (1 + x)^2;

You can use a space instead of * if you are using 2D input

 

Octagon...

dharr 8672
December 17 2024
4 0
> 

restart
> 

_local(O); with(Student:-MultivariateCalculus)

Vectors of adjacent sides

> 

A := [a__1, a__2]; B := [b__1, b__2]

Origin and far point C

> 

O := [0, 0]; C := A+B

Midpoints

> 

OA := (O+A)*(1/2); OB := (O+B)*(1/2); AC := (A+C)*(1/2); BC := (B+C)*(1/2)

Area of polygon formula. Accepts a sequence of points in cyclic order.

> 

Area := proc () local n, V, i; n := nargs; V := Array(0 .. n-1, [args]); simplify((1/2)*add(-V[`mod`(i+1, n)][1]*V[i][2]+V[`mod`(i+1, n)][2]*V[i][1], i = 0 .. n-1)) end proc

Parallelogram area

> 

p_area := Area(O, B, C, A)

-a__1*b__2+a__2*b__1

Lines forming octagon

> 

lines := [Line(O, BC), Line(OB, C), Line(B, AC), Line(BC, A), Line(C, OA), Line(AC, O), Line(A, OB), Line(OA, B), Line(O, BC)]
> 

intersections := seq(GetIntersection(lines[i], lines[i+1]), i = 1 .. 8)
> 

o_area := Area(intersections)

-(1/6)*b__2*a__1+(1/6)*b__1*a__2

> 

simplify(o_area/p_area)

1/6

> 

NULL

Download Parallelogram.mw

Edit: With plot: Parallelogram2.mw

minor issues...

dharr 8672
December 16 2024
1 0

See changes in red

> 

test := simplify((int((x^2+x)*(x-a+I*b)/abs(x-a-I*b)^2, x = -5 .. 5))/(I*Pi))

(1/2)*((-I*a^2+(-I+2*b)*a+I*b^2+b)*ln(a^2+b^2-10*a+25)+(I*a^2+(I-2*b)*a-I*b^2-b)*ln(a^2+b^2+10*a+25)+(-2*a^2+(-(4*I)*b-2)*a-2*b*(I-b))*arctan((a-5)/b)+(2*a^2+((4*I)*b+2)*a+2*b*(I-b))*arctan((a+5)/b)-(20*I)*a-20*I+20*b)/Pi

Above: the "full" definite integral.  The integral is being taken with respect to x

test1 := simplify(int((x^2+x)*(x-a+I*b)/abs(x-a+I*b)^2, x))NULL

(1/2)*(a^2+(1+(2*I)*b)*a+b*(I-b))*ln(a^2-2*a*x+b^2+x^2)+arctan((-x+a)/b)*(-I*a^2+(-I+2*b)*a+I*b^2+b)+x*(I*b+a+(1/2)*x+1)

The indefinite integralNULL

> 

test2 := simplify(subs(x = 5, test1)-subs(x = -5, test1))

10+(1/2)*(a^2+(1+(2*I)*b)*a+b*(I-b))*ln(a^2+b^2-10*a+25)-(1/2)*(a^2+(1+(2*I)*b)*a+b*(I-b))*ln(a^2+b^2+10*a+25)+arctan((a-5)/b)*(-I*a^2+(-I+2*b)*a+I*b^2+b)+(I*a^2+(I-2*b)*a-I*b^2-b)*arctan((a+5)/b)+(10*I)*b+10*a

"Hard coding" the definite integral

> 

simplify(I*test*Pi-test2)

0

 

The above result is now zero.

Download Complex_valued_integrals2.mw

ListTools:-Classify...

dharr 8672
December 15 2024
1 3

Maybe you only need to classify, and not order the sublists after classification? If so you only need LL and not the sort line in the procedure. I'm not too familiar with these orderings and I'm not sure which one you actually want.

> 

restart;
> 

F := [x^2 - y - 1, y^2 - y, y^3 - y, x^2 - 4, y^3 - 1, z^2 - 4, 2*x^2, x^2*y];

[x^2-y-1, y^2-y, y^3-y, x^2-4, y^3-1, z^2-4, 2*x^2, x^2*y]

> 

monomial_classify := proc(L::list(polynom), ord)
  local LL, Lmon, xx;
  Lmon := xx -> Groebner:-LeadingMonomial(xx, ord);
  LL := map(convert,convert(ListTools:-Classify(Lmon, L),'list'),'list');
  sort(LL, (a, b) -> Groebner:-TestOrder(Lmon(a[1]), Lmon(b[1]), ord));
end proc:
> 

monomial_classify(F, grlex(z,y,x));

[[2*x^2, x^2-4, x^2-y-1], [y^2-y], [z^2-4], [x^2*y], [y^3-1, y^3-y]]

> 

 

Download sort2.mw

Explore...

dharr 8672
December 14 2024
0 0

I'm assuming you want two plots, of Cb(t) and Cm(t) vs time, and not Cb(t) vs Cm(t). You seem to have cut and pasted things, so I set up to solve the problem from live commands. You may have to adjust some things to get exactly what you want. (The plot doesn't show here in Mapleprimes, but does if you download the document.) I had to make a guess about the meaning of c0 and c; they may be the wrong way around.

restart

 

 

 

 

Absorption from stomach

de1 := diff(Cm(t), t) = -Ka*Cm(t)                            

 

Glucose concentration in blood

de2 := diff(Cb(t), t) = Ka*Cm(t)-Ku*Cb(t)

 

Initial conditions

 

ins := Cb(0) = c0, Cm(0) = c 

 

Solve

solution := dsolve({de1, de2, ins})

{Cb(t) = -(-Ka+Ku)*(Ka*c+Ka*c0-Ku*c0)*exp(-Ku*t)/(Ka-Ku)^2-Ka*exp(-Ka*t)*c/(Ka-Ku), Cm(t) = c*exp(-Ka*t)}

Ceqns := eval([Cb(t), Cm(t)], solution)

[-(-Ka+Ku)*(Ka*c+Ka*c0-Ku*c0)*exp(-Ku*t)/(Ka-Ku)^2-Ka*exp(-Ka*t)*c/(Ka-Ku), c*exp(-Ka*t)]

Explore(plot(Ceqns, t = 0 .. 5, 0 .. 30), Ka = 1 .. 10, Ku = 1 .. 10, c0 = 0 .. 20, c = 0 .. 20, initialvalues = [Ka = 2, Ku = 1, c0 = 10, c = 10])

 

NULL

NULL

``

NULL

NULL

Download Explore_plot.mw

sfloat-hfloat conversion...

dharr 8672
December 11 2024
3 0

See the ?evalhf help page, There is overhead converting from software floats to hardware floats and back again. The conversion to hfloats is presumably still present below, but it doesn't have to convert back.

> 

restart;
> 

UseHardwareFloats:=false;

false

> 

CodeTools:-Usage( for i from 1 to 100 by 0.001 do exp(i) end do):
CodeTools:-Usage( for i from 1 to 100 by 0.001 do evalhf[hfloat](exp(i)) end do):

memory used=0.87GiB, alloc change=109.00MiB, cpu time=2.31s, real time=2.42s, gc time=93.75ms
memory used=12.84MiB, alloc change=0 bytes, cpu time=78.00ms, real time=68.00ms, gc time=0ns

 

 

Download evalhf.mw

First 11 12 13 14 15 16 17 Last Page 13 of 84