1
$\begingroup$

I am working on a Mathematica project where I want to plot the group delay. The problem is this code is running too slow. Is there any possibility of speeding up the code? Here is the original code:

Clear["Global`*"]
plotset = {FrameStyle -> Directive[Thickness[0.004]], 
      TicksStyle -> Directive[Black, 18]}; 

plotset2 = FrameTicksStyle -> {{Directive[Black, 18], Directive[FontOpacity -> 0, FontSize -> 0]}, {Directive[Black, 18], Directive[FontOpacity -> 0, FontSize -> 0]}};

G = 0.5 k1; \[Theta] = 0;    

\[CapitalPhi]m = 0;

 B = 5*10^-10;

\[Gamma] = 2*\[Pi]*28*10^9;

 N1 = Sqrt[3*10^16];

Em = Sqrt[5]/4*(\[Gamma])*N1*B;

\[Omega]c = 2*\[Pi]*7.86*10^9;

\[Omega]b = 2*\[Pi]*11.42*10^6;

\[CapitalDelta]a = \[Omega]b;

k1 = \[Pi]*3.35*10^6;

k2 = k1;

J = 0.5 k1;

kb = 300*\[Pi]; 

km1 = \[Pi]*1.12*10^6; 

km2 = \[Pi]*1.12*10^6;

\[CapitalDelta]m1 = \[Omega]b;

\[CapitalDelta]m2 = \[Omega]b; 

\[HBar] = 1.054*10^-34;

Subscript[g, 1] = 2 k1;  

Subscript[g, 2] = 2 k1;

wl = 2*\[Pi]*7.9*10^9;

sss = Sqrt[p1/(\[HBar]*wl)];

E1 = sss;

Ep  = 0.0001 E1;

\[Eta] = 0.5 ; 

gmb = 2*Pi;

h1 = -I*\[CapitalDelta]a + I*\[Delta] - k1;

SuperStar[h1] = I*\[CapitalDelta]a - I*\[Delta] - k1;

h2 = -I*\[CapitalDelta]a - I*\[Delta] - k1;

SuperStar[h2] = I*\[CapitalDelta]a + I*\[Delta] - k1;

h3 = -I*\[CapitalDelta]a + I*\[Delta] + k2;

SuperStar[h3] = I*\[CapitalDelta]a - I*\[Delta] + k2;

h4 = -I*\[CapitalDelta]a - I*\[Delta] + k2;

SuperStar[h4] = I*\[CapitalDelta]a + I*\[Delta] + k2;

h5 = -I*\[Omega]b + I*\[Delta] - kb;

SuperStar[h5] = I*\[Omega]b - I*\[Delta] - kb;

h6 = -I*\[Omega]b - I*\[Delta] - kb;

SuperStar[h6] = I*\[Omega]b + I*\[Delta] - kb;

h7 = -I*\[CapitalDelta]m1 + I*\[Delta] - km1;

SuperStar[h7] = I*\[CapitalDelta]m1 - I*\[Delta] - km1;

h8 = -I*\[CapitalDelta]m1 - I*\[Delta] - km1;

SuperStar[h8] = I*\[CapitalDelta]m1 + I*\[Delta] - km1;

h9 = -I*\[CapitalDelta]m2 + I*\[Delta] - km2;

SuperStar[h9] = I*\[CapitalDelta]m2 - I*\[Delta] - km2;

h10 = -I*\[CapitalDelta]m2 - I*\[Delta] - km2;

SuperStar[h10] = I*\[CapitalDelta]m2 + I*\[Delta] - km2;

\[Alpha] = ((k1 + I \[CapitalDelta]a)*((k2 - I \[CapitalDelta]a)) - 
    J^2)/(k2 - I \[CapitalDelta]a);

\[Beta] = (k1 - I \[CapitalDelta]a) + Subscript[g, 1]^2/(-I \[CapitalDelta]m1 + km1) + Subscript[g,2]^2/(-I \[CapitalDelta]m2 + km2) -  
   J^2/(k2 + I \[CapitalDelta]a);

