# Maximizing over an integral with a single parameter

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.

(*  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}} *)


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 ?

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

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}]


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}} *)
`
• +1...let me wait & see if someone comes up with an alternative that does not involve interpolation... – thils Dec 1 '15 at 9:45
• 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 :-) – Jason B. Dec 1 '15 at 9:51