@vv Thanks, I just assumed that was the default order. As you undoubtedly know, the exact order isn't important, just that it is consistent for the left and right eigenvectors.

[Fixed and comment added.]

@Christopher2222 Yes, you are right. Sorry for the confusion. The toggle input/output display works for plots and also if the output is on a different line because it was generated with "enter". But for equations where the input and output are on the same line, it seems you have to use view show/hidecontents and check/uncheck input or output. 

@Christopher2222 I only have 2024 right now but I think it is the same. After "plot(x^2)" then "ctrl=" you should have the plot on the same line, to the right of the command. Put the cursor in the plot command so that you are in the same document block (selecting the plot itself should also work). Then Edit->Document blocks->toggle input/output display makes the command disappear.

@mmcdara local gamma; evalf(gamma); returns gamma and not the number in Maple 2024, i.e., the bug is fixed.

In general, use simplify(A-B) to ckeck if A and B are the same. But your worksheet has several problems. You cannot use N and N[1] as separate variables; once you assign to N[1], N no longer has the value you set previously. You can avoid this issue by using different names, which can be done by making subscripts by, for example, N__1. The same problem occurs with x and x[01], y and y[01] and maybe other places.

Do not use Digits:=1, which means to calculate with one digit accuracy. If you just want to display 1 digit, use Tools->Options->precision->round screen display.

You should not compare two indefinite integrals, since the arbritary constant could be different in the two cases.

@MaPal93 I didn't go through it all - just two places you asked questions about.

derivatives_final.mw

Two second derivatives change sign. For each, I find the zero with solve() but it returns two zeros instead of just one as I expected. Why?

I don't know - perhaps it is a solution on another branch of Lambda. Use fsolve if you just want a numerical answer.

I find minima and maxima with [DirectSearch]GlobalOptima, which is most of the time accurate but quite slow. Any alternative that is more efficient? I write "most of the time" because occasionally it gives me some very random results, completely off from the peaks and troughs I can eyeball from the plots.

GlobalOptima is overkill here - you only want a local maximum - Search is fast here, but needs an initialpoint. But here I would just use fsolve on the derivative with a range or initial guess.

My axes are always gamma*sigma_d^2 vs gamma*sigma_d*sigma_v. Are there any other param combinations that I can use? Would it make sense to scale the plots together with the actual axes, e.g., dividing both by gamma*sigma_d so that my new y-axis and x-axis are, respectively, sigma_d and sigma_v? How to do so? By doing so perhaps the plots are easier to interpret...

Yes, I wondered about that. Some people would non-dimensionalize the problem right at the beginning and then everything would be in terms of diminsionless quantities and the reader would be left to convert back to the real quantities. Back at the beginning you had f1 =Lambda*(sigma__v/sigma__d), so let's write a = sigma__v/sigma__d. and you also have Gamma = gamma*sigma__v*sigma__d, let's write b = gamma*sigma__v*sigma__d. So all quantities you write have to just have combinations of a and b. Here your axes are b/a and b. So for a universal plot with  axes sigma__d and sigma__v, those quantities would have to be expressible in terms of a and b, which I don't think is possible (you can't get rid of gamma unless you just take combinations of a alone, which can't be separated into sigma__v and sigma__d). Another way to think about this is to suppose you had such a plot, and you chose a point on the line for certain values of sigma__d and sigma__v. Does that result apply for all gamma? So the value of gamma seems essential to the problem, which may make sense to you in terms on the original problem.

@sand15 I usually use it in combinatorial applications where being zero for b>0 makes mathematical sense. Wikipedia gives a "definition" that it is the coefficient in the expansion (1 + x)^alpha = sum(binomial(alpha,k)*x^k,k=0..infinity), for any complex alpha and |x|<1. Then  (1+x)^2 = 1*1 + 2*x + 1*x^2 + 0*x^3 + ...

DLMF defines it for complex alpha (=z) by 1.2.6, so binomial(2,3) = 2(2-1)..(2-3+1)/3! = 2(2-1)(2-2)/3! = 0

You are right about the conversions, but I think it is usual when converting between special functions that some cases might be exceptions.

@MaPal93 Looking nicer? You know what they say about beauty.

derivatives_2.mw

1. A plot is worth a thousand lines of code. It shows that you are correct and is is wrong (a bug, though my expectations for is on complicated expressions are low).

2. Your 2nd derivatives look correct. Unless the functions are pathological, the symmetry holds.

3. You forgot with(plots); display and textplot are in this package.

@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.

derivatives.mw

@MaPal93 

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

working_example2.mw

• 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:

something.mw

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:

@MaPal93 

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.

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.