Dr. David Harrington

2830 Reputation

17 Badges

17 years, 273 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a professor of chemistry at the University of Victoria, BC, Canada, where 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

plot(F(R), R = 0 .. 100) is correct and works for me. You seemed to have some extra character in that line.


The discriminant is negative for arbitrary values of the a[i,j], so the roots are usually complex (edit: for real a[i,j]). But for values that make the disciminant zero, you will get two equal real roots. (@tomleslie left off a term of your expression).

z := k^2*a[1, 1]^2+k^2*b[1, 1]^2+4*k*a[0, 2]*a[1, 1]+4*k*b[0, 2]*b[1, 1]+4*a[0, 2]^2+4*b[0, 2]^2

k^2*a[1, 1]^2+k^2*b[1, 1]^2+4*k*a[0, 2]*a[1, 1]+4*k*b[0, 2]*b[1, 1]+4*a[0, 2]^2+4*b[0, 2]^2

ans := [solve(z, k)];

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

The discriminant "b^2-4*a*c" can be factored, so the square root disappears. But it is negative, so the roots are complex, unless the discriminant becomes zero

discrim(z, k);

-16*a[0, 2]^2*b[1, 1]^2+32*a[0, 2]*a[1, 1]*b[0, 2]*b[1, 1]-16*a[1, 1]^2*b[0, 2]^2

-16*(a[0, 2]*b[1, 1]-a[1, 1]*b[0, 2])^2

realvals := {a[0, 2] = 1, a[1, 1] = 1, b[0, 2] = 2, b[1, 1] = 2}

{a[0, 2] = 1, a[1, 1] = 1, b[0, 2] = 2, b[1, 1] = 2}

eval(discrim(z, k), realvals);


eval(ans, realvals);

[-2, -2]



Download discrim.mw

The syntax is solve({eq1,eq2,...},{var1,var2,...}). You don't seem to have collected the equations and variables in sets (or lists). Use fsolve if you want a numerical solution. In general, uploading your worksheet (green up-arrow) is helpful for us to diagnose your problem.

If dat is your matrix, then plot(dat) will work. You can add various plot options for label fonts etc.



Need specific N value to make progress.

N := 5


eq := add(alpha[n]*kappa^n*n, n = 0 .. N)^2+add(beta[n]*kappa^n*n, n = 0 .. N)^2


Symbolically can find the zero root, but will need specific values of alpha and beta to get more

solve(eq, kappa);

0, 0, RootOf((25*alpha[5]^2+25*beta[5]^2)*_Z^8+(40*alpha[4]*alpha[5]+40*beta[4]*beta[5])*_Z^7+(30*alpha[3]*alpha[5]+16*alpha[4]^2+30*beta[3]*beta[5]+16*beta[4]^2)*_Z^6+(20*alpha[2]*alpha[5]+24*alpha[3]*alpha[4]+20*beta[2]*beta[5]+24*beta[3]*beta[4])*_Z^5+(10*alpha[1]*alpha[5]+16*alpha[2]*alpha[4]+9*alpha[3]^2+10*beta[1]*beta[5]+16*beta[2]*beta[4]+9*beta[3]^2)*_Z^4+(8*alpha[1]*alpha[4]+12*alpha[2]*alpha[3]+8*beta[1]*beta[4]+12*beta[2]*beta[3])*_Z^3+(6*alpha[1]*alpha[3]+4*alpha[2]^2+6*beta[1]*beta[3]+4*beta[2]^2)*_Z^2+(4*alpha[1]*alpha[2]+4*beta[1]*beta[2])*_Z+alpha[1]^2+beta[1]^2)

Choose some specific values

vals := zip(`=`, [seq(alpha[n], n = 0 .. N), seq(beta[n], n = 0 .. N)], [seq((rand(-3 .. 3))(), i = 1 .. 2*N+2)])

[alpha[0] = -1, alpha[1] = 2, alpha[2] = 2, alpha[3] = 0, alpha[4] = 2, alpha[5] = 2, beta[0] = 2, beta[1] = -3, beta[2] = 0, beta[3] = -2, beta[4] = -3, beta[5] = 1]

eq2 := eval(eq, vals)


ans := solve(eq2, kappa);

0, 0, RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 1), RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 2), RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 3), RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 4), RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 5), RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 6), RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 7), RootOf(125*_Z^8+40*_Z^7+148*_Z^6+224*_Z^5+110*_Z^4+104*_Z^3+52*_Z^2+16*_Z+13, index = 8)


[0., 0., .2359930256+.5296052631*I, .4481955872+1.141969704*I, -.1028659271+.6003394077*I, -.7413226857+0.6729615150e-1*I, -.7413226857-0.6729615150e-1*I, -.1028659271-.6003394077*I, .4481955872-1.141969704*I, .2359930256-.5296052631*I]

fsolve(eq2, complex);

-.741322685659808-0.672961515020724e-1*I, -.741322685659808+0.672961515020724e-1*I, -.102865927127570-.600339407715825*I, -.102865927127570+.600339407715825*I, 0., 0., .235993025621495-.529605263094995*I, .235993025621495+.529605263094995*I, .448195587165883-1.14196970387710*I, .448195587165883+1.14196970387710*I


