1
$\begingroup$

hope you're having a nice day!

Context: I'm trying to constrain parameters from a model of Orbital Resonance equilibrium. In short, I know a parameter "$\Delta_{eq}$" which is derived from observations. I know a set of parameters {$m_1$,$m_2$,$M_*$,$q$,$f_1$,$f_{2p}$} and I know (not precisely yet) physical boundaries for variables $\alpha$, $hr_1$, $hr_2$, $\Sigma_1$, $\Sigma_2$, given $\Delta_{eq}$. Explicitly

$$\Delta_{eq}=\sqrt{2} \sqrt{\frac{m_1 m_2 \left(\alpha m_2 (q+1)+m_1 q\right) \left(\text{hr}_1^4 q^2 \Sigma _2 f_{2 p}^2+\alpha ^{3/2} f_1^2 \text{hr}_2^4 (q+1)^2 \Sigma _1\right)}{\text{hr}_1^2 \text{hr}_2^2 M_*^2 q^2 (q+1)^2 \left(\text{hr}_1^2 m_2 \Sigma _2-\sqrt{\alpha } \text{hr}_2^2 m_1 \Sigma _1\right)}}$$

My try: Evaluating known values (observables) I can explore the parameters I'm interested in, first plotting the real part of

$$0.000492702 \sqrt{\frac{(0.000704487 \alpha +0.000933589) \left(5.6644 \alpha ^{3/2} \text{hr}_2^4 \Sigma _1+2.84934 \text{hr}_1^4 \Sigma _2\right)}{\text{hr}_1^2 \text{hr}_2^2 \left(0.000352244 \text{hr}_1^2 \Sigma _2-0.000933589 \sqrt{\alpha } \text{hr}_2^2 \Sigma _1\right)}}$$

and Manipulating

Manipulate[Plot[Re[0.0004927023775077934* Sqrt[((0.0009335890199999999 + 0.00070448742*\[Alpha])* (5.6644*\[Alpha]^(3/2)*Subscript[hr, 2]^4*Subscript[\[CapitalSigma], 1] + 2.849344*Subscript[hr, 1]^4*Subscript[\[CapitalSigma], 2]))/ (Subscript[hr, 1]^2*Subscript[hr, 2]^2*(-0.0009335890199999999*Sqrt[\[Alpha]]* Subscript[hr, 2]^2*Subscript[\[CapitalSigma], 1] + 0.00035224371*Subscript[hr, 1]^2* Subscript[\[CapitalSigma], 2]))]], {\[Alpha], 0, 0.5}, PlotRange -> All], {Subscript[hr, 1], 0.01, 0.05}, {Subscript[hr, 2], 0.01, 0.05}, {Subscript[\[CapitalSigma], 1], 10^(-6), 10^(-3)}, {Subscript[\[CapitalSigma], 2], 10^(-6), 10^(-3)}]

What, at the time of this post, I've not been able to study is the constrained space of parameters that satisfy $\Delta_{eq} \approx 0.035$, of my particular interest. That is to say, solving:

sol1 = Reduce[ Re[0.0004927023775077934*Sqrt[((0.00070448742*\[Alpha] + 0.0009335890199999999)* (5.6644*\[Alpha]^(3/2)*Subscript[hr, 2]^4*Subscript[\[CapitalSigma], 1] + 2.849344*Subscript[hr, 1]^4*Subscript[\[CapitalSigma], 2]))/(Subscript[hr, 1]^2* Subscript[hr, 2]^2*(0.00035224371*Subscript[hr, 1]^2* Subscript[\[CapitalSigma], 2] - 0.0009335890199999999*Sqrt[\[Alpha]]*Subscript[hr, 2]^2* Subscript[\[CapitalSigma], 1]))]] == 0.035 && 1/10^6 < Subscript[\[CapitalSigma], 2] < Subscript[\[CapitalSigma], 1] < 1/10^3 && 0 < \[Alpha] < 1 && 0.009 < Subscript[hr, 1] < 0.1 && 0.009 < Subscript[hr, 2] < 0.1, {\[Alpha], Subscript[\[CapitalSigma], 1], Subscript[\[CapitalSigma], 2], Subscript[hr, 1], Subscript[hr, 2]}, Reals]

