Dr. David Harrington

6078 Reputation

21 Badges

19 years, 320 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 replies submitted by dharr

@dharr Here are the derivatives, cast in a form that easily shows the combination of parameters that can be used to scale to make them dimensionless.



Most looks good. See comments on the worksheet as well. (I redid some parts where you did it "by hand", just to check.)


  • However, I have sigma_d3/sigma_d on the y-axis instead of sigma_d/sigma_d3 as for the lambda1-lambda3 comparison you did. Does it make sense to have an inverted y-axis or should I try to have the same y-axis to facilitate interpretation?

It's up to you; easy to take the reciprocal if you want. Sometimes one way up seems more natural than the other, depending on the problem.

  • For the alpha-alpha_s comparison, I found that alpha > alpha_s regardless of the parameter values.

agreed (for the absolute values)

  • However, I am having a tough time doing the beta-alpha and beta-alpha_s comparisons. Can you help?

They have different "units" so I wouldn't normally do a comparison, but it is possible - just plot K vs Gamma.

I will take a look at the derivatives in a while. I don't have a strong opinion about moving to a new thread; I can respond here.

@mary120 OK, so now using the initial condition x(phimax/2) = 0 works for the original function:


Notice that it doesn't go far in the +x direction because your function is not well behaved near zero, and it doesn't go far in the -x direction because the values become complex. Both of these issues are to do with original issue, which is the accuracy of your function and its mix of very large and very small numbers. If you have a parameter close to 1e-14 and another near 4000, then you need to give all your constants to at least 17 digit accuracy, or better, redefine your equation in terms of different variables that don't lead to numbers many orders of magnitude away from 1.

For the two functions in your last post, your plot of F(phi) does not show a region with a minimum or maximum, so fsolve cannot find phimax, where the derivative is zero. So find a region were the plot looks the way you expect, then proceed from there.

@dharr This version is more efficient to generate the whole sequence up to some value, and uses the q-series expansions directly.

congruentlist := proc(nmax::posint)
	local mgmax, ngmax, m2max, m4max, nn, g, theta2, theta4, gth2, gth4,
            an, bn, n, m, i, j, q, c, x;
	description "Find list of congruent numbers up to possible value nmax using Tunnell's theorem";
	# ensure enough terms in the q-series for nmax to be reliable
	ngmax := ceil(sqrt(nmax/8));
	mgmax := ceil((sqrt(nmax) + 1)/4);
	m2max := ceil(sqrt(nmax/2));
	m4max := ceil(sqrt(nmax/8));
	g := add(add((-1)^n*q^((4*m + 1)^2 + 8*n^2), m = -mgmax .. mgmax), n = -ngmax .. ngmax);
	theta2 := 1 + 2*add(q^(2*m^2), m = 1 .. m2max);
	theta4 := 1 + 2*add(q^(4*m^2), m = 1 .. m4max);
	gth2 := expand(g*theta2);
	an := [seq(coeff(gth2, q, i), i = 1 .. nmax, 2)];
	gth4 := expand(g*theta4);
	bn := [seq(coeff(gth4, q, i), i = 1 .. nmax/2)];
	j := 0;
	for n to nmax do
		# divide out largest square and use squarefree residue
		nn:= mul(map(x -> x[1]^(x[2] mod 2), ifactors(n)[2]));
		if nn::odd and (an[(nn+1)/2] = 0) or nn::even and (bn[nn/2] = 0) then
			j := j + 1;
			c[j] := n;
		end if;
	end do;
	convert(c, list)			
end proc:



The method given was tailored for that specific example. It gives a zone diagram for when f1 is greater/less/equal than f3 for the two key combined parameters of the problem, which are plotted on the two axes. At the heart of the problem is the non-dimensionalized function Lambda(Gamma), which is "dressed up" with the additional parameters of the problem. So the manipulations were to find the relevant key parameters. You say that the scaling was just to plot it, but that amounts to the same thing - finding a single plot that encapsulates the data so you don't need a series of plots for different parameters. The key is to solve for when f1 = f3 (still true for scaled versions), which gives a line separating the different zones - then one can test the various regions to see whether they are for f1>f3 or f1<f3. If I undestand what you want, its give a complete solution for what combinations of parameters lead to f1 </>/= f3.

The same philosophy applies to partial derivatives, because these are all dressed up versions of diff(Lambda(Gamma),Gamma) via the chain rule. For combinations of functions, you might have more key parameters and a more complicated situations, or you might have just two.

You ask how to fix complicated_comparison.mw, but really it is giving you what you asked for. The LambertW appears for equations involving the form x = exp(-a*x), which is a key part of your approximate version, with x = f3. There isn't a simpler analytical solution. I could have solved for f3 and made a plot of LambertW functions, but it is simpler to use the parametric plot. (My experience is that LambertW functions are tricky because of the different branches.) But the key point is that this is the line on some plot. Notice that in asking for solve to use the various inequalities you are asking for all cases of these parameters, and not products/ratios of them. If you first figure out the key parameter combinations and call them, say a and b, then this method might be more useful.

complicated_comparison2.mw didn't work because there are still parameters in g__1sc; a more careful analysis/manipulation of the parameters is required.

