Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica: Optimizing NIntegrate with complex integrand

Tags:
  1. Dec 4, 2011 #1
    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:
    [tex]
    A(t)=\frac{1}{2\pi}\int_{-3ω}^{3ω} S(ω,f_1(ω),f_2(ω,z),z) \cdot e^{i \omega t} dω
    [/tex]
    [tex]
    \text{where } S(ω,f_1(ω),f_2(ω))\text{ - is complex function proportional to } \omega,\omega^2, e^{\omega}.\text{ z - real constant}
    [/tex]

    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. jcsd
  3. Dec 4, 2011 #2
    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?
     
  4. Dec 4, 2011 #3
    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?
     

    Attached Files:

  5. Dec 5, 2011 #4
    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.
     

    Attached Files:

    • MWE.nb
      MWE.nb
      File size:
      26.8 KB
      Views:
      68
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Mathematica: Optimizing NIntegrate with complex integrand
Loading...