# Double NIntegrate

The functions I have are

i[d_,\[Tau]_,g_,w_]=-FractionBox["1", "8"] (-1+SqrtBox[RowBox[{"1", "-", SuperscriptBox["d", "2"]}]]) (2 SuperscriptBox["g", "2"]+SuperscriptBox["w", "2"]-SuperscriptBox["w", "2"] Cos[2g*\[Tau]])

a[\[Tau]_,g_,d_]=ArcCos[FractionBox["1", "4"] SqrtBox[RowBox[{"1", "-", "d"}]] (SqrtBox[RowBox[{"1", "-", "d"}]]+SqrtBox[RowBox[{"1", "+", "d"}]]+(SqrtBox[RowBox[{"1", "-", "d"}]]-SqrtBox[RowBox[{"1", "+", "d"}]]) Cos[g *\[Tau]])+FractionBox["1", "4"] SqrtBox[RowBox[{"1", "+", "d"}]] (SqrtBox[RowBox[{"1", "-", "d"}]]+SqrtBox[RowBox[{"1", "+", "d"}]]+(-SqrtBox[RowBox[{"1", "-", "d"}]]+SqrtBox[RowBox[{"1", "+", "d"}]]) Cos[g *\[Tau]])]

b[t_,w_,g_,d_]=Sqrt[-FractionBox["1", "8"] (-1+SqrtBox[RowBox[{"1", "-", SuperscriptBox["d", "2"]}]]) (2 SuperscriptBox["g", "2"]+SuperscriptBox["w", "2"]-SuperscriptBox["w", "2"] Cos[2g*t])]

c[\[Tau]_,w_,g_,d_]:=NIntegrate[b[t,w,g,d],{t,0,\[Tau]}]

e[\[Tau]_,g_,d_,w_]:=a[\[Tau],g,d]*\[Tau]/c[\[Tau],w,g,d]


Where all parameters except t and g are fixed (tau=4, d=0.5, w=5)

(sorry for the format but if you copy it to Mathematica you get the right functions)

Now as you can see my function c is a numerical Integral and my function e contains c as one of the building blocks. As vizaluisation here is the Code that plots my function e (just so you see that in theory a mean value should be existent):

\[Tau]=4
d=0.5
w=5
Plot[{e[\[Tau],g,d,w]},{g,0,5},AxesLabel->{Style["g",Bold,16]},PlotLabels->{Style["\[Tau](QSL)",Bold,13]}]


My Problem now is that I want to plot the Mean of function e over the interval and to do that I would need to calculate something along the Lines of

1/j*NIntegrate[e[\[Tau]_,g_,d_,w_],{g,0,j}]


and to plot this function

Plot[{1/j*NIntegrate[e[\[Tau]_,g_,d_,w_],{g,0,j}]},{j,0,5}]


which is the mean of e over the variable g in [0,5]. But if I enter a code above in Mathematica it gives me a ton of errors. I guess this is due to the fact that mathematica has no explicit expression for the function I want to integrate because it in itself contains a numerical integration.

A double Integral does not work because my first integration of c is the denominator of the function e and I dont want to Integrate over 1/c by just putting

1/j*NIntegrate[{a[\[Tau],g,d]*\[Tau]/(Sqrt[-FractionBox["1", "8"] (-1+SqrtBox[RowBox[{"1", "-", SuperscriptBox["d", "2"]}]]) (2 SuperscriptBox["g", "2"]+SuperscriptBox["w", "2"]-SuperscriptBox["w", "2"] Cos[2g*t])])},{{t,0,\[Tau]}},{g,0,j}]


some help yould be highly apreciated. Thank you in advance!

• Convert your code to InputForm prior to copy and paste into MSE. Jul 15 at 13:47
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jul 15 at 15:59
• You've got a lot of these  in your code. What symbol is this supposed to be? Jul 16 at 9:35

Clear["Global*"]

τ = 4;
d = 1/2; (* use exact values for known constants *)
w = 5;

i[d_, τ_, g_,
w_] = -(1/8) (-1 + Sqrt[1 - d^2]) (2 g^2 + w^2 - w^2 Cos[2 g*τ]) //
Simplify;

a[τ_, g_, d_] =
ArcCos[1/4 Sqrt[
1 -
d] (Sqrt[1 - d] + Sqrt[
1 + d] + (Sqrt[1 - d] - Sqrt[1 + d]) Cos[g*τ]) +
1/4 Sqrt[
1 +
d] (Sqrt[1 - d] + Sqrt[
1 + d] + (-Sqrt[1 - d] + Sqrt[1 + d]) Cos[g*τ])] // Simplify;

b[t_, w_, g_, d_] =
Sqrt[-(1/8) (-1 + Sqrt[1 - d^2]) (2 g^2 + w^2 - w^2 Cos[2 g*t])] //
Simplify;


Since c uses a numeric technique, its arguments should be restricted to being numeric. To control precision, use arbitrary precision rather than machine precision.

c[τ_?NumericQ, w_?NumericQ, g_?NumericQ, d_?NumericQ] :=
Module[{ww, gg, dd, ττ},
{ww, gg, dd, ττ} = SetPrecision[{w, g, d, τ}, 20];
NIntegrate[b[t, ww, gg, dd], {t, 0, ττ}, WorkingPrecision -> 15]]

e[τ_, g_, d_, w_] := a[τ, g, d]*τ/c[τ, w, g, d]

Plot[{e[τ, g, d, w]}, {g, 0, 5},
AxesLabel -> {Style["g", Bold, 16]},
PlotLabels -> {Style["τ(QSL)", Bold, 13]}]


mean[j_?Positive] := Module[
{jj = SetPrecision[j, 15]},
1/jj*NIntegrate[e[τ, g, d, w],
{g, 0, jj},
WorkingPrecision -> 12,


Calculation of mean is slow

(data = Table[{j, mean[j]}, {j, 0.05, 5, 0.075}];) // AbsoluteTiming

(* {240.857, Null} *)

ListLinePlot[data]


EDIT Following the recommendation by Akku14 to redefine c

Clear[c, e, mean]

c[τ_, w_, g_, d_] = Assuming[g >= 0,
Integrate[b[t, w, g, d], {t, 0, τ}]] // Normal

(* 1/4 (-1 + Sqrt[3]) EllipticE[4 g, -(25/g^2)] *)

e[τ_, g_, d_, w_] := a[τ, g, d]*τ/c[τ, w, g, d]

mean[j_?Positive] :=
Module[{jj = SetPrecision[j, 15]},
1/jj*NIntegrate[e[τ, g, d, w], {g, 0, jj}, WorkingPrecision -> 12]]

Plot[mean[j], {j, 1/100, 5}]


• Thanks a lot for your answer. I did not even know that there are such possibilities to modify a calculation. I will try to calculate it in my mathematica but I think my calculation time is limited to 1 min. Thank you! Jul 15 at 15:19
• In this case c  can be done with Integrate  analyticaly, which makes things easier and faster. Jul 15 at 19:35
• @Akku14 - Thanks. Edited to use Integrate` Jul 15 at 22:53