(I didn't look closely at partial_derivatives.mw)

@dharr Probably there is a more efficient way but this works.

Calculate the sides of  triangle having the congruent number n.
Procedure congruent is in the startup code.


  local a,b,c,i;
  if not congruent(n) then error "%1 is not a congruent number", n end if;
  for i do
  until c::rational;
end proc:


3, 4, 5


3/2, 20/3, 41/6


Download sides.mw

@JAMET congruent(3) returns false, which is correct; neither Wikipedia nor OEIS list 3 as congruent. For the sides calculation see below.

@jud I am using version 2023, which has it, but 2022 doesn't.

In response to your email request, here is my version.

Test if n is a congruent number by Tunnell's theorem.
(Assumes the Birch and Swinnerton-Dyer conjecture is true.)
Procedure congruent is in startup code.


We calculate the number of integer solutions An, Bn, Cn, Dn for the following equations.
For a squarefree integer:
If n is odd then it is congruent if 2*An = Bn
If n is even then it is congruent if 2*Cn = Dn
If n is not square-free, divide out largest square first.

Aeq := n = 2*x^2+y^2+32*z^2; Beq := n = 2*x^2+y^2+8*z^2; Ceq := n = 8*x^2+2*y^2+64*z^2; Deq := n = 8*x^2+2*y^2+16*z^2

n = 2*x^2+y^2+32*z^2

n = 2*x^2+y^2+8*z^2

n = 8*x^2+2*y^2+64*z^2

n = 8*x^2+2*y^2+16*z^2

From Wikipedia, the first congruent numbers are

L := [5, 6, 7, 13, 14, 15, 20, 21, 22, 23, 24, 28, 29, 30, 31, 34, 37, 38, 39, 41, 45, 46, 47, 52, 53, 54, 55, 56, 60, 61, 62, 63, 65, 69, 70, 71, 77, 78, 79, 80, 84, 85, 86, 87, 88, 92, 93, 94, 95, 96, 101, 102, 103, 109, 110, 111, 112, 116, 117, 118, 119, 120]

L2 := select(congruent, [`$`(1 .. L[-1])])

EqualEntries(L, L2)



Download congruent.mw

The code in the worksheet should work in older versions of Maple. In 2023 it can be done a little shorter as:

congruent := proc(n::posint)
	local nn, An, Bn, c1, c2, c3, c4, x, y, p, q, mx, my, mz;
	description "Determine if n is a congruent number using Tunnell's theorem";
	# divide out largest square and use squarefree residue
	nn := mul(map(x -> x[1]^(x[2] mod 2), ifactors(n)[2]));
	if nn::odd then
  		c1, c2, c3, c4 := 2, 1, 32, 8
		c1, c2, c3, c4 := 8, 2, 64, 16
	end if;
	An := 0;
	Bn := 0;
	for x from 0 while (p := nn - c1*x^2) >= 0 do
		for y from 0 while (q := p - c2*y^2) >= 0 do
			if issqr(q/c3) then
				# count all solutions e.g., for x nonzero there two possibilities of opposite sign
				mx := signum(x) + 1;
				my := signum(y) + 1;
				mz := signum(q) + 1;
				An := An + mx*my*mz;
			end if;
			if issqr(q/c4) then
				mx := signum(x) + 1;
				my := signum(y) + 1;
				mz := signum(q) + 1;
				Bn := Bn + mx*my*mz;
   			end if;
		end do
	end do:
	evalb(2*An = Bn);
end proc;



1- How can solve these two ODE(s) and  plot phi vs x in one phi-x coordinate?

Just solve them separately and generate two plots, in variables p1 and p2. Then display(p1,p2) will put them on the same plot.

2- How can I plot phi-x in a symmetric interval of x (for example, from x=-100..100).

The solver can go forwards and backwards from the initial point. But you didn't explain the significance of x=0 so it is not clear how the x axis is defined.

3- How can I export data of plot of these ODEs in the form of two distinct ascii table? (I couldn’t export data of the last figure in your last uploaded worksheet in the form of a ascii table)

I don't know why the export didn't work. If you upload a worksheet containing the error I can take a look. 

@mary120 if you assign the odeplot to a variable p, then M:=plottools:-getdata(p)[3]; gives a Matrix with the data. The matrix can then be exported with ExportMatrix to varrous  file types, e.g., target = csv will be comma separated. 

@mary120 To plot phi vs x just use


This goes through the origin because of the initial condition x(0)=0. But your plot doesn't go through the origin so I don't understand what you really want. Note the slope is wrong, so you probably need diff(x(phi),phi)=-integrand


Here is something, but it is unclear how the x axis is constructed - what does x=0 mean on your plot?


@Carl Love My original example was a graph with some randomly chosen edges. Later @Anthrazit wanted it to be a complete graph, and so I gave a version that does the complete graph for the same points, which does give the same path and distance as you find, see the worksheet TravelingSalesmanCompleteGraph.mw above.

I was also assuming that the start and end vertices were prescribed, and were to be given as the first and last in the list (since @Anthrazit mentioned something about nearest and outermost bolts), rather than just the shortest of all Hamiltonian paths. 

@Saha   You have gamma(Theta - Theta[a]) instead of gamma*(Theta - Theta[a])

@sand15 Yes! That syntax error can be avoided by diff(..., [eta$(m-2)]) , which returns the undifferentiated function. But probably not what the OP intended.

5 6 7 8 9 10 11 Last Page 7 of 62