1
$\begingroup$

I'm trying to calculate the $A$ and $F$ values by using this code,

       AA[z_, B_] = -a z^2 - d B^2 z^5;
    
    G1dd[zh_, 
       B_] = {{-((E^(2 AA[z, B]) L^2 g[z, zh, B])/z^2), 0, 0, 0, 
        0}, {0, (E^(2 AA[z, B]) L^2)/(z^2 g[z, zh, B]), 0, 0, 0}, {0, 
        0, (E^(2 AA[z, B]) L^2)/z^2, 0, 0}, {0, 0, 
        0, (E^(B^2 z^2 + 2 AA[z, B]) L^2)/z^2, 2 B π αp}, {0, 
        0, 0, -2 B π αp, (E^(B^2 z^2 + 2 AA[z, B]) L^2)/z^2}};
    G1uu[zh_, 
       B_] = {{-((E^(-2 AA[z, B]) z^2)/(L^2 g[z, zh, B])), 0, 0, 0, 
        0}, {0, (E^(-2 AA[z, B]) z^2 g[z, zh, B])/L^2, 0, 0, 0}, {0, 
        0, (E^(-2 AA[z, B]) z^2)/L^2, 0, 0}, {0, 0, 
        0, (E^(B^2 z^2 + 
              2 AA[z, B]) L^2 z^2)/(E^(2 B^2 z^2 + 4 AA[z, B]) L^4 + 
           4 B^2 π^2 z^4 αp^2), -((2 B π z^4 \
    αp)/(E^(2 B^2 z^2 + 4 AA[z, B]) L^4 + 
             4 B^2 π^2 z^4 αp^2))}, {0, 0, 
        0, (2 B π z^4 αp)/(E^(2 B^2 z^2 + 4 AA[z, B]) L^4 + 
           4 B^2 π^2 z^4 αp^2), (E^(B^2 z^2 + 
              2 AA[z, B]) L^2 z^2)/(E^(2 B^2 z^2 + 4 AA[z, B]) L^4 + 
           4 B^2 π^2 z^4 αp^2)}};
    g[z_?NumericQ, zh_?NumericQ, B_?NumericQ] := 
     1 - NIntegrate[ζ^3 Exp[-B^2 ζ^2 - 
           3 AA[ζ, B]], {ζ, 0, z}]/
       NIntegrate[ζ^3 Exp[-B^2 ζ^2 - 
           3 AA[ζ, B]], {ζ, 0, zh}]
    Temp[zh_?NumericQ, 
      B_?NumericQ] := (E^(-B^2 zh^2 - 
           3 AA[zh, B]) zh^3)/(4 π NIntegrate[
         E^(-B^2 ζ^2 - 3 AA[ζ, B]) ζ^3, {ζ, 0, 
          zh}])
    
    ff[zh_, B_] = -(1/z^2) 2 E^(2 (B^2 z^2 + AA[z, B])) L^2 Sqrt[
        1 + (B^2 E^(-2 B^2 z^2 - 4 AA[z, B]) z^4)/(bb^2 L^4)] (g[z, zh, 
           B] (-2 + 2 B^2 z^2 + 3 z D[AA[z, B], z]) + z D[g[z, zh, B], z]);
    term[z_, zh_, B_] = 
      Sqrt[-Det[G1dd[zh, B]]] G1uu[zh, B][[2, 2]] G1uu[zh, B][[3, 3]] ff[
        zh, B];
    αvalue[zh_, B_, ω_] = (I ω)/(4 π Temp[zh, B]);
    ϵ = 0.0001; L = 1; ls = 2.702; αp = ls^2; bb = 
     1/(2 π αp); a = 0.15; d = 0.013;
    zhBpt2 = 1.342;
    sol[zh_, B_, ω_] := 
     NDSolve[{f''[z] + 
         D[Log[term[z, zhBpt2, 0.2]], z] f'[
           z] - ω^2 ((G1uu[zhBpt2, 0.2][[1, 1]]/
             G1uu[zhBpt2, 0.2][[2, 2]])) f[z] == 0, 
       f[ϵ] == 
        1 + ω^2/
          2 ϵ^2 - (1/2 0.2^2 ω^2 + ω^4/
            16) Log[ϵ] ϵ^4, 
       f'[ϵ] == ω^2/
          2 ϵ - (1/2 0.2^2 ω^2 + ω^4/
            16) ϵ^3 - (2 0.2^2 ω^2 + ω^4/
            4) Log[ϵ] ϵ^3}, 
      f, {z, ϵ, zhBpt2 - ϵ}]
    v1[zh_, B_, ω_] := 
     Part[Simplify[(A + 
             I F) (1 - (zhBpt2 - 0.0001)/zhBpt2)^-αvalue[zhBpt2, 
              0.2, ω] + (A - 
             I F) (1 - (zhBpt2 - 0.0001)/zhBpt2)^αvalue[zhBpt2, 
             0.2, ω] == Evaluate[f[zhBpt2 - 0.0001]] /. 
        sol[zh, B, ω]], 1]
    v2[zh_, B_, ω_] := 
     Part[Simplify[(A + 
             I F) (1 - (zhBpt2 - 0.0002)/zhBpt2)^-αvalue[zhBpt2, 
              0.2, ω] + (A - 
             I F) (1 - (zhBpt2 - 0.0002)/zhBpt2)^αvalue[zhBpt2, 
             0.2, ω] == Evaluate[f[zhBpt2 - 0.0002]] /. 
        sol[zh, B, ω]], 1]
    s[zh_, B_, ω_] := {A, F} /. 
      Solve[v1[zh, B, ω] && v2[zh, B, ω], {A, F}]
    AnFvalues = 
     Flatten[Table[s[zhBpt2, 0.2, ω], {ω, 0.01, 3, 0.01}], 
      1]

