@MaPal93 Another simple relationship, which may be obvious to you from the formulation of the problem is

where the X__i are the values at the maximum.

MatrixNew8.mw

@dharr Elsewhere (and at the bottom of the attached worksheet), the calibration case was shown not to have any positive solutions, thus providing a counterexample for the original conjecture. What can be shown more generally?

I worked up the rest of the problem into a matrix based form.

MatrixNew7.mw

Ultimately it turns into a problem about the signs of diagonal elements that doesn't have clear-cut answers. However, along the way (assuming the lambda's were positive) I derived a (to me) surprising simple result:

So if there is some reason to believe the betas are positive, then the lambdas can't be positive and we have arrived at a contradiction, disproving the conjecture that there exists a positive solution (for all cases). If the betas can have any sign then one has to work harder. Looking at the key expressions for the three lambdas:

The denominators are quadratic forms of positive semidefinite matrices and therefore are non-negative. So for positive lambdas we need positive numerators. For the calibration case with rho__v zero and sigma__v1 = sigma__v2, this immediately leads to beta__11 positive, beta__22 positive and beta__31+beta__32 positive. Perhaps this can be worked out more fully into some contradiction between the signs of betas and lambdas, or some cases allowing positive solutions.

For real solutions other than positive ones, the matrix approach doesn't have much to say, since positive lambdas were key to making Q positive definite. For that reason, it is not clear to me how the original conjecture could be modified.

Edit: thnk more about the cases where the denominators are zero.

@delvin I don't know much about these functions, but I presume you need JacobiSN(zeta[1]*phi,ell) etc where ell is the modulus for your situation.

@MaPal93 The conditions on lambda__2 arise fairly naturally from just making sure that there are no zero denominators in the solutions, so I think you just need to first find the solutions and then just look at how the conditions relate to the solutions. Here's my analysis, which leads only to complex solutions  :(

MatrixNew3.mw

@MaPal93 I would suggest complete workup with your calibration system without putting in numerical values for what I called parameters (@mmcdara's "remain" variables). Since the restrictions on lambda__2 are not that onerous, just find when the quadratic equation in Edit: corrected: lambda__3 has real, positive solutions, remove any that don't obey the conditions on lambda__2, and then work out lambda__1 for the OK lambda__2's, and see if/when they are positive.

Maybe this case is simple enough that solve/parametric can do this for you.

Once you understand the calibration system, then you can try some less restrictive ones, where the conditions and components might be more complicated.

@mmcdara I followed only this thread4 (different from the other 4!), and worked extensively on a matrix formulation, but got no response so I lost interest. I know very little statistics, but think that sometimes a matrix formulation can be easier than working with polynomials. Like you, I somehow expect that there must be an easier way, and I am slowly grasping the interrelationships between the different quantities.

This time the matrix Q has a lot fewer lambda's in, and so there might be a better chance of some solutions that are tractable. So I was willing to work a bit more on it.

@MaPal93 Perhaps you can provide @mmcdara with the notes you sent me, or upload them here.

@MaPal93 In the earlier problem, (separate thread), I recall that PolynomialSystems sometimes produced three components, so you had to solve the third, then the second, then the first, whereas above there are fewer "restrictions" (badly chosen word) since you can choose lambda__2 and lambda__3 arbitrarily. Since you are still working on a simpler calibration system, perhaps there will be more components for other cases. Anyway, I'll shortly post a separate answer that solves the calibration system without the denominators problem (and has two components).

@acer Fantastic. Vote up.

@emersondiaz I assume you are talking about the dsolve(DE,...) solution. Your worksheet doesn't reproduce the error you mention (I instead get a message about a singulaity) so I don't know how you got that. So here are some comments:

1. The default method automatically adjust the stepsize, so you should not use the stepsize option.

2. If your equations occur with many diffferent timescales (or rscales here), then you may need a stiff solver - add stiff=true (I routinely use this).

3. Terms in the equations and the numbers you are interested in for r and R are many orders of magnitude away from 1 (10^6 and 10^45 appear with 10^(-31). For example (r^2 - 1.690735000*10^37*r + 0.81). So the equations need to be rescaled (non-dimensionalize would be typical, or use a different unit solution) to avoid these numbers.

Point 3 is the most important point.

@jalal 

Like this:

Download PEMDASTest.mw

@NIMA112 You want now to solve a PDE on two regions with matching conditions at the common boundary, Eqs 6.30-6.33. As far as I know, Maple's numerical PDE solver only handles single rectangular regions. A Tool like COMSOL is probably better suted to this. Or you could follow the procedure given to find the complicated series solution; I didn't try to understand that in detail.

@mmcdara Thanks. Of course I first tried the obvious elimination of all of v,a,b,R,T, which doesn't work. v turns into Z to make the cubic, so I tried that first and made progress. Then since A is a nondimensionalized version of a, and similarly for b I added them and the R and T disappeared because they "had to" as they weren't in the final result.

I use solve more than I use eliminate, but here I would have got a RootOf solving for Z: solve({eqA, eqB, eqP, eqZ}, {Z, a, b, v}).

@mmcdara You analysis appeared as I was working on the numerical solution. Vote up. 

@Carl Love Actually, I wasn't relying on that in the procedure. If it does remember (probably most times), then the loop executes faster, but otherwise it is just slower.

I added a note about it in the worksheet.

@ider I think @sursumCorda's point is that by using a truncation order of zero you only get the principal part and the O() term, which you can remove by convert/polynom. So the following should give the principal part:

convert(series(BesselK(4, x), x = 0, 0),polynom);

You suggested it sometimes stops earlier, which I didn't understand. Do you have an example of that?

(The documentation for laurent suggests some subtlety when there are an infinite number of negative powers, which I didn't follow.)