-3
$\begingroup$

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

enter image description here

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!

$\endgroup$
3
  • 4
    $\begingroup$ Convert your code to InputForm prior to copy and paste into MSE. $\endgroup$
    – Bob Hanlon
    Jul 15 at 13:47
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Jul 15 at 15:59
  • $\begingroup$ You've got a lot of these  in your code. What symbol is this supposed to be? $\endgroup$
    – flinty
    Jul 16 at 9:35
0
$\begingroup$
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]}]

enter image description here

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

Calculation of mean is slow

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

(* {240.857, Null} *)

ListLinePlot[data]

enter image description here

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

enter image description here

$\endgroup$
3
  • $\begingroup$ 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! $\endgroup$ Jul 15 at 15:19
  • $\begingroup$ In this case c can be done with Integrate analyticaly, which makes things easier and faster. $\endgroup$
    – Akku14
    Jul 15 at 19:35
  • $\begingroup$ @Akku14 - Thanks. Edited to use Integrate $\endgroup$
    – Bob Hanlon
    Jul 15 at 22:53

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.