It has been running for the last 5 hours without any result till now.

In the code, when I defined AA[z] as,

 AA[z_] = -a z^2;

and used Integrate in place of NIntegrate in g[z, zh, B] and in Temp[zh, B] as shown in the code below,

        AA[z_] = -a z^2;
    
    G1dd[zh_, 
       B_] = {{-((E^(2 AA[z]) L^2 g[z, zh, B])/z^2), 0, 0, 0, 
        0}, {0, (E^(2 AA[z]) L^2)/(z^2 g[z, zh, B]), 0, 0, 0}, {0, 
        0, (E^(2 AA[z]) L^2)/z^2, 0, 0}, {0, 0, 
        0, (E^(B^2 z^2 + 2 AA[z]) L^2)/z^2, 2 B π αp}, {0, 0, 
        0, -2 B π αp, (E^(B^2 z^2 + 2 AA[z]) L^2)/z^2}};
    G1uu[zh_, 
       B_] = {{-((E^(-2 AA[z]) z^2)/(L^2 g[z, zh, B])), 0, 0, 0, 
        0}, {0, (E^(-2 AA[z]) z^2 g[z, zh, B])/L^2, 0, 0, 0}, {0, 
        0, (E^(-2 AA[z]) z^2)/L^2, 0, 0}, {0, 0, 
        0, (E^(B^2 z^2 + 2 AA[z]) L^2 z^2)/(E^(2 B^2 z^2 + 4 AA[z]) L^4 + 
           4 B^2 π^2 z^4 αp^2), -((2 B π z^4 \
    αp)/(E^(2 B^2 z^2 + 4 AA[z]) L^4 + 
             4 B^2 π^2 z^4 αp^2))}, {0, 0, 
        0, (2 B π z^4 αp)/(E^(2 B^2 z^2 + 4 AA[z]) L^4 + 
           4 B^2 π^2 z^4 αp^2), (E^(B^2 z^2 + 
              2 AA[z]) L^2 z^2)/(E^(2 B^2 z^2 + 4 AA[z]) L^4 + 
           4 B^2 π^2 z^4 αp^2)}};
    
    g[z_, zh_, B_] := 
     1 - Integrate[ζ^3 Exp[-B^2 ζ^2 - 
           3 AA[ζ]], {ζ, 0, z}]/
       Integrate[ζ^3 Exp[-B^2 ζ^2 - 3 AA[ζ]], {ζ, 
         0, zh}]
    
    Temp[zh_, 
      B_] := (E^(-B^2 zh^2 - 3 AA[zh]) zh^3)/(4 π Integrate[
         E^(-B^2 ζ^2 - 3 AA[ζ]) ζ^3, {ζ, 0, zh}])
    
    
    ff[zh_, B_] = -(1/z^2) 2 E^(2 (B^2 z^2 + AA[z])) L^2 Sqrt[
        1 + (B^2 E^(-2 B^2 z^2 - 4 AA[z]) z^4)/(bb^2 L^4)] (g[z, zh, 
           B] (-2 + 2 B^2 z^2 + 3 z D[AA[z], z]) + z D[g[z, zh, B], z]);
    
    term[z_, zh_, B_] = 
      Sqrt[-Det[G1dd[zh, B]]] G1uu[zh, B][[2, 2]] G1uu[zh, B][[3, 3]] ff[
        zh, B];
    
    αvalue[zh_, B_, ω_] = (I ω)/(4 π Temp[zh, B]);
    
    ϵ = 0.0001; L = 1; ls = 2.702; αp = ls^2; bb = 
     1/(2 π αp); a = 0.15; d = 0.013;
    zhBpt2 = 1.342;
    
    sol[zh_, B_, ω_] := 
     NDSolve[{f''[z] + 
         D[Log[term[z, zhBpt2, 0.2]], z] f'[
           z] - ω^2 ((G1uu[zhBpt2, 0.2][[1, 1]]/
             G1uu[zhBpt2, 0.2][[2, 2]])) f[z] == 0, 
       f[ϵ] == 
        1 + ω^2/
          2 ϵ^2 - (1/2 0.2^2 ω^2 + ω^4/
            16) Log[ϵ] ϵ^4, 
       f'[ϵ] == ω^2/
          2 ϵ - (1/2 0.2^2 ω^2 + ω^4/
            16) ϵ^3 - (2 0.2^2 ω^2 + ω^4/
            4) Log[ϵ] ϵ^3}, 
      f, {z, ϵ, zhBpt2 - ϵ}]
    
    v1[zh_, B_, ω_] := 
     Part[Simplify[(A + 
             I F) (1 - (zhBpt2 - 0.0001)/zhBpt2)^-αvalue[zhBpt2, 
              0.2, ω] + (A - 
             I F) (1 - (zhBpt2 - 0.0001)/zhBpt2)^αvalue[zhBpt2, 
             0.2, ω] == Evaluate[f[zhBpt2 - 0.0001]] /. 
        sol[zh, B, ω]], 1]
    
    v2[zh_, B_, ω_] := 
     Part[Simplify[(A + 
             I F) (1 - (zhBpt2 - 0.0002)/zhBpt2)^-αvalue[zhBpt2, 
              0.2, ω] + (A - 
             I F) (1 - (zhBpt2 - 0.0002)/zhBpt2)^αvalue[zhBpt2, 
             0.2, ω] == Evaluate[f[zhBpt2 - 0.0002]] /. 
        sol[zh, B, ω]], 1]
    
    s[zh_, B_, ω_] := {A, F} /. 
      Solve[v1[zh, B, ω] && v2[zh, B, ω], {A, F}]
    
    AnFvalues = 
     Flatten[Table[s[zhBpt2, 0, ω], {ω, 0.01, 3, 0.01}], 1]

