# Mathematica: Optimizing NIntegrate with complex integrand

Tags:
1. Dec 4, 2011

### KarolisK

Hi, first of all I am new to Mathematica.(Mathematica v8.0.0)
The Problem:
I have been having some issues with one particular integral lately:
I need to perform InverseFourier over a region of certain frequency and set the result as a function of time for later use:
$$A(t)=\frac{1}{2\pi}\int_{-3ω}^{3ω} S(ω,f_1(ω),f_2(ω,z),z) \cdot e^{i \omega t} dω$$
$$\text{where } S(ω,f_1(ω),f_2(ω))\text{ - is complex function proportional to } \omega,\omega^2, e^{\omega}.\text{ z - real constant}$$

I managed to get the solution by writing a line:
Code (Text):
A[t_?ArrayQ] := (1/(2*\[Pi]))*NIntegrate[S[w, f1[w], f2[-w, 3540], 3540]*Exp[I*w*t ], {w, -3*ws, 3*ws},
Method -> "GlobalAdaptive", WorkingPrecision -> 4];
As seen in the code input t will be an array of real numbers. Later on since A(t) is complex I square it and plug in t as an array(z=A[tarray];), however it takes some time to calculate the result for an array of 62 values. The problem is I am quite convinced this is not the best programming form of performing the integral calculation. I think it reintegrates the integral with every new t value. So is there anyone who can help me with this problem? I would be very grateful for your help!
As a comparison Mathcad calculates the result in ~1s

2. Dec 4, 2011

### KarolisK

Ok, so its probably because the precision is too big. Now im confused how to make the calculations with a precision to 0,001 parts of a number. So in case of 0.0019, mathematica would return 0.002?

3. Dec 4, 2011

### KarolisK

Ok, It seems I need to use AccuracyGoal for this, however it does not work for me, or am I doing something wrong? I would imagine AccuracyGoal would set the result accuracy to 4 digits, no?

File size:
67.7 KB
Views:
81
4. Dec 5, 2011

### KarolisK

Ok, since noone is replying and it is kinda akward to write one more reply to myself, I am not sure if I have stated the problem clearly. Therefore, I have made a minimum working example ready for compilation. So, as mentioned before I am stuck because the evaluation of NIntegrate takes too long (comments in the MWE). I need to reduce the precision and accuracy of the NIntegrate or rewrite it in some way? Hopefully, I have made it clear this time :) . Thank you in advance.

MWE:
Code (Text):

ClearAll["Global`*"];

(*---Constants---*)
tau = 20/Sqrt[2*Log[2.]]; gam = 0; Subscript[A, 0] = 1;
(*---Constants---*)

S0[w_] := (Sqrt[\[Pi]]*tau*Subscript[A, 0]/Sqrt[1 + I*gam])*
Exp[-1*(w*tau)^2/(4*(1 + I*gam))];
(*-------------------------------------------------*)

(*---Constants---*)
Subscript[\[Nu], 13] = 0.081; Subscript[\[Nu], 23] = 0.093;
gk1 = 0.066; gk2 = -0.915;
\[CapitalDelta]k = 1.257*10^(-3);
(*---Constants---*)

D1[w_] :=
Subscript[\[Nu], 13]*w + (1/2)*gk1*w^2 + 0*\[CapitalDelta]k/2;
D2[w_] :=
Subscript[\[Nu], 23]*w + (1/2)*gk2*w^2 - 0*\[CapitalDelta]k/2;
(*-------------------------------------------------*)
\[CapitalDelta]1[w_, sug1_] := D1[w] - I*sug1;(************)
\[CapitalDelta]2[w_, sug2_] := D2[w] - I*sug2;(************)
(*-------------------------------------------------*)
(*---Constants---*)
Subscript[\[Alpha], 10] =
0*0.01823; Subscript[\[Beta], 1] = -0.03*0; Subscript[\[Epsilon], 1] \
= 0.630*0; Subscript[\[Epsilon], 12] = -8*0; Subscript[\[Epsilon], \
13] = 22*0;
(*---Constants---*)
\[Alpha]1[w_] :=
Subscript[\[Alpha], 10] + Subscript[\[Beta], 1]*w +
Subscript[\[Epsilon], 1]*
w^2 + (Subscript[\[Epsilon], 12]*w^2 +
Subscript[\[Epsilon], 13]*w^4);
appsug[w_] := 0.6 + (1/(3 + Exp[-82*(w - 0.02)]));
\[Alpha]2[w_, Lopa_] := -2*Log[appsug[w]]/Lopa;
(*-------------------------------------------------*)
Ga[w_, sug1_, sug2_] :=
1/2*(\[CapitalDelta]1[w, sug1] +
Conjugate[\[CapitalDelta]2[-w, sug2]]);
(*-------------------------------------------------*)
(*---Constants---*)
Subscript[\[CapitalGamma], 0] = Sqrt[14.5]/3120;
(*---Constants---*)
\[CapitalGamma]a[w_, sug1_, sug2_] :=
Sqrt[Subscript[\[CapitalGamma], 0]^2 - Ga[w, sug1, sug1]^2];
\[CapitalGamma]1[w_, sug1_,
sug2_] := -I/
2*(\[CapitalDelta]1[w, sug1] -
Conjugate[\[CapitalDelta]2[-w, sug2]]) + \[CapitalGamma]a[w,
sug1, sug2]; \[CapitalGamma]2[w_, sug1_,
sug2_] := -I/
2*(\[CapitalDelta]1[w, sug1] -
Conjugate[\[CapitalDelta]2[-w, sug2]]) - \[CapitalGamma]a[w,
sug1, sug2];
(*-------------------------------------------------*)
S1a[w_, sug1_, sug2_,
Lopa_] := (S0[w]/2)*(1 -
I*(Ga[w, sug1, sug2]/\[CapitalGamma]a[w, sug1, sug2]))*
Exp[\[CapitalGamma]1[w, sug1, sug2]*Lopa] + (S0[w]/2)*(1 +
I*(Ga[w, sug1, sug2]/\[CapitalGamma]a[w, sug1, sug2]))*
Exp[\[CapitalGamma]2[w, sug1, sug2]*Lopa];
(*-------------------------------------------------*)
P1[w_, Lopa_] :=
Abs[S1a[w, \[Alpha]1[w], \[Alpha]2[-w, Lopa], Lopa]]^2;
(*-------------------------------------------------*)
ws = 2/tau; imax = 164; \[CapitalDelta]w = (2/imax)*ws; wxi =
Array[-ws + \[CapitalDelta]w*(# + 1) &, {imax - 1}];
tmax = 60; tbs = 30*20; ddt = (2/tmax)*tbs; txi =
Array[-tbs + ddt*(1 + #) &, {tmax - 1}];
(*-------------------------------------------------*)
(*-------------------------------------------------*)

(*---Constants---*)
Lopa1 = 3540;
(*---Constants---*)
Aac[t_?ArrayQ] := (1/(2*\[Pi]))*
NIntegrate[
S1a[w, \[Alpha]1[w], \[Alpha]2[-w, Lopa1], Lopa1]*
Exp[I*w*t], {w, -3*ws, 3*ws}];
envlp = Abs[
Aac[txi]]^2;(*----Takes too much time to evaluate, such a
precision of values is not needed-----*)
(*-------------------------------------------------*)
data = Transpose[{txi, envlp}];
ListLinePlot[data, PlotRange -> {{0, 600}, {0, 120}}, Frame -> True]

I Have also attached the notebook version if needed.

File size:
26.8 KB
Views:
65