\[Gamma]1 = ((\[Alpha]*\[Beta] - 4 G^2)/\[Beta]) + (Subscript[g,1]^2/(I \[CapitalDelta]m1 + km1)) + (Subscript[g, 2]^2/(I \[CapitalDelta]m2 + km2));
\[Psi]1 = (1/\[Gamma]1)*(E1 Sqrt[2*k1* \[Eta]] + ((2*G*Exp[I*\[Theta]])/\[Beta]*(E1 Sqrt[2*k1* \[Eta]])));

\[Psi]2 = (-I*Subscript[g, 2])/(km1 + I \[CapitalDelta]m1);

ms = \[Psi]1*\[Psi]2; 

F = ms*gmb;

\[Alpha]1 = h7*SuperStar[h6]* h5 + F^2*(SuperStar[h6] - h5) // Simplify;

\[Alpha]2 =    SuperStar[h2] + J^2/SuperStar[h4] + Subscript[g,1]^2/SuperStar[h10] // Simplify;

\[Alpha]3 = \[Alpha]2*SuperStar[h8] + Subscript[g, 2]^2 // Simplify;

\[Alpha]4 = 2*G*Exp[-I*\[Theta]]*I*Subscript[g, 2] + 
       I*Subscript[g, 2]*\[Alpha]2 // Simplify;

\[Alpha]5 = (-I*Subscript[g, 2]*\[Alpha]4 - 
     2*G*Exp[-I*\[Theta]]*\[Alpha]3) // Simplify;

\[Alpha]6 = I*Subscript[g, 2]*h7*\[Alpha]2 // Simplify;

\[Alpha]7 = -I*Subscript[g, 2]*\[Alpha]5 + 
    I*Subscript[g, 2]*\[Alpha]2*\[Alpha]3 // Simplify;

\[Alpha]8 = 
  I*Subscript[g, 2]*\[Alpha]6 + \[Alpha]2*\[Alpha]3*h7 // Simplify; 

\[Alpha]9 = SuperStar[h8]*\[Alpha]2*\[Alpha]3 // Simplify;

\[Xi]1 = (Subscript[g, 2]^2*\[Alpha]2 - \[Alpha]3*\[Alpha]2) // 
   Simplify;

\[Xi]2 = F^2*(SuperStar[h6] - h5)*\[Xi]1 - 
    h5*SuperStar[h6]*\[Alpha]9 // Simplify;

\[Alpha]10 = \[Alpha]1*\[Alpha]9 - 
    F^2*(SuperStar[h6] - h5)*\[Alpha]8 // Simplify;

\[Alpha]11 = (I*Subscript[g, 2]* h5*SuperStar[h6]*\[Alpha]9 - 
     F^2*(SuperStar[h6] - h5)*\[Alpha]7) // Simplify;

\[Xi]3 = ((\[Alpha]6*\[Xi]2 )/(\[Alpha]2*\[Alpha]3*\[Alpha]10) + (
     I*Subscript[g, 2]*\[Alpha]2)/(\[Alpha]2*\[Alpha]3)) // Simplify;

\[Alpha]12 = ((\[Alpha]5*\[Alpha]10 +
\[Alpha]6*\[Alpha]11)/(\[Alpha]10*\[Alpha]2*\[Alpha]3)) // Simplify;

\[Alpha]13 =    h1 + J^2/h3 - (I*Subscript[g, 1])/h9 - (I*Subscript[g, 2]*\[Alpha]11)/\[Alpha]10 + 2*G*Exp[-I*\[Theta]]*\[Alpha]12 // Simplify;

\[Xi]4 = ((-I*Subscript[g, 2]*\[Xi]2 )/\[Alpha]10 +    2*G*Exp[-I*\[Theta]]*\[Xi]3) // Simplify;

A1 = (\[Xi]4*Em*Exp[I*\[CapitalPhi]m] -    Sqrt[2*\[Eta]*k1]*Ep)/\[Alpha]13 // Simplify;

Tp = (Sqrt[2*\[Eta]*k1]*A1)/Ep // Simplify;

rr = Abs[(1 - ( Sqrt[2*\[Eta]*k1]*(\[Xi]4*Em*Exp[I*\[CapitalPhi]m] - 
         Sqrt[2*\[Eta]*k1]*Ep))/(Ep *\[Alpha]13))]^2 // Simplify;

v = ComplexExpand[Im[Tp]] // Simplify;

R = ComplexExpand[Re[Tp]] // Simplify;

\[Phi]12 = ArcTan[v/R] // Simplify;