it took around 30 seconds to compute the values, but with the present form of AA[z] as,

AA[z_, B_] = -a z^2 - d B^2 z^5;

and using NIntegrate, it keeps running.

Can anyone explain why it takes so much computation time and how to overcome this?

$\endgroup$
6
  • $\begingroup$ Why do you want to use NIntegerate if Integrate works? $\endgroup$ Commented Apr 17 at 15:45
  • $\begingroup$ When I'm using the second AA[z] i.e. $- a z^2 -d B^2 z^5$, Integrate is not working for finite values of $B$, that's why I'm using NIntegrate. For first AA[z] Integrate is working for finite $B$. $\endgroup$ Commented Apr 18 at 3:35
  • $\begingroup$ Please see the edited codes $\endgroup$ Commented Apr 18 at 4:32
  • $\begingroup$ If you make the all dependencies explicit, for example g[a_?NumericQ, d_?NumericQ, z_?NumericQ, zh_?NumericQ, B_?NumericQ] := ..., then it becomes easier to test the code step-by-step and find out where the bottleneck is, and start optimizing the steps. As it is, your code is too messy to analyze. $\endgroup$
    – Roman
    Commented Apr 18 at 7:18
  • $\begingroup$ Derivatives like D[g[z, zh, B], z] should probably be defined as separate functions: dgdz[a_?NumericQ, d_?NumericQ, z_?NumericQ, zh_?NumericQ, B_?NumericQ] := ... with an explicit formulation of the derivative (remember the fundamental theorem of calculus). $\endgroup$
    – Roman
    Commented Apr 18 at 7:20

0

Your Answer

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