@dharr thank you for the useful note.
Given the Jacobian, is this correct? compact_derivatives_putting_all_together.mw

@dharr thanks for clarifying. Yes Veil[] works good here before differentiation: compact_derivatives_Veil.mw

EDIT: there must be some smarter usage of alias(), but the script above is a start...

@dharr "Following RootOf is in the derivatives, so alias it.": where is it? isn't it the case that it's in the derivatives of X__2 only because it's in the derivatives of X__1 but which we haven't specified yet? Why would we care about any RootOf() in the derivatives of X__2 before even specifying X__1? Isn't X__1 the only "source" of any RootOf()?

I was hoping to get something nicer looking with

***with(LargeExpressions):
compact := collect(dX2dy1, [X__1(y__1,y__2),Diff(X__1(y__1,y__2),y__1)], [Veil[A],Veil[B]] );

but no output...

***this command is maybe wrong: I don't know how and I never used collect(,veil[]) with two wrappers. I was thinking of A[i] as the coefficients on X__1(y__1,y__2)-terms and B[i] as the coefficients on Diff(X__1(y__1,y__2),y__1)-terms. (and similarly for dX2dy2.)

@janhardo mmm, not too clear to me. Would you mind explain? Do you mean I need to insert diff(Lambda_1,Gamma_1) and diff(Lambda_1,Gamma_2) somewhere in diff(Lambda_2,Gamma_1) and diff(Lambda_2,Gamma_2) for the latter diff() to "capture"?

it would be useful if you could work with my code. Thanks!

@dharr understood, thanks.

@dharr with yours I encounter Error, Got internal error in Typesetting:-Parse : "non-algebraic exponent: rho = 0"

this is what works for me instead:

InertForm:-Display(`%>`(%limit(X[i]^`ρ=0`, alpha = infinity), 0), inert = false)

@dharr thank you for your quick reply.

@dharr how are the sample points in m:-SamplePoints chosen? In particular, why the sample points for regions 2 and 3 are chosen on the y=0 axis? (Why not choose any other point on the respective region?)

@mmcdara thank you for the great insights. I know more about it now and I started to understand why and when one would need a GE. I'll do some research myself but yes of course if you happen to have/come across some useful references feel free to share here...perhaps would be of interest to other readers as well. Thanks again!

@mmcdara I am back on this. Thank you for taking a look so far. I have been thinking more about this thread.

To summarise, fsolve() fragility depends on:

1. Starting points: if set wrong, these cause fsolve() to fail at the onset
2. Search range: if set wrong, these cause fsolve() to fail sometime later when the lower bound of either of the 6 variables is hit (recall that I am only interested in positive solutions)
3. Step: if set wrong, this causes fsolve() to fail sometime later (see point 2.) or/and generate a discontinuity which is most likely an artifact (but not always, see @acer comment below about the residual check precision)

Having said that, while I lowered my expectations about finding acceptable solutions over a wider range of INDEX (I am kind of giving up...), I am still puzzled when I can't even find a solution at the very first iteration: 

• another_example_A.mw
• another_example_B.mw

If you happen to have some time to look at these 2 it would be great. I am puzzled because the equations are literally the same as before (and the fixed parameters are set to the same values) except that I increase linearly sigma_v1 and rho, respectively (instead of sigma_delta). And yet, I can't find a way to get fsolve() "started".

Thanks. 

@mmcdara thank you so much!! This procedure is super useful!! 

Regarding the Kriging, thank you also for the additional work. It's extremely precise indeed, but the interpolation function (last output) is extremely complicated. Then what do you mean by "the interpolator depends on a single parameter named psi"?

Overall, from my current understanding, it seems that Kriging would be an overkill for curves which are essentially linear except for very small values of x, which is exactly the type of contourplot I obtain if I slice my surfaces 2 and 3 at z=0. So perhaps is not worth it for these two cases. Howeber, you show its outstanding accuracy for the Lennard-Jones potential-like curves, which is great, but the interpolation function looks much more convoluted than the simple LJ form...so again: what do you mean by "the interpolator depends on a single parameter named psi"?
 

@mmcdara if I run exactly the same script that you posted I get totally different accuracy and optimal parameters: 

Download display_formula_mmcdara_different_functional_form.mw

Why?

@mmcdara thanks.
 

1. Interpolation:-Kriging. I swear I did have a look at the help page and subpages but I didn't find anything useful. I would need a command like Spline/LeastSquares, that is, taking up the pts points as a list of lists and the independent variable (Gamma or Gamma_1 depending on my surface) and generating the functional form as some relatively simple model with optimized parameter values. Can you please help me? (I didn't find something like this). I am really interested in learning about Kriging, but if that's not feasible after all, what CurveFitting alternatives more sophisticated than LeastSquares would you suggest? Note that it doesn't have to be extremely accurate but at least I would like to capture the nonlinearity for very small values on the X-axis...feel free to suggest an "ad hoc" functional form like you did with the LJ potential...perhaps that would the simplest if you already have in mind a form which is mostly a line except at very small positive values on the X-axis...
2. grid=[N, N]. I confirm I do generate more points with a finer grid: finer_grid.mw, but other than a better looking plot (with visibly more interpolation points) the coefficients of the linear form from LeastSquare are mostly unaltered. Perhaps a finer grid would come more useful if I opt for a more flexible form (see point 1. above).
3. "Verify". Ok, I understand now.
4. From your other comment (in the "textplot question" of mine about displaying expressions as typed):

   model := a[1] + a[2]*((a[3]/Gamma)^a[4] - (a[5]/Gamma)^a[6]);

   (1)

   >	K := unapply(model, Gamma); J := add((map2(op, 2, pts) -~ K~(map2(op, 1, pts)))^~2): opt := evalf[6]( Optimization:-NLPSolve(J, assume=nonnegative) );

   Warning, limiting number of major iterations has been reached

   
   Is opt[1] the error/accuracy of the procedure? How did you achieve a better fitting (thus more optimal parameter values) compared to the original opt[1]=0.00164888? (given that the model and everything else is exactly the same as before except the naming of the 6 parameters)
   
   Again, thanks a lot for your time!

@Carl Love are you okay with the new "version for you"? 

@mmcdara it seems like all I had to do was to remove the last element in the list of lists: pts := select(L->L<>RGB,pts).

I tried to apply your script to the second surface that I mentioned in the question plus an additional one. This third surface is kind of specular (wrt y=x) to the second one...anyway both Spline and LeastSquare look good for these latter two surfaces.

Questions:

1. If I am not satisfied with the linear fit of LeastSquare (especially at the lower ends) but I still want to have a closed form expression for the approximation of my curves (like the L-J potential for the first curve), which functional form would you recommend for the second and third surface? (or which function in CurveFitting would you use?)
2. How to control the numer of interpolation points pts? (I would like to have more points for more granularity)
3. Why the estimations in the verification blocks have the wrong signs with respect to their true values?
    

Thanks a lot for your kind help: slices_of_three_surfaces.mw