2
$\begingroup$

There is probably a neat approach to solve this problem...but can't get to it at the moment. How do you maximize an integral with respect to a single parameter? My code below produces error messages. Please note I have amended the Wf[a_,b_] to a simpler form to make it more explicit.

 Wf[a_, b1_] := 1/2 Erfc[(-a + 0.1)/(Sqrt[2] b1)]


 Wfunc[al_, r_, s1_, s2_, m_] := 
 16*(Wf[s1*Sqrt[4 + (2*al*r*m)^2], s2]/(1 + (al*r*m)^2)^(3/2) - 
     Wf[s1*Sqrt[4 + (2*al*r*(1 - m))^2], 
       s2]/(1 + (al*r*(1 - m))^2)^(3/2))^2 

ListPlot[Table[{rm, 
   NIntegrate[
    y^2*Wfunc[y, rm, 0.015, 0.005, 0.71], {y, 0.01, 5}]}, {rm, 0.01, 
   2, 0.01}]]

 Maximize[
 NIntegrate[y^2*Wfunc[y, rm, 0.015, 0.005, 0.31], {y, 0.01, 5}], rm]

What I want is the value of rm that maximizes the Integral in the last line without error msgs.

enter image description here

(*  The integrand 16 y^2 (Erfc[141.421 (0.1 +Times[<<2>>])]/(2 (1+0.0961 Power[<<2>>] Power[<<2>>])^(3/2))-Erfc[141.421 (0.1 +Times[<<2>>])]/(2 (1+0.4761 Power[<<2>>] Power[<<2>>])^(3/2)))^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{0.01`,5}}. >> *)

(* {0.109851, {rm -> 1.20352}} *)
$\endgroup$
5
$\begingroup$

I would create another function for the integral as

int[rm_?NumericQ] := 
 NIntegrate[y^2*Wfunc[y, rm, 0.015, 0.005, 0.31], {y, 0.01, 5}]

This should prevent any non-numerical arguments getting into Nintegrate. Now we get the following data

Table[{rm, int[rm]}, {rm, 0.01, 2, 0.01}]

Also

In[10]:= NMaximize[int[rm], rm]

Out[10]= {0.109851, {rm -> 1.20352}}

Did you get a different rm value by using Maximize ?

$\endgroup$
  • $\begingroup$ Lotus.........I wonder why u missed out on the right rm value by a small amt? $\endgroup$ – thils Dec 1 '15 at 9:44
  • $\begingroup$ I am surprised myself. Especially since the Maximum value seems to be correct. $\endgroup$ – Lotus Dec 1 '15 at 9:46
  • $\begingroup$ Me too! I wonder what's happening here. $\endgroup$ – thils Dec 1 '15 at 9:47
  • $\begingroup$ It was a simple transcription error - you had .71 instead of .31, fixed it for you (hope you don't mind) $\endgroup$ – Jason B. Dec 1 '15 at 9:48
  • $\begingroup$ I am glad that this is resolved. $\endgroup$ – Lotus Dec 1 '15 at 9:51
3
$\begingroup$

You got the right answer, but it gave error messages because it tried to evaluate the integral - so take the integral out of the optimization step by building an interpolation function.

list = Table[{rm, 
   NIntegrate[
    y^2*Wfunc[y, rm, 0.015, 0.005, 0.31], {y, 0.01, 5}]}, {rm, 0.01, 
   10, 0.01}];
func = Interpolation[list];
Plot[func[x], {x, .1, 9}]

enter image description here

All of these will give the right answer,

Maximize[{func[x], 0 < x < 2}, x]
NMaximize[{func[x], 0 < x < 2}, x]
NMaximize[{func[x], 0.5 < x < 10}, x]
FindMaximum[func[x], x]
FindMaximum[func[x], {x, 1}]
(* {0.109851, {x -> 1.20353}} *)
(* {0.109851, {x -> 1.20353}} *)
(* {0.109851, {x -> 1.20353}} *)
(* {0.109851, {x -> 1.20353}} *)
(* {0.109851, {x -> 1.20353}} *)

But you have to visually inspect the plot to make sure you don't get the second, local maximum

NMaximize[{func[x], 0 < x < 10}, x]
FindMaximum[func[x], {x, 2.1}]
(* {0.104022, {x -> 2.60866}} *)
(* {0.104022, {x -> 2.60866}} *)
$\endgroup$
  • $\begingroup$ +1...let me wait & see if someone comes up with an alternative that does not involve interpolation... $\endgroup$ – thils Dec 1 '15 at 9:45
  • 1
    $\begingroup$ Yeah, Lotus's method is more correct I reckon, though not any faster than the interpolating method - sometimes I go for the fast and dirty method :-) $\endgroup$ – Jason B. Dec 1 '15 at 9:51

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.