Download solve.mw

Correct syntax is GroupOrder(ds[2]);

I agree; it's a bug. The coefficients are the same since simplify(c1(n)-c2(n)) assuming n::integer; returns zero. The fhat1 sum isn't evaluated correctly, but if sum is replaced by Sum, then (at least in Maple 2015), evalf(eval(fhat1,x=0.9)); returns 0.451026812 as expected.

Order is important, so use lists, not sets


[[a[0], b[0], c[0]], [a[1], b[1], c[1]], [a[2], b[2], c[2]], [a[3], b[3], c[3]], [a[4], b[4], c[4]]]


[i, j, k]


[i*a[0]+j*b[0]+k*c[0], i*a[1]+j*b[1]+k*c[1], i*a[2]+j*b[2]+k*c[2], i*a[3]+j*b[3]+k*c[3], i*a[4]+j*b[4]+k*c[4]]


Download vecs.mw

The output of dsolve is an equation. You only want to plot the right-hand side. There may be other solutions.


ODE := u(x)*(diff(u(x), x))-(1/10)*(diff(u(x), `$`(x, 2))) = 0;

u(x)*(diff(u(x), x))-(1/10)*(diff(diff(u(x), x), x)) = 0

Bcs := u(0) = 1, u(1) = 0;

u(0) = 1, u(1) = 0

ans := rhs(dsolve({Bcs, ODE}));


plot(ans, x = 0 .. 1);




Download solve_ode.mw

I think this is what you want, but perhaps you want a progressive animation, rather than a single point moving? Or something else altogether...



xcoords:=[$1..100]: #x coordinates
ycoords:=[seq(i^2,i=1..nops(xcoords))]:   #y coordinates




Download animate.mw

The error message is because you are assigning to kt_poly, which is not allowed, because it is passed into the procedure (and already has a value). The easiest fix is just to have the procedure return the value of Quotient as the last line of the procedure. Now you get an error about Quotient not allowing multivariate polynomials, which is because the X passed into the procedure in cyclo_poly[i] is diiferent from the local X. The simplest fix is to pass the variable you want (X) as a third argument to the procedure.


Programm zur Brechnung der Kreisteilungspolynome_2022-04-05 Ki



restart; with(Algebraic); with(NumberTheory)

Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received NumberTheory

n := 6;



cyclo_poly := Vector[row](1 .. n, 0)

cyclo_poly := Vector[row](6, {(1) = X-1, (2) = X+1, (3) = 0, (4) = 0, (5) = 0, (6) = 0})

i := 1;



c_poly := proc (i, kt_poly, X) local j, hz1, hz2; j := i+1; hz1 := X^j-1; hz2 := kt_poly; Algebraic:-Quotient(hz1, hz2, X) end proc:

cyclo_poly[2] := c_poly(i, cyclo_poly[i], X);


i := 1;



j := i+1;





Download Programm_zur_Brechnung_der_Kreisteilungspolynome_2022-04-05_Ki.mw

s := series(1/tan(x), x, 9);


s2 := convert(s, polynom)+`...`



Download dots.mw

A plot of the indefinite integral shows discontinuities at 0, Pi, 2*Pi. If you take the limits from the left or right as appropriate, and avoid the discontinuities, you can get the right answer.

restart; with(Student[Calculus1])

Ellipse := [3*cos(theta), 2*sin(theta)]:


Why do the answers from evaluating EllALintegral differ from those provided by the ArcLength command?

evalf(ArcLength(Ellipse, theta = 0 .. 2*Pi));


evalf(Student:-Calculus1:-ArcLength(Ellipse, theta = 0 .. (1/4)*Pi));


Student:-Calculus1:-ArcLength(Ellipse, theta = 0 .. 2*Pi, output = integral);

Int((9*sin(theta)^2+4*cos(theta)^2)^(1/2), theta = 0 .. 2*Pi)


plot(sqrt(9*sin(theta)^2+4*cos(theta)^2), theta = 0 .. 2*Pi);

EllALintegral := int(sqrt(9*sin(theta)^2+4*cos(theta)^2), theta);

-3*((5*cos(theta)^2-9)*(cos(theta)^2-1))^(1/2)*(1-cos(theta)^2)^(1/2)*EllipticE(cos(theta), (1/3)*5^(1/2))/((5*cos(theta)^4-14*cos(theta)^2+9)^(1/2)*sin(theta))

plot(EllALintegral, theta = 0 .. 2*Pi);

limit(EllALintegral, theta = 2*Pi, left);





2*(limit(EllALintegral, theta = Pi, left)-(limit(EllALintegral, theta = 0, right)));




Download ArcLength_of_Ellipse.mw

I used eval to substitute C1 and C2 into psi, which had square brackets around that I removed. Then the limit works. You should be usung LinearAlgebra rather than linalg.


I usually force the hardware type I want, here integer[8]. Then you get an explicit error about the overflow.

Note that tail recursion gives a massive speedup - as the documentation says: "During compilation, tail-recursive calls are detected and transformed into iteration. This means that tail-recursive procedures compile to very fast code."

Download fibonacci.mw

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