Question 1: What can I do to efficently solve the problem? I've also tried with NSolve and constraining only $\Sigma$s.

Question 2: solving only for $\Sigma$s I've found that the solution is $\{\}$ empty for guessed $hr$s so, in the case that there's no solution for $\Delta_{eq}=0.035$, how can I find the minimum $\Delta_{eq}$? Also, is there any way to solve for $\Delta_{eq} \pm \delta\Delta$ given $\delta\{m_1,...,f_{2p}\}$ uncertainties?

Disclaimer: Tried to copy the code as Raw Input following this question but couldn't do better. I apologize if the code is hard to read.

And by the way, I appreciate any suggestion/alternative regarding on the problem

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ Do you have multiple observations? If you have many you might determine the variables hr1, ... by optimization (e.g., =NonlinearModelFit) $\endgroup$ Commented Oct 30, 2021 at 13:55
  • $\begingroup$ @JackLaVigne This is not the case but it would be useful to consider it when possible. Thanks! $\endgroup$ Commented Oct 30, 2021 at 14:05

2 Answers 2

Reset to default
3
$\begingroup$

I am going to change the name of some of your variables, removing the subscripts. You can find many questions and answers that show that using variables with subscripts can create problems, see this answer.

Here is a function using your variables where I have, for example, substituted Subscript[hr, 1] with hr1, ...

fun[α_, hr1_, hr2_, Σ1_, Σ2_] := 0.0004927023775077934* Sqrt[((0.0009335890199999999 + 0.00070448742*α)*(5.6644*α^(3/2)* hr2^4*Σ1 + 2.849344*hr1^4*Σ2))/(hr1^2* hr2^2*(-0.0009335890199999999*Sqrt[α]* hr2^2*Σ1 + 0.00035224371*hr1^2*Σ2))]

Here is the Traditional form

TraditionalForm[fun[α, hr1, hr2, Σ1, Σ2]]

Mathematica graphics

Now to get variables in the ball park try this approach. Create a table iterating over the variables with the ranges from your Manipulate.

data = Table[ {α, hr1, hr2, Σ1, Σ2, fun[α, hr1, hr2, Σ1, Σ2]}, {α, 0, 0.5, 0.1}, {hr1, 0.01, 0.05, 0.01}, {hr2, 0.01, 0.05, 0.01}, {Σ1, 10^-6, 10^-3, (10^-3 - 10^-6)/5}, {Σ2, 10^-6, 10^-3, (10^-3 - 10^-6)/5} ];

data will have the dimensions

Dimensions@data (* {6, 5, 5, 6, 6, 6} *)

Flatten it

data = Flatten[data, 4];

This is a list of 5400 sub lists.

Select only those whose 6th element is Real

