Mathematica Mathematica: Optimizing NIntegrate with complex integrand

AI Thread Summary
The discussion revolves around optimizing the NIntegrate function in Mathematica for a complex integral that involves an Inverse Fourier transform. The user is experiencing slow computation times when evaluating the integral for an array of time values, suspecting that the integral is being recalculated for each time point. They seek advice on reducing the precision and accuracy of NIntegrate to speed up the calculations, as they compare the performance with Mathcad, which computes the result significantly faster. The user is unsure about the correct usage of the AccuracyGoal option to achieve the desired precision. Overall, the thread highlights the need for more efficient coding practices in Mathematica to handle complex integrals.
KarolisK
Messages
8
Reaction score
0
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:
<br /> A(t)=\frac{1}{2\pi}\int_{-3ω}^{3ω} S(ω,f_1(ω),f_2(ω,z),z) \cdot e^{i \omega t} dω <br />
<br /> \text{where } S(ω,f_1(ω),f_2(ω))\text{ - is complex function proportional to } \omega,\omega^2, e^{\omega}.\text{ z - real constant}<br />

I managed to get the solution by writing a line:
Code:
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
 
Physics news on Phys.org
Ok, so its probably because the precision is too big. Now I am 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?
 
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?
 

Attachments

  • Untitled.png
    Untitled.png
    58.1 KB · Views: 627
Ok, since no one 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:
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.
 

Attachments

Similar threads

Replies
1
Views
4K
Replies
5
Views
3K
Replies
2
Views
3K
Replies
4
Views
9K
Replies
3
Views
2K
Replies
2
Views
3K
Back
Top