dharr

Dr. David Harrington

8195 Reputation

22 Badges

20 years, 331 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 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

There were errors in V. Now it works but is very slow, even with only a 10 x 10 grid. There are common expressions in the integrand that some hand coding could optimize. One could also play with some of the integration options (method, errors).

restart; with(plots); with(LinearAlgebra); with(Student[VectorCalculus]); V := proc (k, x, t) options operator, arrow; -((1/2)*I)*exp(I*k*x-k^2*t)*(1/(k-I)+1/(k+I)-k*(1/(k^2+I)+1/(k^2-I)))/Pi end proc

proc (k, x, t) options operator, arrow; Student:-VectorCalculus:-`-`(Student:-VectorCalculus:-`*`(Student:-VectorCalculus:-`*`(Student:-VectorCalculus:-`*`(I, Student:-VectorCalculus:-`*`(2, Pi)^Student:-VectorCalculus:-`-`(1)), exp(Student:-VectorCalculus:-`+`(Student:-VectorCalculus:-`*`(Student:-VectorCalculus:-`*`(I, k), x), Student:-VectorCalculus:-`-`(Student:-VectorCalculus:-`*`(k^2, t))))), Student:-VectorCalculus:-`+`(Student:-VectorCalculus:-`+`(Student:-VectorCalculus:-`+`(k, Student:-VectorCalculus:-`-`(I))^Student:-VectorCalculus:-`-`(1), Student:-VectorCalculus:-`+`(k, I)^Student:-VectorCalculus:-`-`(1)), Student:-VectorCalculus:-`-`(Student:-VectorCalculus:-`*`(k, Student:-VectorCalculus:-`+`(Student:-VectorCalculus:-`+`(k^2, I)^Student:-VectorCalculus:-`-`(1), Student:-VectorCalculus:-`+`(k^2, Student:-VectorCalculus:-`-`(I))^Student:-VectorCalculus:-`-`(1))))))) end proc

V(k, x, t)

-((1/2)*I)*exp(I*k*x-k^2*t)*(1/(k-I)+1/(k+I)-k*(1/(k^2+I)+1/(k^2-I)))/Pi

phi1 := (1/8)*Pi; phi2 := 7*Pi*(1/8); k1 := proc (r) options operator, arrow; r*exp(I*phi1) end proc; k2 := proc (r) options operator, arrow; r*exp(I*phi2) end proc; dk1 := proc (r) options operator, arrow; diff(k1(r), r) end proc; dk2 := proc (r) options operator, arrow; diff(k2(r), r) end proc

Integrand now real

integr := `assuming`([simplify(evalc(V(k1(r), x, t)*dk1(r)-V(k2(r), x, t)*dk2(r)))], [r > 0])

r*exp(-(1/2)*(r*t*2^(1/2)+2*cos((3/8)*Pi)*x)*r)*(((r^8-2*r^4+1)*2^(1/2)-2*r^6-2*r^2)*cos((1/2)*(r*t*2^(1/2)-2*cos((1/8)*Pi)*x)*r)+sin((1/2)*(r*t*2^(1/2)-2*cos((1/8)*Pi)*x)*r)*(r^8*2^(1/2)+2*r^6-2*r^2-2^(1/2)))/((r^8+1)*(1+r^4+r^2*2^(1/2))*Pi)

evalf(eval(integr, {r = 2, t = 1, x = 1}))

0.1505152135e-2

u1 := unapply(Int(integr, r = 0 .. infinity), x, t)

proc (x, t) options operator, arrow; Int(r*exp(-(1/2)*(r*t*2^(1/2)+2*cos((3/8)*Pi)*x)*r)*(((r^8-2*r^4+1)*2^(1/2)-2*r^6-2*r^2)*cos((1/2)*(r*t*2^(1/2)-2*cos((1/8)*Pi)*x)*r)+sin((1/2)*(r*t*2^(1/2)-2*cos((1/8)*Pi)*x)*r)*(r^8*2^(1/2)+2*r^6-2*r^2-2^(1/2)))/((r^8+1)*(1+r^4+r^2*2^(1/2))*Pi), r = 0 .. infinity) end proc

Now real - but slow

evalf(u1(2, 1))