pp1 = N[D[\[Phi]12, \[Delta]]] // Simplify;

\[Delta] = \[Omega]b;

Plot[Evaluate[pp1], {p1, 0.001, 0.01}, Frame -> True, 
 PlotLegends -> Placed[LineLegend[{""}, LegendLayout -> {"Column", 1}, 
    LegendMarkerSize -> {{30, 20}}], {Right, 0.90}, Pane[#, 450, Alignment -> Right] &], ImageSize -> 450, PlotStyle -> {Blue, Dashing[Large], Thickness[0.007]}, GridLines -> Automatic, 
 FrameLabel -> {Style["", 18, Bold], Style["", 18, Bold]}, 
 Evaluate@plotset, Evaluate@plotset2, Axes -> True, PlotRange -> All]

This is the link to my group delay formula.

enter image description here

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ If you try to scrape-n-paste your own code back into Mathematica you will discover a long list of problems. Some or most of this can be fixed if you put four spaces in front of each line when you are posting or perhaps when you are editing your code. This should give StackExchange the hint that this is to be displayed as code and not to have all your backslashes and some other characters ruined in your post. The reason some readers will actually want to paste your code back into MMA is so that they can actually test their idea to see if it does speed up your code enough to help you. $\endgroup$
    – Bill
    Commented Aug 3 at 4:39
  • $\begingroup$ plotset and plotset2 are nowhere defined. $\endgroup$
    – Daniel Huber
    Commented Aug 3 at 6:31
  • $\begingroup$ @ Daniel Huber could you please tell me how to simple plot of it. Because when I run it takes very long time. You just plot it and not consider the plotset and plotset2 in plot command. I just want to run if faster. $\endgroup$
    – abdul wahab
    Commented Aug 3 at 6:50
  • $\begingroup$ @ Bill I think now the code is fine. Could you please rewrite it so it can run fast. . Because when I run it takes very long time. $\endgroup$
    – abdul wahab
    Commented Aug 3 at 6:52
  • 2
    $\begingroup$ Most of the Simplify appear to be just taking up time, leaving in the very last one may help. The ComplexExpand may not be doing anything. I might try without that and see if it really changes the result or not. When a Plot is too slow the first thing I try is ListPlot[Table[pp1,{p1,0.001,0.01,0.001}],Joined->True] and if that is too slow then I try p1=0.001;pp1 and see how long that takes. Then I expect Plot will take 100 or 1000 times longer than that to plot the 100 or 1000 points. $\endgroup$
    – Bill
    Commented Aug 3 at 14:16

1 Answer 1

Reset to default
1
$\begingroup$

To speeding up the code we can remove all Simplify and use numerical differentiations instead of analytical one. This code takes less than 0.05s

Clear["Global`*"]
AbsoluteTiming[
 plotset = {FrameStyle -> Directive[Thickness[0.004]], 
   TicksStyle -> Directive[Black, 18]};
 
 plotset2 = 
  FrameTicksStyle -> {{Directive[Black, 18], 
     Directive[FontOpacity -> 0, FontSize -> 0]}, {Directive[Black, 
      18], Directive[FontOpacity -> 0, FontSize -> 0]}};
 
 G = 0.5  k1; \[Theta] = 0;
 
 \[CapitalPhi]m = 0;
 
 B = 5*10^-10;
 
 \[Gamma] = 2*\[Pi]*28*10^9;
 
 N1 = Sqrt[3*10^16];
 
 Em = Sqrt[5]/4*(\[Gamma])*N1*B;
 
 \[Omega]c = 2*\[Pi]*7.86*10^9;
 
 \[Omega]b = 2*\[Pi]*11.42*10^6;
 
 \[CapitalDelta]a = \[Omega]b;
 
 k1 = \[Pi]*3.35*10^6;
 
 k2 = k1;
 
 J = 0.5  k1;
 
 kb = 300*\[Pi];
 
 km1 = \[Pi]*1.12*10^6;
 
 km2 = \[Pi]*1.12*10^6;
 
 \[CapitalDelta]m1 = \[Omega]b;
 
 \[CapitalDelta]m2 = \[Omega]b;
 
 \[HBar] = 1.054*10^-34;
 
 Subscript[g, 1] = 2  k1;
 
 Subscript[g, 2] = 2  k1;
 
 wl = 2*\[Pi]*7.9*10^9;
 
 sss = Sqrt[p1/(\[HBar]*wl)];
 
 E1 = sss;
 
 Ep = 0.0001  E1;
 
 \[Eta] = 0.5;
 
 gmb = 2*Pi;
 
 h1 = -I*\[CapitalDelta]a + I*\[Delta] - k1;
 
 SuperStar[h1] = I*\[CapitalDelta]a - I*\[Delta] - k1;
 
 h2 = -I*\[CapitalDelta]a - I*\[Delta] - k1;
 
 SuperStar[h2] = I*\[CapitalDelta]a + I*\[Delta] - k1;
 
 h3 = -I*\[CapitalDelta]a + I*\[Delta] + k2;
 
 SuperStar[h3] = I*\[CapitalDelta]a - I*\[Delta] + k2;
 
 h4 = -I*\[CapitalDelta]a - I*\[Delta] + k2;
 
 SuperStar[h4] = I*\[CapitalDelta]a + I*\[Delta] + k2;
 
 h5 = -I*\[Omega]b + I*\[Delta] - kb;
 
 SuperStar[h5] = I*\[Omega]b - I*\[Delta] - kb;
 
 h6 = -I*\[Omega]b - I*\[Delta] - kb;
 
 SuperStar[h6] = I*\[Omega]b + I*\[Delta] - kb;
 
 h7 = -I*\[CapitalDelta]m1 + I*\[Delta] - km1;
 
 SuperStar[h7] = I*\[CapitalDelta]m1 - I*\[Delta] - km1;
 
 h8 = -I*\[CapitalDelta]m1 - I*\[Delta] - km1;
 
 SuperStar[h8] = I*\[CapitalDelta]m1 + I*\[Delta] - km1;
 
 h9 = -I*\[CapitalDelta]m2 + I*\[Delta] - km2;
 
 SuperStar[h9] = I*\[CapitalDelta]m2 - I*\[Delta] - km2;
 
 h10 = -I*\[CapitalDelta]m2 - I*\[Delta] - km2;
 
 SuperStar[h10] = I*\[CapitalDelta]m2 + I*\[Delta] - km2;
 
 \[Alpha] = ((k1 + I  \[CapitalDelta]a)*((k2 - I  \[CapitalDelta]a)) -
      J^2)/(k2 - I  \[CapitalDelta]a);
 
 \[Beta] = (k1 - I  \[CapitalDelta]a) + 
   Subscript[g, 1]^2/(-I  \[CapitalDelta]m1 + km1) + 
   Subscript[g, 2]^2/(-I  \[CapitalDelta]m2 + km2) - 
   J^2/(k2 + I  \[CapitalDelta]a);
 
 \[Gamma]1 = ((\[Alpha]*\[Beta] - 4  G^2)/\[Beta]) + (Subscript[g, 
       1]^2/(I  \[CapitalDelta]m1 + km1)) + (Subscript[g, 
       2]^2/(I  \[CapitalDelta]m2 + km2));
 \[Psi]1 = (1/\[Gamma]1)*(E1  Sqrt[
       2*k1*\[Eta]] + ((2*G*Exp[I*\[Theta]])/\[Beta]*(E1  Sqrt[
          2*k1*\[Eta]])));
 
 \[Psi]2 = (-I*Subscript[g, 2])/(km1 + I  \[CapitalDelta]m1);
 
 ms = \[Psi]1*\[Psi]2;
 
 F = ms*gmb;
 
 \[Alpha]1 = h7*SuperStar[h6]*h5 + F^2*(SuperStar[h6] - h5);
 
 \[Alpha]2 = 
  SuperStar[h2] + J^2/SuperStar[h4] + Subscript[g, 1]^2/SuperStar[h10];
 
 \[Alpha]3 = \[Alpha]2*SuperStar[h8] + Subscript[g, 2]^2;
 
 \[Alpha]4 = 
  2*G*Exp[-I*\[Theta]]*I*Subscript[g, 2] + I*Subscript[g, 2]*\[Alpha]2;
 
 \[Alpha]5 = (-I*Subscript[g, 2]*\[Alpha]4 - 
    2*G*Exp[-I*\[Theta]]*\[Alpha]3);
 
 \[Alpha]6 = I*Subscript[g, 2]*h7*\[Alpha]2;
 
 \[Alpha]7 = -I*Subscript[g, 2]*\[Alpha]5 + 
   I*Subscript[g, 2]*\[Alpha]2*\[Alpha]3;
 
 \[Alpha]8 = I*Subscript[g, 2]*\[Alpha]6 + \[Alpha]2*\[Alpha]3*h7;
 
 \[Alpha]9 = SuperStar[h8]*\[Alpha]2*\[Alpha]3;
 
 \[Xi]1 = (Subscript[g, 2]^2*\[Alpha]2 - \[Alpha]3*\[Alpha]2);
 
 \[Xi]2 = F^2*(SuperStar[h6] - h5)*\[Xi]1 - h5*SuperStar[h6]*\[Alpha]9;
 
 \[Alpha]10 = \[Alpha]1*\[Alpha]9 - F^2*(SuperStar[h6] - h5)*\[Alpha]8;
 
 \[Alpha]11 = (I*Subscript[g, 2]*h5*SuperStar[h6]*\[Alpha]9 - 
    F^2*(SuperStar[h6] - h5)*\[Alpha]7);
 
 \[Xi]3 = ((\[Alpha]6*\[Xi]2)/(\[Alpha]2*\[Alpha]3*\[Alpha]10) + (I*
       Subscript[g, 2]*\[Alpha]2)/(\[Alpha]2*\[Alpha]3));
 
 \[Alpha]12 = ((\[Alpha]5*\[Alpha]10 + \
\[Alpha]6*\[Alpha]11)/(\[Alpha]10*\[Alpha]2*\[Alpha]3));
 
 \[Alpha]13 = 
  h1 + J^2/h3 - (I*Subscript[g, 1])/
    h9 - (I*Subscript[g, 2]*\[Alpha]11)/\[Alpha]10 + 
   2*G*Exp[-I*\[Theta]]*\[Alpha]12;
 
 \[Xi]4 = ((-I*Subscript[g, 2]*\[Xi]2)/\[Alpha]10 + 
    2*G*Exp[-I*\[Theta]]*\[Xi]3);
 
 A1 = (\[Xi]4*Em*Exp[I*\[CapitalPhi]m] - 
     Sqrt[2*\[Eta]*k1]*Ep)/\[Alpha]13;
 
 Tp = (Sqrt[2*\[Eta]*k1]*A1)/Ep;
 
 rr = Abs[(1 - (Sqrt[
          2*\[Eta]*k1]*(\[Xi]4*Em*Exp[I*\[CapitalPhi]m] - 
           Sqrt[2*\[Eta]*k1]*Ep))/(Ep*\[Alpha]13))]^2;
 
 (*v=ComplexExpand[Im[Tp]];
 
 R=ComplexExpand[Re[Tp]];*)
 
 
 \[Phi]12 = Arg[Tp];
 
 del = 0.1; dp1 = Range[.001, .01, .0001]; 
 dpp1 = Table[
   Evaluate[\[Phi]12], {\[Delta], \[Omega]b, \[Omega]b + del, 
    del}, {p1, dp1}]; 
 pp1 = Transpose[{dp1, (dpp1[[2]] - dpp1[[1]])/del}];]

Visualization

ListLinePlot[pp1, Frame -> True, 
 PlotLegends -> 
  Placed[LineLegend[{""}, LegendLayout -> {"Column", 1}, 
    LegendMarkerSize -> {{30, 20}}], {Right, 0.90}, 
   Pane[#, 450, Alignment -> Right] &], ImageSize -> 450, 
 PlotStyle -> {Blue, Dashing[Large], Thickness[0.007]}, 
 GridLines -> Automatic, 
 FrameLabel -> {Style["", 18, Bold], Style["", 18, Bold]}, 
 Evaluate@plotset, Evaluate@plotset2, Axes -> True, PlotRange -> All]

Figure 1

Improve this answer
$\endgroup$
2
  • $\begingroup$ thank you so much. Yes, it works very fast. $\endgroup$
    – abdul wahab
    Commented Aug 6 at 2:27
  • 1
    $\begingroup$ @abdulwahab You are welcome! $\endgroup$
    – Alex Trounev
    Commented Aug 6 at 2:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.