0
$\begingroup$

I have to find the values of the parametersω and λ that minimize the function Int[ω, λ] defined as the sum of five quite complicated integrals of a trial function f[r, ω, λ], namely

a = -74.04252664070837;
b = 208.01432471151327;
d = -65.08706834153939;
A = 1.56692098226;
X0 = 3.9468;
X1 = 2.4323;
X2 = 3.6348;
f[r_, ω_, λ_] = r/Sqrt[r^2 + 2] (1 + 4 E^(-λ r) Sin[ω r])
Df[r_, ω_, λ_] = D[f[r, ω, λ], r];
j[x_, r_] = 2 A^2 r x;
g0[x_, r_] =  
  a + b/2 + 3 d/4 + A^2 (b + d) (x^2 + r^2) + d A^4 ((x^2 + r^2)^2 + 4 x^2 r^2);
g1[x_, r_] = j[x, r] (b + 2 d) + 4 d A^4 r x (r^2 + x^2);

Int1[ω_, λ_] = NIntegrate[(Df[r, ω, λ]^2 + f[r, ω, λ]^2/r^2) r, {r, 0, 100}];
Int2[ω_, λ_] = 
  Pi^1.5/A 
    NIntegrate[
      ((1 - f[r, ω, λ]^2) (1 - f[x, ω, λ]^2) r x E^(-A^2 (r^2 + x^2)) 
        (g0[x, r] BesselI[0, j[x, r]] - g1[x, r] BesselI[1, j[x, r]])), 
      {x, 0, Infinity}, {r, 0, Infinity}, 
      Method -> "LevinRule"] 
Int3[ω_, λ_] = 
  -1/2 X0 
    NIntegrate[(1 - f[r, ω, λ]^2)^2 r, {r, 0, Infinity}, Method -> "LevinRule"];
Int4[ω_, λ_] = 
  -1/3 X1 
     NIntegrate[((1 - f[r, ω, λ]^2)^2 (2 + f[r, ω, λ]^2)) r, {r, 0, Infinity}, 
       Method -> "LevinRule"];
Int5[ω_, λ_] = 
  1/4 X2 
    NIntegrate[
      (1 -f[r, ω, λ]^2)^2 (f[r, ω, λ]^4 + 2 f[r, ω, λ]^2 + 3) r, {r, 0, Infinity}, 
      Method -> "LevinRule"] 
Int[ω_, λ_] = Int1[ω, λ] + Int2[ω, λ] + Int3[ω, λ] + Int4[ω, λ] + Int5[ω, λ];

I can't use any minimization algorithm since Mathematica doesn't seem to be able to minimize $Int[\omega,\lambda]$. At least I would like to plot or evaluate (also a table) $Int[\omega,\lambda]$ in a reasonable region for the parameters $0<\omega<2$ and $0<\lambda<2$ but Mathematica returns the values of very slowly $Int[\omega,\lambda]$ and for certain values in this region crashes and my pc with it.

Also, the output shows many errors. Since I'm quite new to Mathematica, I actually don't know how to solve the problem nor the best method of integration for Int2 and Int5 which seem to be the most difficult integrals to evaluate. I used also "LevinRule" because the output of Mathematica seemed to suggest that option.

$\endgroup$
  • 2
    $\begingroup$ In your posted code you have not defined ω and λ before calling NIntegrate. $\endgroup$ – Anton Antonov Jul 8 '18 at 1:18
  • $\begingroup$ @AntonAntonov I don't understand what to you mean: do I have to give a value for \[Omega] and \[Lambda] before calling NIntegrate ? Also in you redefined code are missing expressions for Int2, Int3, Int4,Int5 so you evaluated only Int1 while my main problem are Int2 and Int5 $\endgroup$ – edinorog Jul 8 '18 at 9:26
2
$\begingroup$

As always, you should carefully verify that any result is correct.

On a fairly old slow desktop machine with only 12gb with a fresh clean evaluation of this:

a = -74.04252664070837;
b = 208.01432471151327;
d = -65.08706834153939;
A = 1.56692098226;
X0 = 3.9468;
X1 = 2.4323;
X2 = 3.6348;
Clear[r, ω, λ, f, Df, j, g0, g1];
f[r_, ω_, λ_] = r/Sqrt[r^2 + 2] (1 + 4 E^(-λ r) Sin[ω r]);
Df[r_, ω_, λ_] = D[f[r, ω, λ], r];
j[x_, r_] = 2 A^2 r x;
g0[x_, r_] = a + b/2 + 3 d/4 + A^2 (b + d) (x^2 + r^2) + d A^4 ((x^2 + r^2)^2 + 4 x^2 r^2);
g1[x_, r_] = j[x, r] (b + 2 d) + 4 d A^4 r x (r^2 + x^2);
Clear[Int1, Int2, Int3, Int4, Int5];
exp1 = Simplify[(Df[r, ω, λ]^2 + f[r, ω, λ]^2/r^2) r];
exp2 = Simplify[((1 - f[r, ω, λ]^2) (1 -  f[x, ω, λ]^2)* r x E^(-A^2 (r^2 + x^2))*
  (g0[x, r] BesselI[0, j[x, r]] - g1[x, r]*BesselI[1, j[x, r]]))];
exp3 = Simplify[(1 - f[r, ω, λ]^2)^2 r];
exp4 = Simplify[((1 - f[r, ω, λ]^2)^2 (2 + f[r, ω, λ]^2)) r];
exp5 = Simplify[(1 - f[r, ω, λ]^2)^2 (f[r, ω, λ]^4 + 2 f[r, ω, λ]^2 + 3) r];
bound=10;
Int1[ω_, λ_] := NIntegrate[exp1, {r, 0, bound}, Method -> "MonteCarlo"];
Int2[ω_, λ_] := Pi^1.5/A*NIntegrate[exp2, {x,0,bound}, {r,0,bound}, Method->"MonteCarlo"];
Int3[ω_, λ_] := -1/2 X0*NIntegrate[exp3, {r, 0, bound}, Method -> "MonteCarlo"];
Int4[ω_, λ_] := -1/3 X1*NIntegrate[exp4, {r, 0, bound}, Method -> "MonteCarlo"];
Int5[ω_, λ_] := 1/4 X2*NIntegrate[exp5, {r, 0, bound}, Method -> "MonteCarlo"];
Int[ω_, λ_] := Int1[ω, λ] + Int2[ω, λ] + Int3[ω, λ] + Int4[ω, λ] + Int5[ω, λ];

Look for a minimum near ω==0.5, λ==1.2

SortBy[Flatten[Table[{ω, λ, Int[ω, λ]}, {ω, 1/4, 3/4, 1/16}, {λ, 1, 3/2, 1/16}], 1], Last]
ListPlot3D[%]

enter image description here

and there does not appear to be any clear minimum near your estimate.

Can you do the work to determine what upper bound is sufficient for your integrands to be very close to zero beyond that?

$\endgroup$
  • $\begingroup$ > Can you reproduce exactly what I have done here without making any changes? Yes I can. But if I try to evaluate the other integrals as you did, my computer really slows down and mathematica crashes, now (using your code) without any output. My laptop is not that bad, so I don't really understand what the problem is. You managed to evaluate the sum for more than few points? $\endgroup$ – edinorog Jul 8 '18 at 10:32
  • $\begingroup$ Also using my previous code I was able to find the value of the sum for a few points, but now with your code (adding the definition of parameters and function g0, g1 etc) , can't get nothing at all. $\endgroup$ – edinorog Jul 8 '18 at 10:40
  • $\begingroup$ it took quite a long time but yes I did. the problem now is that I really dont want ` ω=0 ` as a minimum because it would imply that my trial function doesn't oscillate at all which is a thing that I don't expect from theory. Do you think that choosing a different bound for the integral (maybe 10 is too small, specially for Int1) results could be better? I say that I would expect a minimum at about ` ω=0.5, λ=1.2`..Thank you for the efforts. $\endgroup$ – edinorog Jul 8 '18 at 21:21
  • $\begingroup$ Also Int1 diverges as bound$\to\infty$ so the bound is necessary while the other integrals do not need a bound (maybe only to reduce computation time) or we can choose a larger bound.. $\endgroup$ – edinorog Jul 8 '18 at 21:29
  • 1
    $\begingroup$ @Ale Lenoci Edit doubled the speed $\endgroup$ – Bill Jul 9 '18 at 16:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.