0.5826332317e-1

u := proc (x, t) options operator, arrow; exp(-x/sqrt(2))*cos(t-x/sqrt(2))+u1(x, t) end proc

proc (x, t) options operator, arrow; Student:-VectorCalculus:-`+`(Student:-VectorCalculus:-`*`(exp(Student:-VectorCalculus:-`-`(Student:-VectorCalculus:-`*`(x, sqrt(2)^Student:-VectorCalculus:-`-`(1)))), cos(Student:-VectorCalculus:-`+`(t, Student:-VectorCalculus:-`-`(Student:-VectorCalculus:-`*`(x, sqrt(2)^Student:-VectorCalculus:-`-`(1)))))), u1(x, t)) end proc

Very slow but it works

CodeTools:-Usage(plot3d(u(x, t), x = 0 .. 3, t = 0 .. 2*Pi, grid = [10, 10], adaptmesh = false, axes = boxed, shading = zhue, style = surface, labels = ["x", "t", "u(x,t)"]))

memory used=14.20GiB, alloc change=32.00MiB, cpu time=93.38s, real time=80.66s, gc time=22.44s

NULL

Download fokas_method.mw

 

Here's an attempt at a ternary plot. If you want tick marks along the edges of the triangle, that would be tricky to do, but undoubtedly possible. (One edge would be easy because it is along the x-axis.)

restart

with(plots); with(plottools)

prob of m and k successes in N trials with prob p,q per trial

trinomial := proc (m, k, N, p, q) options operator, arrow; factorial(N)*p^m*q^k*(1-p-q)^(N-m-k)/(factorial(m)*factorial(k)*factorial(N-m-k)) end proc

proc (m, k, N, p, q) options operator, arrow; factorial(N)*p^m*q^k*(1-p-q)^(N-m-k)/(factorial(m)*factorial(k)*factorial(N-m-k)) end proc

Ternary coords for m = N at (N,0), k=N at (1/2,sqrt(3)/2)*N, third one =N at the origin. - see  https://en.wikipedia.org/wiki/Ternary_plot

tcoords := proc (m, k) options operator, arrow; (2*m+k)/2., sqrt(3)*k/2. end proc

proc (m, k) options operator, arrow; (2*m+k)/2., sqrt(3)*k/2. end proc

triangle := proc (N) options operator, arrow; polygon([[tcoords(0, N), 0], [tcoords(N, 0), 0], [0, 0, 0]], color = yellow) end proc

pins := proc (N, p, q) options operator, arrow; [seq(seq(line([tcoords(m, k), 0], [tcoords(m, k), trinomial(m, k, N, p, q)], color = red, thickness = 3), m = 0 .. N-k), k = 0 .. N)] end proc

Warning, (in pins) `k` is implicitly declared local

Warning, (in pins) `m` is implicitly declared local

lbls := proc (N) options operator, arrow; textplot3d({[0, 0, 0, cat("var3 = ", N), align = [below, right]], [(1/2)*N, (1/2)*sqrt(3)*N, 0, cat("var2 = ", N), align = [below, right]], [N, 0, 0, cat("var1 = ", N), align = [below, left]]}) end proc

N := 20

20

display(pins(N, .2, .5), triangle(N), lbls(N), axes = normal, labels = ["", "", prob], axis[1, 2] = [color = white])

NULL

Download trinomial.mw

Here's how I would follow the Youtube simplifications in Maple land. I also show that term1 and term2 are the some parameters. Edit: I mistakenly thought that @Carl Love was using t for x, but I now see he was referring to a different post.
 

restart

We want to show that term 1 = term2. We follow the procedure here

term1 := (2*cos(2^n*x)+1)/(2*cos(x)+1)

(2*cos(2^n*x)+1)/(2*cos(x)+1)

term2 := product(2*cos(2^k*x)-1, k = 0 .. n-1)

product(2*cos(2^k*x)-1, k = 0 .. n-1)

Let's write out the factors in the product up to n=6

f := [op(eval(term2, n = 6))]

[2*cos(x)-1, 2*cos(2*x)-1, 2*cos(4*x)-1, 2*cos(8*x)-1, 2*cos(16*x)-1, 2*cos(32*x)-1]

