Mathematica: Optimizing NIntegrate with complex integrand

Click For Summary

Discussion Overview

The discussion revolves around optimizing the use of NIntegrate in Mathematica for a complex integral involving an Inverse Fourier transform. Participants explore issues related to computation time and precision when evaluating the integral over an array of time values.

Discussion Character

  • Technical explanation
  • Mathematical reasoning
  • Debate/contested

Main Points Raised

  • One participant describes their problem with a specific integral and expresses concern that NIntegrate may be recalculating the integral for each time value, leading to inefficiency.
  • Another participant suggests that the precision settings in Mathematica might be contributing to the long computation times and questions how to adjust the precision to a specific threshold.
  • A different participant proposes using the AccuracyGoal option to manage the precision of the results but expresses confusion about its implementation and effectiveness.
  • One participant shares a minimum working example (MWE) to illustrate their issue, emphasizing the need to reduce the precision and accuracy of NIntegrate to improve performance.

Areas of Agreement / Disagreement

Participants express various opinions on how to optimize the integral calculation, but there is no consensus on the best approach or solution to the problem. The discussion remains unresolved regarding the optimal settings for precision and accuracy in NIntegrate.

Contextual Notes

Participants note limitations related to the precision and accuracy of the integral calculations, but specific assumptions or dependencies on definitions are not fully explored. The discussion highlights the need for further clarification on the use of Mathematica's options for numerical integration.

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: 659
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 ·
Replies
1
Views
4K
  • · Replies 5 ·
Replies
5
Views
4K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 4 ·
Replies
4
Views
9K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 3 ·
Replies
3
Views
2K
Replies
3
Views
2K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 1 ·
Replies
1
Views
2K