## 135 Reputation

1 years, 318 days

## @dharr how are the sample points in...

@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?)

## thanks!! this is a great summary...

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

## 2 more examples...

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:

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.

## miscellaneous...

@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"?

## different optimal parameters?...

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

 > restart;
 > pts := [[1.020408163, .9618194785], [1.021860259, .9591836735], [1.047329717, .9183673469], [1.077147904, .8775510204], [1.112217401, .8367346939], [1.153643677, .7959183673], [1.202786007, .7551020408], [1.224489796, .7394037725], [1.266731737, .7142857143], [1.350499968, .6734693878], [1.428571429, .6426647396], [1.463093584, .6326530612], [1.632653061, .5927592700], [1.638566951, .5918367347], [1.836734694, .5653094517], [2.024654644, .5510204082], [2.040816327, .5499100928], [2.244897959, .5415404104], [2.448979592, .5375577894], [2.653061224, .5364038532], [2.857142857, .5371174931], [3.061224490, .5390909888], [3.265306122, .5419306191], [3.469387755, .5453753581], [3.673469388, .5492483720], [3.759428515, .5510204082], [3.877551020, .5532080969], [4.081632653, .5571718264], [4.285714286, .5612391621], [4.489795918, .5653739340], [4.693877551, .5695506292], [4.897959184, .5737510616], [5.102040816, .5779621428], [5.306122449, .5821743709], [5.510204082, .5863807930], [5.714285714, .5905762856], [5.775983717, .5918367347], [5.918367347, .5943431266], [6.122448980, .5978993957], [6.326530612, .6014214453], [6.530612245, .6049104007], [6.734693878, .6083674439], [6.938775510, .6117937612], [7.142857143, .6151905094], [7.346938776, .6185587956], [7.551020408, .6218996660], [7.755102041, .6252141001], [7.959183673, .6285030098], [8.163265306, .6317672398], [8.219364040, .6326530612], [8.367346939, .6345736540], [8.571428571, .6371910423], [8.775510204, .6397830180], [8.979591837, .6423509843], [9.183673469, .6448962158], [9.387755102, .6474198713], [9.591836735, .6499230055], [9.795918367, .6524065794], [10.00000000, .6548714697]]:
 > 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) );
 (2)
 > alias(STS = StringTools:-Substitute): string_model := convert(model, string); values       := evalf[3](opt[2])
 (3)
 > for n from 1 to 6 do   string_model := STS(string_model, convert(lhs(values[n]), string), convert(rhs(values[n]), string)) end do:
 > plots:-textplot([0, 0, cat("W__LJ = ", string_model)], axes=none,'font'=["helvetica","roman",25], size=[1000, 100])
 >

## "Replies to your replies" to points 1 to...

@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) );

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 &qu...

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

## wrong signs in 'Verify' blocks...

@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

## Am I correct?...

@mmcdara questions:

1. true_Gamma (=1.112) is the asymptote in the last plot?
2. the Lennard-Jones potential is the green curve? That is, the functional form K with optimal params a, b, c, d, f?
3. then, is my summary below correct?
• regardless of rho, for Gamma<true_Gamma my surface is always negative
• for Gamma>true_Gamma, my surface is negative if rho<green curve and positive if rho>green curve

Any misunderstanding?

Why don't you embed the L-J snippet in the worksheet and post your comment as an answer I can upvote/give best answer?

EDIT: encounter error when running it on my machine: Surface1_fsolve_workaround_mmcdara_MaPal.mw

## @acer I don't understand.Suppos...

@acer I don't understand.

Suppose this command:
plot3d(Z(X,Y),X=0.01..10,Y=0.01..10,'colorscheme'=["zgradient",["LightGray", "Gray", "Green"]], labels = [X, Y, Zeta(X,Y)], view=[default,default,0..10]);

What do you mean by

max( lowvalue, min( maxvalue,  expr ) )

around the expression?

## thanks...

@mmcdara it's the same numerical calibration exercise, but this time for sigma_delta2 (ncal3): discontinuity_persists_ncal3.mw

## Other reasons fsolve() breaks...

@mmcdara thank you for your ideas so far.

I don't want to spam this thread with another worksheet but I did find another instance of fsolve() breaking before the end of the for loop (besides the still unsolved and different issue with the discontinuity). Perhaps this is helpful to know for other contributors still on this (@Carl Love) or other readers.

Apparently hitting the lower bound of the chosen range for the solution is not the only reason fsolve() could "break". As such, dynamic adjustment of the lower bound doesn't always solve the issue. For this specific calibration I am referring to, the last valid iteration before the "break" has all its 6 lambdas well above their lower bounds.

"even if you find a good strategy for one calibration case there is no guarantee it still performs well for another one." Agree with you...as you argued in your previous comment perhaps the work required just to make this specific calibration work is not worth the effort.

Let's see if Carl Love has some other ideas in mind.

"Why do you believe positive solutions exist for INDEX > 0.8731? Do you feel they do or are you sure they really do?" I thought about it. I agree with you that it's not necessarily the case. Based on my understanding of the problem at hand, there should be positive and smooth solutions (for all 6 lambdas) for a range of sigma_* (* wildcard) values within (0,x] where x is unknown but surely <1. It seems like you found that such x is 0.8731 for this specific calibration. Of course x may differ for a different calibration.

P.S.: could you please check the update to my comment above? Why does the discontinuity persists for a broader plotting range?  Thanks.

## black-box procedures/operator-form inste...

@acer thanks for describing the inner workings of fsolve().

You write: "Moreover (and I know this sounds facile, but I'm being serious), it can help to pause and try to define the notion of a "root" -- in this float context. Once you have such, you can program for it."
I don't understand what you have in mind.

Also "Using black-box procedures can separate the working precision for the functional evaluations (Digits now set withing the scope of the procedures) from the accuracy requirement in the residual check." So how to use operator-form instead of expressions, for the equations? What exactly does it mean? Are you confident that it's the residual check at INDEX~0.554 the root cause of the singularity?

 1 2 3 4 5 6 7 Last Page 1 of 15
﻿