dataReal = Select[data, Im[#[[6]]] == 0 &];

This reduces it to 2929 sub lists.

Find the sub list whose absolute value is closest to 0.035.

minDataReal = Min[Map[(# - 0.035)^2 &, dataReal[[All, 6]]]] (* 5.30542*10^-9 *) Position[Map[(# - 0.035)^2 &, dataReal[[All, 6]]], minDataReal] (* {{1236}} *)

This occurs at position 1236 in dataReal.

dataReal[[1236]] (* {0.1, 0.04, 0.05, 251/1250000, 4001/5000000, 0.0349272} *)

Now plot it.

Module[ { hr1 = 0.04, hr2 = 0.05, Σ1 = 251/1250000, Σ2 = 4001/5000000 }, Show[ Plot[fun[α, hr1, hr2, Σ1, Σ2], {α, 0, 0.5}, PlotRange -> All, PlotStyle -> Black ], Graphics[{ Red, PointSize -> 0.025, Point[{0.1, 0.0349272}] }] ] ]

Mathematica graphics

Notice that hr2 hits an upper limit while the rest of the variables are within your bounds. This suggests that you might want to raise the upper limit of hr2.

You can fine tune the estimate by adjusting min, max and step of the variables.

If one sets α one is able to successfully run FindInstance.

FindInstance[ { fun[0.1, hr1, hr2, Σ1, Σ2] == 0.035, 0.009 < hr1 < 0.1, 0.009 < hr2 < 0.1, 10^-6 < Σ1 < 10^-3, 10^-6 < Σ2 < 10^-3 }, {hr1, hr2, Σ1, Σ2}, Reals, 5 ] (* {{hr1 -> 0.0754329, hr2 -> 0.0478999, Σ1 -> 0.000213849, Σ2 -> 0.000247237}, {hr1 -> 0.0845794, hr2 -> 0.0433808, Σ1 -> 0.0000638252, Σ2 -> 0.0000988564}, {hr1 -> 0.0845794, hr2 -> 0.0465941, Σ1 -> 0.000702261, Σ2 -> 0.000702261}, {hr1 -> 0.0775397, hr2 -> 0.0673684, Σ1 -> 0.000731934, Σ2 -> 0.000731934}, {hr1 -> 0.0911567, hr2 -> 0.0991313, Σ1 -> 0.000616756, Σ2 -> 0.000741676}} *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks a lot! @JackLaVigne Could you please give me a hint on why Reduce & NSolve didn't solve it? Also, do you think I could find a set of sublists with abs. values between Delta+- err easily? I'm not familiarized with Select, Map & Position. $\endgroup$ Commented Oct 30, 2021 at 22:59
  • 1
    $\begingroup$ You can get values that are within a particular error range via: Select[dataReal, (#[[6]] - 0.035)^2 < 10^-6 &]. Here I am using 10e-06 as the upper error value. You are supplying one equation with 6 unknowns and 4 two sided constraints. NSolve needs a determined system and Reduce runs without stopping on my system. It appears from the tabular data that you have the Sigma1 and Sigma2 constraints wired backwards. $\endgroup$ Commented Oct 31, 2021 at 0:09
  • 1
    $\begingroup$ FindInstance works if one sets \[Alpha]. You might want to try setting another variable while holding \[Alpha] constant. $\endgroup$ Commented Oct 31, 2021 at 0:50
1
$\begingroup$

Using Jack's improved function fun, you could also try NMinimize to solve your problem and use the "RandomSearch" method and varying the seed to produce a large list of solution candidates. From here you can then study the histograms of individual variables from the good solutions. This should help you gauge what the solution space looks like, though note there will be bias in the sampling which might favour certain parts of the solution space :

fun[α_, hr1_, hr2_, Σ1_, Σ2_] := 0.0004927023775077934* Sqrt[((0.0009335890199999999 + 0.00070448742*α)*(5.6644*α^(3/2)* hr2^4*Σ1 + 2.849344*hr1^4*Σ2))/(hr1^2* hr2^2*(-0.0009335890199999999*Sqrt[α]* hr2^2*Σ1 + 0.00035224371*hr1^2*Σ2))] target = 0.035; solutionCandidates = Table[ NMinimize[{Abs[(fun[α, hr1, hr2, Σ1, Σ2] - target)^2], 0 < α < 0.5, 0.01 < hr1 < 0.05, 0.01 < hr2 < 0.05, 10^-6 < Σ1 < 10^-3, 10^-6 < Σ2 < 10^-3}, {α, hr1, hr2, Σ1, Σ2}, Method -> {"RandomSearch", "RandomSeed" -> i, "SearchPoints" -> 10}] , {i, 30}]; (* select the best solutions with the lowest error *) solutions = Select[solutionCandidates, #[[1]] < 10^-10 &][[All, 2]]; (* study the histogram of a particular parameter *) Histogram[Σ2 /. solutions]
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.