@sursumCorda If you know the answer, you can cheat. Both convert(z*hypergeom([1/4, 1/2], [5/4], z**4), FormalPOwerSeries) and convert(EllipticF(z, I), FormalPowerSeries) give the same power series, verifying they are the same.

They are converted to apparently different expressions with integrals with convert(..., Int), so interconverting might be possible by rearranging, but it wasn't obvious to me.

@acer Nice solution! I'm a bit surprised the Iterator method is so much slower than @sursumCorda's more-or-less equivalent method with combinat routines.

Just to add to @mmcdara's comments. I especially agree with his last comment. Please clarify that you want to find each of  the 12 variables n1,n,2,n3,n4,4,a1,a2,a3,a4,s1,s2,s3,s4 in terms of the 4 parameters b,q,r,theta. So as noted you need 12 equations for your 12 variables. 

The solutions for b=0 have arbitrary values of n3,n4 and s4. Do you expect this to be true in the general case? A dimensional analysis might help you to reduce the number of variables or parameters.

The reason it works for b=0 is that it can be cast as a polynomial system for that case. That can be done also for the full system, which might make it feasible to solve. But it would be better to simplify first, and we definitely need number of eqns = number of variables. Are the ai, ni and si components of 4-vectors? If so, perhaps the original equation(s) in matrix vector form would be useful.

@dharr And here's how to do it with an event.

Download Event.mw

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