We multiply the first factor by g (the first factor but with a sign change); just (a-b)*(a+b) = a^2-b^2

g := 2*cos(x)+1; f1 := expand(g*f[1])

2*cos(x)+1

4*cos(x)^2-1

We show it is equal to the second factor but with a sign change.

f1 := combine(f1)

1+2*cos(2*x)

So similarly multiplying by the second factor and then combining we get the third factor with a sign change

f1*f[2]; f2 := combine(%)

(1+2*cos(2*x))*(2*cos(2*x)-1)

2*cos(4*x)+1

So continuing up to the last term we will get

fn_minus_1 := 2*cos(2^n*x)+1

2*cos(2^n*x)+1

which we need to divide by the g we multiplied by originally to get the final result

result := fn_minus_1/g

(2*cos(2^n*x)+1)/(2*cos(x)+1)

Maple gives this form immediately if we evaluate term2 for a given n, but can't do the case of general n.

eval(term1, n = 6)

(2*cos(64*x)+1)/(2*cos(x)+1)

Check for some specific params

params := {n = 1, x = (1/6)*Pi}

{n = 1, x = (1/6)*Pi}

rationalize(eval(term1, params)); eval(term2, params)

3^(1/2)-1

3^(1/2)-1

NULL

Download MapleLand.mw

OK, I now understand you want the plot of u(y) to look like the one you just posted, with a slope of zero at y=0 and going down to u(1) = 1.
Now psi = D(u), but in the odeplot that you labelled u(x,y) you asked to plot D(psi), which is not u. u is instead the integral of psi.

pltN[j] := odeplot(
        dsol[j], [y, D(psi)(y)], y = 0..1,
        labels = [typeset(y), typeset(u(x, y))], labelfont = ["TIMES", "italic", 18],

So just add another ODE to your system

diff(u(y),y) = psi(y), u(1)=1 

and plot u(y).

Edit: Of course you will have to change your boundary conditions that assumed D(psi) = u instead of psi = D(u)

Yes, that seems like a problem. Using LinearSolve is faster. (If you can use a float calculation NullSpace is much faster.)

restart; M := LinearAlgebra:-RandomMatrix(200, 500, generator = -50 .. 50)

st := time(); linalg[kernel](M); nops(%); time()-st

32.109

st := time(); X := LinearAlgebra:-LinearSolve(M, `<,>`(`$`(0, 200))); time()-st

24.984

``

nops(indets(X))

300

NULL

NULL

NULL

Download NullSpace_vs_kernel.mw

Edit: If you want it in the same format as NullSpace:

NullSpace.mw

What do you want to plot? You could use plot3d to plot abs(FF) vs x and P1 or Re(FF) vs x and P1 or Im(FF) vs x and P1 or argument(FF) vs x and P1. Here I use complexplot 3d to plot abs(FF) coloured by the phase. The plot is not very informative, probably because your function has numbers with wildly different exponents, e.g. 10^(-2750) and 10^5 but these numbers have only about 5 digit accuracy. Probably you should go back to a version in symbolic form and then evaluate to numeric form with a high setting of Digits so you don't lose accuracy.

PLOT11.mw

Edit: If you really want to find where the function is zero with implicit plot, you should explore more limited regions - the plots above can be a guide for this. I suspect again the quality of your function is making it hard for the internal solver. 

In general this would just be a matter of adjusting the view option, giving both a horizontal and vertical range.

But I think here you want to see the intersection point between a line through the points [30000,0.2838019625*10^8], [60000,0.2122282820*10^8] and a line through the points [30000,0.2838047526*10^8], [60000,0.2126554523*10^8].

The y values are identical for the first few decimal places, so it will not be possible to visually dstinguish them.

In your first attempt you pass an expression to plot3d with the inverse fourier transform not evaluated, so at each point in the plot there will be an attempt to find it, which will fail. The solution is just to evaluate it before the plot, which can be done by assuming t>0. For inverse transforms you almost always need assumptions. (I added with(inttrans).)

ft1.mw

I don't quite understand point 1. Perhaps you want to be able to choose other values of S__1 and S__2? I would set

params := {S__1 = -1, S__2 = -1};

then use eval for the things you want to use these parameters in, e.g.,

eval(P, params);

This way, you only need to change the params line to run the worksheet with your new parameters.

For point 2, I'm not sure if you are asking how to find the conserved quantity, which I've shown how to do in the attached. If you are asking about your empty plot, that is because you need wider contour intervals, which I adjusted. You may want to try without the "contours =" option to see which contours Maple comes up with first, then if necessary make a finer adjustment by adding the contours option.

L1C.mw

It would be helpful if you used the same notation as the paper. Here's my attempt, but I don't get the result in the paper.

restart

with(PDEtools)

undeclare(prime, quiet)

declare(y(delta), quiet)

_local(gamma)

gamma

E9 := 2*zeta[2]*y(delta)^3+(-gamma[1]^2*zeta[2]-2*gamma[1]*gamma[2]*zeta[1]+zeta[3])*y(delta)+zeta[1]^2*zeta[2]*(diff(diff(y(delta), delta), delta))

2*zeta[2]*y(delta)^3+(-gamma[1]^2*zeta[2]-2*gamma[1]*gamma[2]*zeta[1]+zeta[3])*y(delta)+zeta[1]^2*zeta[2]*(diff(diff(y(delta), delta), delta))

E8 := (4*gamma[1]*zeta[2]/zeta[1]+2*gamma[2])*y(delta)^3+(-gamma[1]^2*gamma[2]+gamma[3])*y(delta)+(2*gamma[1]*zeta[1]*zeta[2]+gamma[2]*zeta[1]^2)*(diff(diff(y(delta), delta), delta))

(4*gamma[1]*zeta[2]/zeta[1]+2*gamma[2])*y(delta)^3+(-gamma[1]^2*gamma[2]+gamma[3])*y(delta)+(2*gamma[1]*zeta[1]*zeta[2]+gamma[2]*zeta[1]^2)*(diff(diff(y(delta), delta), delta))

tr := {diff(y(delta), `$`(delta, 2)) = X, y(delta) = Y}

{diff(diff(y(delta), delta), delta) = X, y(delta) = Y}

E8a := eval(E8, tr); E9a := eval(E9, tr)

(4*gamma[1]*zeta[2]/zeta[1]+2*gamma[2])*Y^3+(-gamma[1]^2*gamma[2]+gamma[3])*Y+(2*gamma[1]*zeta[1]*zeta[2]+gamma[2]*zeta[1]^2)*X

2*zeta[2]*Y^3+(-gamma[1]^2*zeta[2]-2*gamma[1]*gamma[2]*zeta[1]+zeta[3])*Y+zeta[1]^2*zeta[2]*X

E := eliminate({E8a, E9a}, X)

[{X = -Y*(2*Y^2*zeta[2]-gamma[1]^2*zeta[2]-2*gamma[1]*gamma[2]*zeta[1]+zeta[3])/(zeta[1]^2*zeta[2])}, {Y*(2*gamma[1]^3*zeta[2]^2+4*gamma[1]^2*gamma[2]*zeta[1]*zeta[2]+2*gamma[1]*gamma[2]^2*zeta[1]^2-2*gamma[1]*zeta[2]*zeta[3]-gamma[2]*zeta[1]*zeta[3]+gamma[3]*zeta[1]*zeta[2])}]

But Y <> 0 so

zeta[3] = solve(E[2][]/Y, zeta[3])

zeta[3] = (2*gamma[1]^3*zeta[2]^2+4*gamma[1]^2*gamma[2]*zeta[1]*zeta[2]+2*gamma[1]*gamma[2]^2*zeta[1]^2+gamma[3]*zeta[1]*zeta[2])/(2*gamma[1]*zeta[2]+gamma[2]*zeta[1])

Or, giving the same result,

solve({E8a, E9a}, {X, zeta[3]})

{X = -Y*(4*Y^2*gamma[1]*zeta[2]+2*Y^2*gamma[2]*zeta[1]-gamma[1]^2*gamma[2]*zeta[1]+gamma[3]*zeta[1])/(zeta[1]^2*(2*gamma[1]*zeta[2]+gamma[2]*zeta[1])), zeta[3] = (2*gamma[1]^3*zeta[2]^2+4*gamma[1]^2*gamma[2]*zeta[1]*zeta[2]+2*gamma[1]*gamma[2]^2*zeta[1]^2+gamma[3]*zeta[1]*zeta[2])/(2*gamma[1]*zeta[2]+gamma[2]*zeta[1])}

NULL

Download C1.mw

In this case it's probably easier to do the phase portrait with DEplot.

restart;

with(plots):

f := (u,v) -> -u+u^3;
g := (u,v) -> -2*v;

proc (u, v) options operator, arrow; -u+u^3 end proc

proc (u, v) options operator, arrow; -2*v end proc

The differential equations

de1 := diff(u(t),t) = f(u(t),v(t));
de2 := diff(v(t),t) = g(u(t),v(t));
 

diff(u(t), t) = -u(t)+u(t)^3

diff(v(t), t) = -2*v(t)

Phase portrait

DEtools:-DEplot({de1,de2},[u(t),v(t)],t=0..5, u=-2..2, v=-2..2);

Conserved quantities - Maple gives two functions, but they involve t separately from u(t) and v(t)
(In the earlier example one was an integral, which wasn't of interest.)

PDEtools:-ConservedCurrents({de1, de2}, [u(t), v(t)]);
cons:=op(rhs(op(%)));

[_J[t](t, u(t), v(t)) = f__1((-u(t)^2+1)*exp(-2*t)/u(t)^2, v(t)*exp(2*t))]

(-u(t)^2+1)*exp(-2*t)/u(t)^2, v(t)*exp(2*t)

The des are uncoupled and can easily be solved. For both solutions, v(t) decays exponentially and so nothing about v(t) alone is converved. In the second conserved quantity, Maple is telling us that if we multiplied v(t) by exp(2*t), the result would be independent of time, and we see that is would just be the constant v(0). Similarly for u(t).

desols:=dsolve({de1,de2,u(0)=u__0,v(0)=v__0},{u(t),v(t)});

{u(t) = 1/(-exp(2*t)*(u__0^2-1)/u__0^2+1)^(1/2), v(t) = v__0*exp(-2*t)}, {u(t) = -1/(-exp(2*t)*(u__0^2-1)/u__0^2+1)^(1/2), v(t) = v__0*exp(-2*t)}

But we are interested in functions of u(t) and v(t) alone that are conserved. If we look at the two conserved quantities, we see that the exp time dependence cancels out in their product:

con1:=simplify(cons[1]*cons[2]);

v(t)*(-u(t)^2+1)/u(t)^2

So if we put the first or second solution of the des in this expression we can check that it is indeed constant

simplify(eval(con1, desols[1]));
simplify(eval(con1, desols[2]));

v__0*(-u__0^2+1)/u__0^2

v__0*(-u__0^2+1)/u__0^2

Get rid of the t's so we can do the contour plot

con2 := eval(con1, {u(t)=u, v(t) = v});

v*(-u^2+1)/u^2

contourplot(con2, u=-2..2, v=-2..2, contours = [seq(-10..10,0.5)], colorbar = false);

NULL

Download conserved.mw

dchange from the PDEtools package can transform from one set of variables to another.

restart;

with(plots): with(PDEtools):

NULL

de1 := diff(x(t),t) = -y(t)+a*x(t)*(x(t)^2+y(t)^2);
de2 := diff(y(t),t) = x(t)+a*y(t)*(x(t)^2+y(t)^2);

diff(x(t), t) = -y(t)+a*x(t)*(x(t)^2+y(t)^2)

diff(y(t), t) = x(t)+a*y(t)*(x(t)^2+y(t)^2)

tr := {x(t) = r(t)*cos(theta(t)), y(t) = r(t)*sin(theta(t))};

{x(t) = r(t)*cos(theta(t)), y(t) = r(t)*sin(theta(t))}

newdes := dchange(tr, {de1,de2}, [r(t), theta(t)]);

{(diff(r(t), t))*cos(theta(t))-r(t)*(diff(theta(t), t))*sin(theta(t)) = -r(t)*sin(theta(t))+a*r(t)*cos(theta(t))*(r(t)^2*cos(theta(t))^2+r(t)^2*sin(theta(t))^2), (diff(r(t), t))*sin(theta(t))+r(t)*(diff(theta(t), t))*cos(theta(t)) = r(t)*cos(theta(t))+a*r(t)*sin(theta(t))*(r(t)^2*cos(theta(t))^2+r(t)^2*sin(theta(t))^2)}

solve(newdes,{diff(r(t),t), diff(theta(t),t)});
simplify(%, {sin(theta(t))^2 + cos(theta(t))^2=1});

{diff(r(t), t) = r(t)^3*a*(sin(theta(t))^2+cos(theta(t))^2), diff(theta(t), t) = 1}

{diff(r(t), t) = r(t)^3*a, diff(theta(t), t) = 1}

NULL

Download polar.mw

I agree it is a bug; it is still present in Maple 2025.1.

Edit: SCR submitted.

A simple workaround is just to substitute a non-indexed name before carrying out the command, and then reverse the substitution on the result; if you have more complicated requirements, I am sure there can be a more sophisticated solution.

restart;

with(PolynomialTools):

p:=(z[1]^2+3*z[1]+1)*f;

(z[1]^2+3*z[1]+1)*f

PolynomialToPDE([p], [z[1]], [f]);

[diff(diff(f(z[1]), z[1]), z[1])+4*f(z[1])]

p2 := eval(p, z[1] = z);
de2 := PolynomialToPDE([p2], [z], [f]);
de := eval(de2, z = z[1]);
 

(z^2+3*z+1)*f

[diff(diff(f(z), z), z)+3*(diff(f(z), z))+f(z)]

[diff(diff(f(z[1]), z[1]), z[1])+3*(diff(f(z[1]), z[1]))+f(z[1])]

 

NULL

Download PolynomialToPDE.mw

When you wrote puzzle(n) := etc you defined a procedure - notice that the output says puzzle := n -> etc. But you did not invoke that procedure - if you do you immediately get the answer. (I would define  puzzle by puzzle := n -> etc, since then the defining syntax is distinct from the calling syntax, but modern Maple allows both.)

restart

"puzzle(n):=(&sum;)(cos(k*Pi/(2*n+1))^())^(4);"

proc (n) options operator, arrow, function_assign; sum(cos(k*Pi/(2*n+1))^4, k = 1 .. n) end proc

puzzle(n)

(3/8)*n-(1/2)*cos((n+1)*Pi/(2*n+1))^4+((1/4)*cot(Pi/(2*n+1))+(1/8)*sec(Pi/(2*n+1))*csc(Pi/(2*n+1)))*cos((n+1)*Pi/(2*n+1))*sin((n+1)*Pi/(2*n+1))+((1/2)*cot(Pi/(2*n+1))-(1/4)*sec(Pi/(2*n+1))*csc(Pi/(2*n+1)))*cos((n+1)*Pi/(2*n+1))^3*sin((n+1)*Pi/(2*n+1))+(1/2)*cos(Pi/(2*n+1))^4-((1/4)*cot(Pi/(2*n+1))+(1/8)*sec(Pi/(2*n+1))*csc(Pi/(2*n+1)))*cos(Pi/(2*n+1))*sin(Pi/(2*n+1))-((1/2)*cot(Pi/(2*n+1))-(1/4)*sec(Pi/(2*n+1))*csc(Pi/(2*n+1)))*cos(Pi/(2*n+1))^3*sin(Pi/(2*n+1))

puzzle(2)

cos((1/5)*Pi)^4+cos((2/5)*Pi)^4

NULL

Download test.mw

@kencom1 In Maple 2025, your worksheet does not run as you show it; the notation f[0](1) when setting the initial conditions leads to an error. You are having problems using the same name for procedures and indexed names, which become tables and overwrite the procedures when you assign to them. The solution is to have different names, e.g. f for your procedure and F for your initial conditions and odes;

restart;

f := x -> x^2;

proc (x) options operator, arrow; x^2 end proc

type(f, procedure);

true

Unexpected

f[1](0) = 2;

0 = 2

Assigning to an indexed f  makes a table called f

f[1]:=3;

3

type(f, table);
type(f, procedure);

true

false

Now your procedure f no longer exists

f(x);

f(x)

NULL

Download notes.mw

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